makefiles/mpileup.mk

   1 OUTMAIN := mpileup
   2 DATESTR := $(shell date +%Y%m%d%Z%H%M%S)
   3 OUTP := $(OUTMAIN)_$(DATESTR)
   4 OUT := $(OUTP).bcf
   5 VCF := $(OUTP).vcf
   6 SAMTOOLS := samtools
   7 BCFTOOLS := bcftools
   8
   9 REF := /bak/seqdata/genomes/Felis_catus_80_masked/Felis_catus80_chr.fa
  10 CMD := $(SAMTOOLS) mpileup -g -d 1000 -t AD,ADF,DP,DPR,DV,DP4,SP -f $(REF)
  11
  12 #BAMS := pti096_clean_aln_pe_rmdup.bam pti183_clean_aln_pe_rmdup.bam pti301_clean_aln_pe_rmdup.bam pti332_clean_aln_pe_rmdup.bam
  13 BAMS := $(sort $(wildcard *.bam))
  14 #CHRS := scaffold1000 scaffold979 scaffold982
  15 RAWCHRS := $(sort $(shell $(SAMTOOLS) view -H $(firstword $(BAMS)) | sed -n 's/^@SQ\tSN:\([^\t]*\).*/\1/p'))
  16
  17 p+ = $(subst |,+,$1)
  18 +p = $(subst +,|,$1)
  19
  20 CHRS := $(call p+,$(RAWCHRS))
  21
  22 BYCHR = $(addsuffix .bcf,$(addprefix bychr/_,$(CHRS)))
  23 #BYCHRINX = $(addsuffix .csi,$(BYCHR))
  24 BAIS = $(addsuffix .bai,$(BAMS))
  25 #GATKBAIS = $(BAMS:.bam=.bai)
  26
  27 define newline
  28
  29 endef
  30
  31 NOOP :=
  32 space = $(NOOP) $(NOOP)
  33
  34 .PHONY: all clean cleanall list
  35
  36 all: list $(OUT) bcfbychr.lst $(VCF)
  37         $(BCFTOOLS) index $(VCF).gz &
  38         $(BCFTOOLS) stats -F $(REF) -s - -d 0,1500,1 $(VCF).gz >$(VCF).stats
  39         $(BCFTOOLS) norm -Df $(REF) -c e -m+both $(VCF).gz | $(BCFTOOLS) filter -sLowQual -e'%QUAL<10' | $(BCFTOOLS) filter -m+ -sDepthHigh -e'DP>650' | $(BCFTOOLS) filter -m+ -sDepthLow -e'DP<2' | $(BCFTOOLS) filter -m+ -sBadSites -e'%QUAL<10 && RPB<0.1' | tee >($(BCFTOOLS) view -Oz -o $(VCF).filtering.gz) | $(BCFTOOLS) view -f .,PASS -Oz -o $(VCF).filtered.gz
  40         $(BCFTOOLS) index $(VCF).filtered.gz &
  41         $(BCFTOOLS) index $(VCF).filtering.gz
  42
  43 list:
  44         @echo -e '$(subst $(space),\n,$(BAMS))' >_bams.lst
  45         @echo -e '$(subst $(space),\n,$(CHRS))' >_chrs.lst
  46         @echo -e '\n$(subst $(space),\n,$(RAWCHRS))' >>_chrs.lst
  47         @echo -e '$(subst $(space),\n,$(BYCHR))' >_bcfbychr.lst
  48
  49 bai: $(BAIS)
  50
  51 $(OUT): $(BYCHR) bcfbychr.lst
  52         ulimit -n 2048
  53         $(BCFTOOLS) concat -a -O b -f bcfbychr.lst -o $(OUT)
  54
  55 $(VCF): $(OUT)
  56         $(BCFTOOLS) call -vmO z -f GQ,GP -o $(@).gz $<
  57
  58 bychr/:
  59         mkdir bychr
  60
  61 bcfbychr.lst: $(BYCHR)
  62         @echo dbg2 [$@] [$*] [$^] [$|]
  63         echo -e '$(subst $(space),\n,$(BYCHR))' >bcfbychr.lst
  64
  65 %.bai: $*
  66         $(SAMTOOLS) index $*
  67
  68 bychr/_%.bcf: $(BAMS) bai | bychr/
  69         $(CMD) -r "$(call +p,$(*))" $(BAMS) >$@
  70         $(BCFTOOLS) index $@
  71
  72 clean:
  73         -rm -fr bychr
  74
  75 cleanall:
  76         -rm -fr bychr
  77         -rm bcfbychr.lst
  78         -rm $(OUT)