|
| 1 | +configfile: "config/config.yaml" |
| 2 | + |
| 3 | + |
| 4 | +rule all: |
| 5 | + input: |
| 6 | + "calls/all.vcf", |
| 7 | + "calls/positions.png", |
| 8 | + "calls/quals.svg", |
| 9 | + |
| 10 | + |
| 11 | +def get_bwa_map_input_fastqs(wildcards): |
| 12 | + return config["samples"][wildcards.sample] |
| 13 | + |
| 14 | + |
| 15 | +rule bwa_mem: |
| 16 | + input: |
| 17 | + reads=get_bwa_map_input_fastqs, |
| 18 | + idx=multiext("data/genome.fa", ".amb", ".ann", ".pac", ".sa"), |
| 19 | + output: |
| 20 | + "sorted_reads/{sample}.bam", |
| 21 | + log: |
| 22 | + "logs/bwa_mem/{sample}.log", |
| 23 | + params: |
| 24 | + extra=r"-R '@RG\tID:{sample}\tSM:{sample}'", |
| 25 | + sorting="samtools", # Can be 'none', 'samtools' or 'picard'. |
| 26 | + sort_order="coordinate", # Can be 'queryname' or 'coordinate'. |
| 27 | + sort_extra="", # Extra args for samtools/picard. |
| 28 | + threads: 8 |
| 29 | + wrapper: |
| 30 | + "v5.7.0/bio/bwa/mem" |
| 31 | + |
| 32 | + |
| 33 | +rule samtools_index: |
| 34 | + input: |
| 35 | + "sorted_reads/{sample}.bam", |
| 36 | + output: |
| 37 | + "sorted_reads/{sample}.bam.bai", |
| 38 | + log: |
| 39 | + "logs/samtools_index/{sample}.log", |
| 40 | + params: |
| 41 | + extra="", # optional params string |
| 42 | + threads: 4 # This value - 1 will be sent to -@ |
| 43 | + wrapper: |
| 44 | + "v5.7.0/bio/samtools/index" |
| 45 | + |
| 46 | +rule bcftools_call: |
| 47 | + input: |
| 48 | + fa="data/genome.fa", |
| 49 | + bam=expand("sorted_reads/{sample}.bam", sample=config["samples"]), |
| 50 | + bai=expand("sorted_reads/{sample}.bam.bai", sample=config["samples"]), |
| 51 | + output: |
| 52 | + "calls/all.vcf", |
| 53 | + log: |
| 54 | + "logs/bcftools_call/bcftools.log", |
| 55 | + shell: |
| 56 | + "(bcftools mpileup -f {input.fa} {input.bam} | " |
| 57 | + "bcftools call -mv - > {output}) 2> {log}" |
| 58 | + |
| 59 | + |
| 60 | +rule plot_positions: |
| 61 | + input: |
| 62 | + rules.bcftools_call.output, |
| 63 | + output: |
| 64 | + "calls/positions.png", |
| 65 | + script: |
| 66 | + "scripts/plot-positions.py" |
| 67 | + |
| 68 | + |
| 69 | +rule plot_quals: |
| 70 | + input: |
| 71 | + "calls/all.vcf", |
| 72 | + output: |
| 73 | + "calls/quals.svg", |
| 74 | + script: |
| 75 | + "scripts/plot-quals.py" |
0 commit comments