3
|
1 import glob
|
|
2 import os
|
|
3 import shutil
|
|
4
|
|
5 import yaml
|
|
6 configfile: "config.yaml"
|
|
7
|
|
8
|
|
9 SAMPLES = []
|
|
10 with open("config.yaml", "r") as yaml_file:
|
|
11 genome_data = yaml.safe_load(yaml_file)
|
|
12 if "ids" in genome_data.keys():
|
|
13 for id in genome_data["ids"]:
|
|
14 SAMPLES.append(id)
|
|
15 if "input_genbanks" in genome_data.keys():
|
|
16 for gb_path in genome_data["input_genbanks"]:
|
|
17 cmd = "grep 'ACCESSION' "+gb_path
|
|
18 returned_value = subprocess.getoutput(cmd)
|
|
19 words = returned_value.split()
|
|
20 SAMPLES.append(words[1])
|
|
21
|
|
22
|
|
23
|
|
24 rule final:
|
|
25 input:
|
|
26 "outputs/GCskew.txt",
|
|
27 "outputs/cactus_outdir/output_cactus.sv.gfa.gz",
|
|
28 "outputs/cactus_outdir/full.gfa",
|
|
29 "outputs/cactus_outdir/full.gfa.png",
|
|
30 "outputs/cactus_outdir/output_cactus.sv.gfa.pangenome.fasta",
|
|
31 "outputs/pav_matrix.tsv",
|
|
32 "outputs/heatmap.svg.gz",
|
|
33 "outputs/cog_output.txt",
|
|
34 "outputs/rarefaction_curves.txt"
|
|
35
|
|
36
|
|
37 rule ncbi_datasets:
|
|
38 input:
|
|
39 "config.yaml"
|
|
40 output:
|
|
41 expand("outputs/genomes/{sample}.fasta", sample=SAMPLES),
|
|
42 expand("outputs/genomes/{sample}.gb", sample=SAMPLES),
|
|
43 expand("outputs/genomes/{sample}.prt", sample=SAMPLES),
|
|
44 expand("outputs/genomes/{sample}.nuc", sample=SAMPLES),
|
|
45 genomes="outputs/genomes/genomes.txt",
|
|
46 strains="outputs/genomes/strains.txt"
|
|
47 shell:
|
|
48 """
|
|
49 perl $PANEX_PATH/Perl/get_data.pl {input} outputs/genomes
|
|
50 """
|
|
51
|
|
52
|
|
53 rule gcskew:
|
|
54 input:
|
|
55 "outputs/genomes/{sample}.fasta"
|
|
56 output:
|
|
57 "outputs/genomes/{sample}.fasta.gcskew.txt"
|
|
58 shell:
|
|
59 """
|
|
60 python3 $PANEX_PATH/SkewIT/src/gcskew.py -i {input} -o {input}.gcskew.txt -k 1000 -w 1000
|
|
61 """
|
|
62
|
|
63 rule concat_gcskew:
|
|
64 input:
|
|
65 expand("outputs/genomes/{sample}.fasta.gcskew.txt", sample=SAMPLES)
|
|
66 output:
|
|
67 out2="outputs/GCskew.txt"
|
|
68 shell:
|
|
69 """
|
|
70 cat {input} >>{output.out2}
|
|
71 """
|
|
72
|
|
73 rule genbank2gff3:
|
|
74 input:
|
|
75 "outputs/genomes/{sample}.gb"
|
|
76 output:
|
|
77 gff1="outputs/genomes/{sample}.gb.gff",
|
|
78 shell:
|
|
79 """
|
|
80 perl $PANEX_PATH/Perl/bp_genbank2gff3.pl -o outputs/genomes {input}
|
|
81 """
|
|
82
|
|
83 rule cactus:
|
|
84 input:
|
|
85 expand("outputs/genomes/{sample}.fasta", sample=SAMPLES),
|
|
86 expand("outputs/genomes/{sample}.gb.gff", sample=SAMPLES),
|
|
87 params:
|
|
88 reference=SAMPLES[0]
|
|
89 output:
|
|
90 gfasv="outputs/cactus_outdir/output_cactus.sv.gfa.gz",
|
|
91 gfafull="outputs/cactus_outdir/full.gfa"
|
|
92 shell:
|
|
93 """
|
|
94 mkdir outputs/cactus_work
|
|
95 cactus-pangenome outputs/jobstore outputs/genomes/seqfile --outDir outputs/cactus_outdir --outName output_cactus --reference {params.reference} --workDir outputs/cactus_work --vcf --mapCores 12 --gfa clip full --viz
|
|
96 zcat outputs/cactus_outdir/output_cactus.full.gfa.gz >{output.gfafull}
|
|
97 """
|
|
98
|
|
99
|
|
100
|
|
101
|
|
102 rule gfatools:
|
|
103 input:
|
|
104 "outputs/cactus_outdir/output_cactus.sv.gfa.gz"
|
|
105 output:
|
|
106 "outputs/cactus_outdir/output_cactus.sv.gfa.pangenome.fasta"
|
|
107 shell:
|
|
108 """
|
|
109 gfatools gfa2fa -s {input} > {output}
|
|
110 """
|
|
111
|
|
112 rule minigraph:
|
|
113 input:
|
|
114 cds_fasta="outputs/genomes/{sample}.nuc",
|
|
115 gfa="outputs/cactus_outdir/output_cactus.sv.gfa.gz",
|
|
116 output:
|
|
117 gaf="outputs/genomes/{sample}.nuc.pangenome.gaf",
|
|
118 shell:
|
|
119 """
|
|
120 minigraph -x lr {input.gfa} {input.cds_fasta} > {output.gaf}
|
|
121 """
|
|
122
|
|
123 rule odgi:
|
|
124 input:
|
|
125 "outputs/cactus_outdir/full.gfa"
|
|
126 output:
|
|
127 og="outputs/cactus_outdir/full.gfa.og",
|
|
128 png="outputs/cactus_outdir/full.gfa.png",
|
|
129 shell:
|
|
130 """
|
|
131 odgi build -g {input} -o {output.og}
|
|
132 odgi viz -i {output.og} -o {output.png}
|
|
133 """
|
|
134
|
|
135 rule gaf2bed:
|
|
136 input:
|
|
137 "outputs/genomes/{sample}.nuc.pangenome.gaf"
|
|
138 output:
|
|
139 "outputs/genomes/{sample}.nuc.pangenome.gaf.bed",
|
|
140 shell:
|
|
141 """
|
|
142 awk '{{print $6"\\t"$8"\\t"$9"\\t"$1}}' {input} >{output}
|
|
143 """
|
|
144
|
|
145 rule bedtools_intersect:
|
|
146 input:
|
|
147 strains="outputs/genomes/strains.txt",
|
|
148 bedfiles=expand("outputs/genomes/{sample}.nuc.pangenome.gaf.bed", sample=SAMPLES)
|
|
149 output:
|
|
150 "outputs/pav_matrix.tsv",
|
|
151 shell:
|
|
152 """
|
|
153 perl $PANEX_PATH/Perl/GeneratePAVfromBed.pl {input.strains} outputs/genomes {output}
|
|
154 """
|
|
155
|
|
156
|
|
157 rule cog:
|
|
158 input:
|
|
159 pav="outputs/pav_matrix.tsv"
|
|
160 output:
|
|
161 cog="outputs/cog_output.txt",
|
|
162 cogstat="outputs/cog_stats.txt",
|
|
163 cogstat2="outputs/cog_stats2.txt",
|
|
164 cogofclusters="outputs/cog_of_clusters.txt"
|
|
165 shell:
|
|
166 """
|
|
167 perl $PANEX_PATH/Perl/GetCogOfCluster.pl {input} outputs/genomes {output.cog} {output.cogstat} {output.cogstat2} {output.cogofclusters} outputs/genomes/strains.txt
|
|
168 """
|
|
169
|
|
170 rule heatmap_upset:
|
|
171 input:
|
|
172 pav="outputs/pav_matrix.tsv"
|
|
173 output:
|
|
174 heatmap="outputs/heatmap.svg.gz",
|
|
175 html="outputs/heatmap.svg.heatmap_plotly.html",
|
|
176 upsetr="outputs/upsetr.svg",
|
|
177 binpav="outputs/heatmap.svg.pangenes_01matrix.txt"
|
|
178 shell:
|
|
179 """
|
|
180 perl $PANEX_PATH/Perl/GenerateHeatmapFromPAV.pl {input.pav} outputs/heatmap.svg
|
|
181 mv outputs/heatmap.svg.upsetr.svg {output.upsetr}
|
|
182 """
|
|
183
|
|
184 rule micropan:
|
|
185 input:
|
|
186 binpav="outputs/heatmap.svg.pangenes_01matrix.txt"
|
|
187 output:
|
|
188 txt="outputs/rarefaction_curves.txt",
|
|
189 pdf="outputs/rarefaction_curves.pdf",
|
|
190 svg="outputs/rarefaction_curves.svg",
|
|
191 heaps="outputs/heaps.tsv"
|
|
192 shell:
|
|
193 """
|
|
194 Rscript $PANEX_PATH/R/micropan_rarefaction.R -f {input.binpav} -p {output.pdf} -a {output.heaps} -o {output.txt}
|
|
195 pdf2svg {output.pdf} {output.svg}
|
|
196 """
|