3
|
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 """
|