Mercurial > repos > dereeper > pangenome_explorer
comparison Snakemake_files/Snakefile_wget_orthofinder_heatmap_upset_COG @ 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 SAMPLES = [] | |
9 with open("config.yaml", "r") as yaml_file: | |
10 genome_data = yaml.safe_load(yaml_file) | |
11 if "ids" in genome_data.keys(): | |
12 for id in genome_data["ids"]: | |
13 SAMPLES.append(id) | |
14 if "input_genbanks" in genome_data.keys(): | |
15 for gb_path in genome_data["input_genbanks"]: | |
16 | |
17 cmd = "grep 'ACCESSION' "+gb_path | |
18 returned_value = subprocess.getoutput(cmd) | |
19 words = returned_value.split() | |
20 name = words[1] | |
21 #words = returned_value.split() | |
22 #SAMPLES.append(words[1]) | |
23 words = gb_path.split("/") | |
24 genbank_file_name = words[-1] | |
25 x = genbank_file_name.replace(".", "_") | |
26 cmd = "sed -i 's/ACCESSION "+name+"/ACCESSION "+x+"/g' "+gb_path | |
27 subprocess.getoutput(cmd) | |
28 | |
29 cmd_locus = "grep 'LOCUS' "+gb_path | |
30 returned_value = subprocess.getoutput(cmd_locus) | |
31 words_locus = returned_value.split() | |
32 name_locus = words_locus[1] | |
33 cmd = "sed -i 's/LOCUS "+name_locus+"/LOCUS "+x+"/' "+gb_path | |
34 subprocess.getoutput(cmd) | |
35 | |
36 cmd = "grep 'ACCESSION' "+gb_path | |
37 returned_value = subprocess.getoutput(cmd) | |
38 words = returned_value.split() | |
39 SAMPLES.append(words[1]) | |
40 | |
41 | |
42 rule final: | |
43 input: | |
44 "outputs/pav_matrix.tsv", | |
45 "outputs/heatmap.svg", | |
46 "outputs/rarefaction_curves.txt" | |
47 | |
48 | |
49 rule ncbi_datasets: | |
50 input: | |
51 "config.yaml" | |
52 output: | |
53 expand("outputs/genomes/{sample}.fasta", sample=SAMPLES), | |
54 expand("outputs/genomes/{sample}.gb", sample=SAMPLES), | |
55 expand("outputs/genomes/{sample}.prt", sample=SAMPLES), | |
56 expand("outputs/genomes/{sample}.nuc", sample=SAMPLES), | |
57 expand("outputs/genomes/{sample}.pep", sample=SAMPLES), | |
58 genomes="outputs/genomes/genomes.txt", | |
59 strains="outputs/genomes/strains.txt" | |
60 shell: | |
61 """ | |
62 perl $PANEX_PATH/Perl/get_data.pl {input} outputs/genomes | |
63 """ | |
64 | |
65 rule genbank2gff3: | |
66 input: | |
67 "outputs/genomes/{sample}.gb" | |
68 output: | |
69 gff1="outputs/genomes/{sample}.gb.gff", | |
70 gff2="outputs/genomes/{sample}.gb.rmdup.gff", | |
71 shell: | |
72 """ | |
73 perl $PANEX_PATH/Perl/bp_genbank2gff3.pl -o outputs/genomes {input} | |
74 perl $PANEX_PATH/Perl/remove_duplicates_in_gff.pl {output.gff1} {output.gff2} | |
75 mv outputs/genomes/*.faa outputs/proteomes | |
76 """ | |
77 | |
78 rule orthofinder: | |
79 envmodules: | |
80 "bioinfo/orthofinder/2.5.4" | |
81 input: | |
82 expand("outputs/genomes/{sample}.pep", sample=SAMPLES) | |
83 params: | |
84 identity=80 | |
85 output: | |
86 pav="outputs/Orthogroups.tsv" | |
87 shell: | |
88 """ | |
89 mkdir -p outputs/proteomes | |
90 cp -rf outputs/genomes/*.pep outputs/proteomes | |
91 orthofinder -f outputs/proteomes | |
92 mv outputs/proteomes/OrthoFinder/R*/Orthogroups/Orthogroups.tsv outputs/Orthogroups.tsv | |
93 """ | |
94 | |
95 rule convert_matrix: | |
96 input: | |
97 pav="outputs/Orthogroups.tsv" | |
98 output: | |
99 "outputs/pav_matrix.tsv" | |
100 shell: | |
101 """ | |
102 perl $PANEX_PATH/Perl/ConvertOrthofinderMatrix.pl outputs/genomes {input} {output} outputs/genomes/strains.txt | |
103 """ | |
104 | |
105 rule cog: | |
106 input: | |
107 pav="outputs/pav_matrix.tsv" | |
108 output: | |
109 cog="outputs/cog_output.txt", | |
110 cogstat="outputs/cog_stats.txt", | |
111 cogstat2="outputs/cog_stats2.txt", | |
112 cogofclusters="outputs/cog_of_clusters.txt" | |
113 shell: | |
114 """ | |
115 perl $PANEX_PATH/Perl/GetCogOfCluster.pl {input} outputs/genomes {output.cog} {output.cogstat} {output.cogstat2} {output.cogofclusters} outputs/genomes/strains.txt | |
116 """ | |
117 | |
118 rule heatmap_upset: | |
119 input: | |
120 pav="outputs/pav_matrix.tsv", | |
121 cogofclusters="outputs/cog_of_clusters.txt" | |
122 output: | |
123 heatmap="outputs/heatmap.svg", | |
124 upsetr="outputs/upsetr.svg", | |
125 binpav="outputs/heatmap.svg.pangenes_01matrix.txt" | |
126 shell: | |
127 """ | |
128 perl $PANEX_PATH/Perl/GenerateHeatmapFromPAV.pl {input.pav} {output.heatmap} | |
129 mv {output.heatmap}.upsetr.svg {output.upsetr} | |
130 """ | |
131 | |
132 rule micropan: | |
133 input: | |
134 binpav="outputs/heatmap.svg.pangenes_01matrix.txt" | |
135 output: | |
136 txt="outputs/rarefaction_curves.txt", | |
137 pdf="outputs/rarefaction_curves.pdf", | |
138 svg="outputs/rarefaction_curves.svg", | |
139 heaps="outputs/heaps.tsv" | |
140 shell: | |
141 """ | |
142 Rscript $PANEX_PATH/R/micropan_rarefaction.R -f {input.binpav} -p {output.pdf} -a {output.heaps} -o {output.txt} | |
143 pdf2svg {output.pdf} {output.svg} | |
144 """ |