Mercurial > repos > dereeper > pangenome_explorer
comparison Snakemake_files/Snakefile_wget_cactus_heatmap_upset_COG2 @ 3:e42d30da7a74 draft
Uploaded
author | dereeper |
---|---|
date | Thu, 30 May 2024 11:52:25 +0000 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
2:97e4e3e818b6 | 3:e42d30da7a74 |
---|---|
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 rule final: | |
24 input: | |
25 "outputs/GCskew.txt", | |
26 "outputs/cactus_outdir/output_cactus.sv.gfa.gz", | |
27 "outputs/cactus_outdir/full.gfa", | |
28 "outputs/cactus_outdir/full.gfa.png", | |
29 "outputs/pav_matrix.tsv", | |
30 "outputs/heatmap.svg.gz", | |
31 "outputs/cog_output.txt", | |
32 "outputs/fastani.out", | |
33 "outputs/rarefaction_curves.txt" | |
34 | |
35 rule ncbi_datasets: | |
36 input: | |
37 "config.yaml" | |
38 output: | |
39 expand("outputs/genomes/{sample}.fasta", sample=SAMPLES), | |
40 expand("outputs/genomes/{sample}.gb", sample=SAMPLES), | |
41 expand("outputs/genomes/{sample}.prt", sample=SAMPLES), | |
42 expand("outputs/genomes/{sample}.nuc", sample=SAMPLES), | |
43 genomes="outputs/genomes/genomes.txt", | |
44 strains="outputs/genomes/strains.txt" | |
45 shell: | |
46 """ | |
47 perl $PANEX_PATH/Perl/get_data.pl {input} outputs/genomes | |
48 """ | |
49 | |
50 | |
51 rule gcskew: | |
52 input: | |
53 "outputs/genomes/{sample}.fasta" | |
54 output: | |
55 "outputs/genomes/{sample}.fasta.gcskew.txt" | |
56 shell: | |
57 """ | |
58 python3 $PANEX_PATH/SkewIT/src/gcskew.py -i {input} -o {input}.gcskew.txt -k 1000 -w 1000 | |
59 """ | |
60 | |
61 rule concat_gcskew: | |
62 input: | |
63 expand("outputs/genomes/{sample}.fasta.gcskew.txt", sample=SAMPLES) | |
64 output: | |
65 out2="outputs/GCskew.txt" | |
66 shell: | |
67 """ | |
68 cat {input} >>{output.out2} | |
69 """ | |
70 | |
71 rule genbank2gff3: | |
72 input: | |
73 "outputs/genomes/{sample}.gb" | |
74 output: | |
75 gff1="outputs/genomes/{sample}.gb.gff", | |
76 gff2="outputs/genomes/{sample}.gb.rmdup.gff" | |
77 shell: | |
78 """ | |
79 perl $PANEX_PATH/Perl/bp_genbank2gff3.pl -o outputs/genomes {input} | |
80 perl $PANEX_PATH/Perl/remove_duplicates_in_gff.pl {output.gff1} {output.gff2} | |
81 """ | |
82 | |
83 rule cactus: | |
84 input: | |
85 expand("outputs/genomes/{sample}.fasta", sample=SAMPLES), | |
86 expand("outputs/genomes/{sample}.gb.rmdup.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 rule create_gene_paths: | |
100 input: | |
101 gff="outputs/genomes/{sample}.gb.rmdup.gff", | |
102 gfa="outputs/cactus_outdir/full.gfa", | |
103 output: | |
104 basename="outputs/genomes/{sample}.gene_segments", | |
105 gene_length="outputs/genomes/{sample}.gene_segments.gene_length.txt", | |
106 bed="outputs/genomes/{sample}.gene_segments.bed", | |
107 shell: | |
108 """ | |
109 perl $PANEX_PATH/Perl/CreateGenePathsFromGFA.pl {input.gfa} {input.gff} {output.basename} | |
110 """ | |
111 | |
112 rule odgi: | |
113 input: | |
114 "outputs/cactus_outdir/full.gfa", | |
115 output: | |
116 og="outputs/cactus_outdir/full.gfa.og", | |
117 png="outputs/cactus_outdir/full.gfa.png", | |
118 shell: | |
119 """ | |
120 odgi build -g {input} -o {output.og} | |
121 odgi viz -i {output.og} -o {output.png} | |
122 """ | |
123 | |
124 rule bedtools_intersect: | |
125 input: | |
126 strains="outputs/genomes/strains.txt", | |
127 bedfiles=expand("outputs/genomes/{sample}.gene_segments.bed", sample=SAMPLES) | |
128 output: | |
129 "outputs/pav_matrix.tsv", | |
130 shell: | |
131 """ | |
132 perl $PANEX_PATH/Perl/GeneratePAVfromBed.pl {input.strains} outputs/genomes {output} | |
133 """ | |
134 | |
135 | |
136 rule cog: | |
137 input: | |
138 pav="outputs/pav_matrix.tsv" | |
139 output: | |
140 cog="outputs/cog_output.txt", | |
141 cogstat="outputs/cog_stats.txt", | |
142 cogstat2="outputs/cog_stats2.txt", | |
143 cogofclusters="outputs/cog_of_clusters.txt" | |
144 shell: | |
145 """ | |
146 perl $PANEX_PATH/Perl/GetCogOfCluster.pl {input} outputs/genomes {output.cog} {output.cogstat} {output.cogstat2} {output.cogofclusters} outputs/genomes/strains.txt | |
147 """ | |
148 | |
149 rule heatmap_upset: | |
150 input: | |
151 pav="outputs/pav_matrix.tsv" | |
152 output: | |
153 heatmap="outputs/heatmap.svg.gz", | |
154 html="outputs/heatmap.svg.heatmap_plotly.html", | |
155 upsetr="outputs/upsetr.svg", | |
156 binpav="outputs/heatmap.svg.pangenes_01matrix.txt" | |
157 shell: | |
158 """ | |
159 perl $PANEX_PATH/Perl/GenerateHeatmapFromPAV.pl {input.pav} outputs/heatmap.svg | |
160 mv outputs/heatmap.svg.upsetr.svg {output.upsetr} | |
161 """ | |
162 | |
163 rule micropan: | |
164 input: | |
165 binpav="outputs/heatmap.svg.pangenes_01matrix.txt" | |
166 output: | |
167 txt="outputs/rarefaction_curves.txt", | |
168 pdf="outputs/rarefaction_curves.pdf", | |
169 svg="outputs/rarefaction_curves.svg", | |
170 heaps="outputs/heaps.tsv" | |
171 shell: | |
172 """ | |
173 Rscript $PANEX_PATH/R/micropan_rarefaction.R -f {input.binpav} -p {output.pdf} -a {output.heaps} -o {output.txt} | |
174 pdf2svg {output.pdf} {output.svg} | |
175 """ | |
176 | |
177 rule fastani: | |
178 input: | |
179 "outputs/genomes/genomes.txt" | |
180 output: | |
181 out="outputs/fastani.out", | |
182 matrix="outputs/fastani.out.matrix", | |
183 completematrix="outputs/fastani.out.matrix.complete", | |
184 pdf="outputs/fastani.out.pdf", | |
185 svg="outputs/fastani.out.svg" | |
186 shell: | |
187 """ | |
188 fastANI --rl {input} --ql {input} -o {output.out} -t 4 --matrix | |
189 perl $PANEX_PATH/Perl/convertANI.pl {output.matrix} outputs/genomes/genomes2.txt >{output.completematrix} | |
190 Rscript $PANEX_PATH/R/heatmap_ani.R -f {output.completematrix} -o {output.pdf} | |
191 pdf2svg {output.pdf} {output.svg} | |
192 """ |