SeminarStuff/SnakeFile at main · MedaFour/SeminarStuff · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
configfile: "config.yaml"
SAMPLES = config["edirect"]["samples"]  # Load the sample dictionary

rule all:
    input:
       expand("results/{sample}_mut_mspe.snp", sample=SAMPLES.keys()),

rule fetch_data:
    output:
       'data/{sample}.fna'
    params:
        id=lambda wildcards: SAMPLES[wildcards.sample]
    conda:
       'envs/edirect.yaml'
    shell:
       'efetch -db nuccore -format {config[edirect][format]} '
       '-id {params.id} |'
       'head -n {config[edirect][subset_size]} > {output}'

rule mutate_sequence:
    input:
       rules.fetch_data.output
    output:
       'results/{sample}_mut.fna'
    conda:
       'envs/emboss.yaml'
    shell:
       'msbar -sequence {input} -outseq {output} '
       '-codon {config[emboss][codon]} '
       '-count {config[emboss][count]} '
       '-point {config[emboss][point]} '
       '-block {config[emboss][block]} '

rule simulate_reads:
    input:
       rules.mutate_sequence.output
    output:
       'results/{sample}_mut_mspe1.fq',
       'results/{sample}_mut_mspe2.fq'
    conda:
       'envs/art.yaml'
    shell:
       'art_illumina -i {input} -ss {config[art][seq_sys]} '
       '-l {config[art][read_length]} '
       '-f {config[art][coverage]} '
       '-m {config[art][fragment_mean]} '
       '-s {config[art][fragment_std]} '
       '-o results/{wildcards.sample}_mut_mspe -na'

rule mapping_reads:
    input:
       read0 = rules.simulate_reads.output[0],
       read1 = rules.simulate_reads.output[1],
       seq = rules.fetch_data.output
    output:
       'results/{sample}_mut_mspe.sam'
    conda:
       'envs/bwa.yaml'
    shell:
       '''
       bwa index {input.seq} &&
       bwa mem {input.seq} {input.read0} {input.read1} >{output}
       '''
rule sam_to_bam:
    input:
       rules.mapping_reads.output
    output:
       sortedbam ='results/{sample}_mut_mspe.sortedbam',
       bam ='results/{sample}_mut_mspe.bam'
    conda:
       'envs/samtools.yaml'
    shell:
       '''
       samtools view -b {input} -o {output.bam} &&
       samtools sort -o {output.sortedbam} {output.bam}
       '''
rule pileup:
    input:
       sortedbam = rules.sam_to_bam.output.sortedbam,
       seq = rules.fetch_data.output
    output:
       'results/{sample}_mut_mspe.mpileup'
    conda:
       'envs/samtools.yaml'
    shell:
       '''
       samtools index {input.sortedbam} &&
       samtools mpileup -f {input.seq} -C {config[pileup][adjustment]} -o {output} {input.sortedbam}
       '''

rule pileup_to_snp:
    input:
       rules.pileup.output
    output:
       'results/{sample}_mut_mspe.snp'
    conda:
       'envs/varscan.yaml'
    shell:
       'varscan mpileup2snp {input} > {output}'