view Snakemake_files/Snakefile_wget_PGAP_heatmap_upset_COG @ 3:e42d30da7a74 draft

Uploaded

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()
            SAMPLES.append(words[1])


rule final:
 input:
  "outputs/pav_matrix.csv",
  "outputs/GCskew.txt",
  "outputs/pav_matrix.tsv",
  "outputs/heatmap.svg",
  "outputs/cog_output.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}.pep", 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 pgap:
    input:
        expand("outputs/genomes/{sample}.pep", sample=SAMPLES)
    params:
        identity=80
    output:
        pav="outputs/pav_matrix.csv",
        newick="outputs/accessory_binary_genes.fa.newick"
    shell:
        """
        listforpgap=$(cat genbank_ids | tr '\n' '+')
        alias consense='/usr/bin/phylip consense';alias dnaml='phylip /usr/bin/dnaml' 
        cp -rf /usr/local/bin/PGAP-1.2.1/Blast_Filter.pl .
        perl /usr/local/bin/PGAP-1.2.1/PGAP.pl --strains $listforpgap --cluster --evolution --method GF --input outputs/genomes --output outputs/pgap_outdir
        cp -rf outputs/pgap_outdir/1.Orthologs_Cluster.txt {output.pav}
        cp -rf outputs/pgap_outdir/4.PanBased.Neighbor-joining.tree {output.newick}
        """

rule convert_matrix:
    input:
        pav="outputs/pav_matrix.csv"
    output:
        "outputs/pav_matrix.tsv"
    shell:
        """
        perl $PANEX_PATH/Perl/ConvertPGAPMatrix.pl outputs/genomes outputs/pav_matrix.csv outputs/pav_matrix.tsv outputs/genomes/strains.txt
        """

rule heatmap_upset:
    input:
        pav="outputs/pav_matrix.tsv"
    output:
        heatmap="outputs/heatmap.svg",
        upsetr="outputs/upsetr.svg"
    shell:
        """
        perl $PANEX_PATH/Perl/GenerateHeatmapFromPAV.pl {input.pav} {output.heatmap}
        mv {output.heatmap}.upsetr.svg {output.upsetr}
        """

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