-
Notifications
You must be signed in to change notification settings - Fork 0
/
snp_mask.smk
106 lines (94 loc) · 3.04 KB
/
snp_mask.smk
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
100
101
102
103
104
105
106
import os
configfile: "config/config.yaml"
cluster_config: "config/cluster.yaml"
include: "helpers.py"
#make sure the output folder exists before running anything
hisat_outdir = get_output_dir(config["project_top_level"], config['histat3n_output_folder'])
os.system("mkdir -p {0}".format(hisat_outdir))
SAMPLES = pd.read_csv(config["sampleCSVpath"], sep = ",")
SAMPLES = SAMPLES.replace(np.nan, '', regex=True)
SAMPLE_NAMES = SAMPLES['sample_name'].tolist()
GENOME_DIR = config['GENOME_DIR']
GENOME_FA = config['fasta']
chr_list = ['chr1', 'chr2', 'chr3', 'chr4', 'chr5',
'chr6', 'chr7', 'chr8', 'chr9', 'chr10', 'chr11', 'chr12',
'chr13', 'chr14', 'chr15', 'chr16', 'chr17', 'chr18', 'chr19',
'chr20', 'chr21', 'chr22', 'chrX', 'chrY', 'chrM']
rule all_call:
input:
expand(hisat_outdir + "{name}.snpmasked.bam", name = SAMPLE_NAMES),
expand(hisat_outdir + "{name}.snpmasked.bam.bai", name = SAMPLE_NAMES)
rule call_snp_mask:
wildcard_constraints:
name="|".join(SAMPLE_NAMES)
input:
hisat_outdir + "{name}.sorted.bam"
output:
temp(expand(hisat_outdir + "{{name}}.sorted{chr}.snpmasked.bam",chr=chr_list))
params:
vcf = config['vcf_path'] #TODO this could be altered to take a VCF per sample in the future
threads:
4
shell:
"""
echo "Processing {input} file..."
python3 scripts/snp_mask.py --bam {input} --vcf {params.vcf} --cpu {threads}
"""
rule merge_chromosomes:
wildcard_constraints:
name="|".join(SAMPLE_NAMES),
chr="|".join(chr_list)
input:
expand(hisat_outdir + "{{name}}.sorted{chr}.snpmasked.bam",chr=chr_list)
output:
hisat_outdir + "{name}.snpmasked.bam"
shell:
"""
samtools merge {output} {input}
"""
rule index_merged:
wildcard_constraints:
name="|".join(SAMPLE_NAMES),
input:
hisat_outdir + "{name}.snpmasked.bam"
output:
hisat_outdir + "{name}.snpmasked.bam.bai"
shell:
"""
samtools index {input}
"""
# rule snp_mask_bams:
# input:
# bam = hisat_outdir + "{name}.sorted.bam",
# bai = hisat_outdir + "{name}.sorted.bam.bai"
# output:
# hisat_outdir + "{name}.sorted.tagged.bam"
# params:
# pickled = GENOME_FA + '.pickle'
# shell:
# """
# python3 scripts/slamdunk_taggers.py -b {input.bam} -p {params.pickled}
# """
# rule tag_bams:
# input:
# bam = hisat_outdir + "{name}.sorted.bam",
# bai = hisat_outdir + "{name}.sorted.bam.bai"
# output:
# hisat_outdir + "{name}.sorted.tagged.bam"
# params:
# pickled = GENOME_FA + '.pickle'
# shell:
# """
# python3 scripts/slamdunk_taggers.py -b {input.bam} -p {params.pickled}
# """
# rule index_tagged_bams:
# input:
# hisat_outdir + "{name}.sorted.tagged.bam"
# output:
# hisat_outdir + "{name}.sorted.tagged.bam.bai"
# threads:
# 4
# shell:
# """
# samtools index {input}
# """