diff Snakemake_files/Snakefile_wget_panacota_heatmap_upset_COG @ 3:e42d30da7a74 draft

Uploaded

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Snakemake_files/Snakefile_wget_panacota_heatmap_upset_COG	Thu May 30 11:52:25 2024 +0000
@@ -0,0 +1,180 @@
+import glob
+import os
+import shutil
+
+import yaml
+configfile: "config.yaml"
+
+
+SAMPLES = []
+with open("config.yaml", "r") as yaml_file:
+    genome_data = yaml.safe_load(yaml_file)
+    if "ids" in genome_data.keys():
+        for id in genome_data["ids"]:
+            SAMPLES.append(id)
+    if "input_genbanks" in genome_data.keys():
+        for gb_path in genome_data["input_genbanks"]:
+            cmd = "grep 'ACCESSION' "+gb_path
+            returned_value = subprocess.getoutput(cmd)
+            words = returned_value.split()
+            name = words[1]
+            #words = returned_value.split()
+            #SAMPLES.append(words[1])
+            words = gb_path.split("/")
+            genbank_file_name = words[-1]
+            x = genbank_file_name.replace(".", "_")
+            cmd = "sed -i 's/ACCESSION   "+name+"/ACCESSION   "+x+"/g' "+gb_path
+            subprocess.getoutput(cmd)
+
+            cmd_locus = "grep 'LOCUS' "+gb_path
+            returned_value = subprocess.getoutput(cmd_locus)
+            words_locus = returned_value.split()
+            name_locus = words_locus[1]
+            cmd = "sed -i 's/LOCUS       "+name_locus+"/LOCUS       "+x+"/' "+gb_path
+            subprocess.getoutput(cmd)
+
+            cmd = "grep 'ACCESSION' "+gb_path
+            returned_value = subprocess.getoutput(cmd)
+            words = returned_value.split()
+            SAMPLES.append(words[1])
+
+
+rule final:
+ input:
+  "outputs/GCskew.txt",
+  "outputs/pav_matrix.tsv",
+  "outputs/heatmap.svg.gz",
+  "outputs/cog_output.txt",
+  "outputs/fastani.out",
+  "outputs/rarefaction_curves.txt"
+
+
+rule ncbi_datasets:
+    input:
+        "config.yaml"
+    output:
+        expand("outputs/genomes/{sample}.fasta", sample=SAMPLES),
+        expand("outputs/genomes/{sample}.gb", sample=SAMPLES),
+        expand("outputs/genomes/{sample}.prt", sample=SAMPLES),
+        genomes="outputs/genomes/genomes.txt"
+    shell:
+        """
+        perl $PANEX_PATH/Perl/get_data.pl {input} outputs/genomes
+        """
+
+
+rule gcskew:
+    input:
+        "outputs/genomes/{sample}.fasta" 
+    output:
+        "outputs/genomes/{sample}.fasta.gcskew.txt"
+    shell:
+        """
+        python3 $PANEX_PATH/SkewIT/src/gcskew.py -i {input} -o {input}.gcskew.txt -k 1000 -w 1000
+        """ 
+
+rule concat_gcskew:
+    input:
+        expand("outputs/genomes/{sample}.fasta.gcskew.txt", sample=SAMPLES)
+    output:
+        out2="outputs/GCskew.txt"
+    shell:
+        """
+         cat {input} >>{output.out2}
+        """
+
+rule genbank2gff3:
+    input:
+        "outputs/genomes/{sample}.gb"
+    output:
+        gff1="outputs/genomes/{sample}.gb.gff",
+        gff2="outputs/genomes/{sample}.gb.rmdup.gff"
+    shell:
+        """
+        perl $PANEX_PATH/Perl/bp_genbank2gff3.pl -o outputs/genomes {input}
+        perl $PANEX_PATH/Perl/remove_duplicates_in_gff.pl {output.gff1} {output.gff2}
+        """
+
+rule panacota:
+    input:
+        expand("outputs/genomes/{sample}.prt", sample=SAMPLES),
+        expand("outputs/genomes/{sample}.gb.gff", sample=SAMPLES)
+    output:
+        pav="outputs/panacota/PanGenome-mydataset.All.prt.lst"
+    params:
+        identity=80
+    shell:
+        """
+        cat outputs/genomes/*prt >outputs/genomes/mydataset.All.prt
+        PanACoTA pangenome -i 0.{params.identity} -l outputs/genomes/mydataset.All.prt -n mydataset -d outputs/genomes -o outputs/panacota
+        mv outputs/panacota/PanGenome-mydataset.All.prt-clust-*-mode1.lst {output}
+        """
+
+rule convert_matrix:
+    input:
+        pav="outputs/panacota/PanGenome-mydataset.All.prt.lst"
+    output:
+        "outputs/pav_matrix.tsv"
+    shell:
+        """
+        perl $PANEX_PATH/Perl/ConvertPanacotaMatrix.pl outputs/genomes {input} {output} outputs/genomes/strains.txt 
+        """
+
+
+rule cog:
+    input:
+        pav="outputs/pav_matrix.tsv"
+    output:
+        cog="outputs/cog_output.txt",
+        cogstat="outputs/cog_stats.txt",
+        cogstat2="outputs/cog_stats2.txt",
+        cogofclusters="outputs/cog_of_clusters.txt"
+    shell:
+        """
+        perl $PANEX_PATH/Perl/GetCogOfCluster.pl {input} outputs/genomes {output.cog} {output.cogstat} {output.cogstat2} {output.cogofclusters} outputs/genomes/strains.txt
+        """
+        
+rule heatmap_upset:
+    input:
+        pav="outputs/pav_matrix.tsv"
+    output:
+        heatmap="outputs/heatmap.svg.gz",
+        html="outputs/heatmap.svg.heatmap_plotly.html",
+        upsetr="outputs/upsetr.svg",
+        binpav="outputs/heatmap.svg.pangenes_01matrix.txt"
+    shell:
+        """
+        perl $PANEX_PATH/Perl/GenerateHeatmapFromPAV.pl {input.pav} outputs/heatmap.svg
+        mv outputs/heatmap.svg.upsetr.svg {output.upsetr}
+        """
+
+rule micropan:
+    input:
+        binpav="outputs/heatmap.svg.pangenes_01matrix.txt"
+    output:
+        txt="outputs/rarefaction_curves.txt",
+        pdf="outputs/rarefaction_curves.pdf",
+        svg="outputs/rarefaction_curves.svg",
+        heaps="outputs/heaps.tsv"
+    shell:
+        """
+        Rscript $PANEX_PATH/R/micropan_rarefaction.R -f {input.binpav} -p {output.pdf} -a {output.heaps} -o {output.txt}
+        pdf2svg {output.pdf} {output.svg}
+        """
+
+rule fastani:
+    input:
+        "outputs/genomes/genomes.txt"
+    output:
+        out="outputs/fastani.out",
+        matrix="outputs/fastani.out.matrix",
+        completematrix="outputs/fastani.out.matrix.complete",
+        pdf="outputs/fastani.out.pdf",
+        svg="outputs/fastani.out.svg"
+    shell:
+        """
+        fastANI --rl {input} --ql {input} -o {output.out} -t 4 --matrix
+        perl $PANEX_PATH/Perl/convertANI.pl {output.matrix} outputs/genomes/genomes2.txt >{output.completematrix}
+        Rscript $PANEX_PATH/R/heatmap_ani.R -f {output.completematrix} -o {output.pdf}
+        pdf2svg {output.pdf} {output.svg}
+        """