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 rule final:
|
|
24 input:
|
|
25 "outputs/pav_matrix.csv",
|
|
26 "outputs/GCskew.txt",
|
|
27 "outputs/pav_matrix.tsv",
|
|
28 "outputs/heatmap.svg",
|
|
29 "outputs/cog_output.txt"
|
|
30
|
|
31
|
|
32 rule ncbi_datasets:
|
|
33 input:
|
|
34 "config.yaml"
|
|
35 output:
|
|
36 expand("outputs/genomes/{sample}.fasta", sample=SAMPLES),
|
|
37 expand("outputs/genomes/{sample}.gb", sample=SAMPLES),
|
|
38 expand("outputs/genomes/{sample}.pep", sample=SAMPLES),
|
|
39 genomes="outputs/genomes/genomes.txt"
|
|
40 shell:
|
|
41 """
|
|
42 perl $PANEX_PATH/Perl/get_data.pl {input} outputs/genomes
|
|
43 """
|
|
44
|
|
45
|
|
46 rule gcskew:
|
|
47 input:
|
|
48 "outputs/genomes/{sample}.fasta"
|
|
49 output:
|
|
50 "outputs/genomes/{sample}.fasta.gcskew.txt"
|
|
51 shell:
|
|
52 """
|
|
53 python3 $PANEX_PATH/SkewIT/src/gcskew.py -i {input} -o {input}.gcskew.txt -k 1000 -w 1000
|
|
54 """
|
|
55
|
|
56 rule concat_gcskew:
|
|
57 input:
|
|
58 expand("outputs/genomes/{sample}.fasta.gcskew.txt", sample=SAMPLES)
|
|
59 output:
|
|
60 out2="outputs/GCskew.txt"
|
|
61 shell:
|
|
62 """
|
|
63 cat {input} >>{output.out2}
|
|
64 """
|
|
65
|
|
66 rule genbank2gff3:
|
|
67 input:
|
|
68 "outputs/genomes/{sample}.gb"
|
|
69 output:
|
|
70 gff1="outputs/genomes/{sample}.gb.gff",
|
|
71 gff2="outputs/genomes/{sample}.gb.rmdup.gff"
|
|
72 shell:
|
|
73 """
|
|
74 perl $PANEX_PATH/Perl/bp_genbank2gff3.pl -o outputs/genomes {input}
|
|
75 perl $PANEX_PATH/Perl/remove_duplicates_in_gff.pl {output.gff1} {output.gff2}
|
|
76 """
|
|
77
|
|
78
|
|
79 rule pgap:
|
|
80 input:
|
|
81 expand("outputs/genomes/{sample}.pep", sample=SAMPLES)
|
|
82 params:
|
|
83 identity=80
|
|
84 output:
|
|
85 pav="outputs/pav_matrix.csv",
|
|
86 newick="outputs/accessory_binary_genes.fa.newick"
|
|
87 shell:
|
|
88 """
|
|
89 listforpgap=$(cat genbank_ids | tr '\n' '+')
|
|
90 alias consense='/usr/bin/phylip consense';alias dnaml='phylip /usr/bin/dnaml'
|
|
91 cp -rf /usr/local/bin/PGAP-1.2.1/Blast_Filter.pl .
|
|
92 perl /usr/local/bin/PGAP-1.2.1/PGAP.pl --strains $listforpgap --cluster --evolution --method GF --input outputs/genomes --output outputs/pgap_outdir
|
|
93 cp -rf outputs/pgap_outdir/1.Orthologs_Cluster.txt {output.pav}
|
|
94 cp -rf outputs/pgap_outdir/4.PanBased.Neighbor-joining.tree {output.newick}
|
|
95 """
|
|
96
|
|
97 rule convert_matrix:
|
|
98 input:
|
|
99 pav="outputs/pav_matrix.csv"
|
|
100 output:
|
|
101 "outputs/pav_matrix.tsv"
|
|
102 shell:
|
|
103 """
|
|
104 perl $PANEX_PATH/Perl/ConvertPGAPMatrix.pl outputs/genomes outputs/pav_matrix.csv outputs/pav_matrix.tsv outputs/genomes/strains.txt
|
|
105 """
|
|
106
|
|
107 rule heatmap_upset:
|
|
108 input:
|
|
109 pav="outputs/pav_matrix.tsv"
|
|
110 output:
|
|
111 heatmap="outputs/heatmap.svg",
|
|
112 upsetr="outputs/upsetr.svg"
|
|
113 shell:
|
|
114 """
|
|
115 perl $PANEX_PATH/Perl/GenerateHeatmapFromPAV.pl {input.pav} {output.heatmap}
|
|
116 mv {output.heatmap}.upsetr.svg {output.upsetr}
|
|
117 """
|
|
118
|
|
119 rule cog:
|
|
120 input:
|
|
121 pav="outputs/pav_matrix.tsv"
|
|
122 output:
|
|
123 cog="outputs/cog_output.txt",
|
|
124 cogstat="outputs/cog_stats.txt",
|
|
125 cogstat2="outputs/cog_stats2.txt",
|
|
126 cogofclusters="outputs/cog_of_clusters.txt"
|
|
127 shell:
|
|
128 """
|
|
129 perl $PANEX_PATH/Perl/GetCogOfCluster.pl {input} outputs/genomes {output.cog} {output.cogstat} {output.cogstat2} {output.cogofclusters} outputs/genomes/strains.txt
|
|
130 """
|
|
131
|