Mercurial > repos > dereeper > pangenome_explorer
view Snakemake_files/Snakefile_wget_cactus_heatmap_upset_COG2 @ 13:152d7c43478b draft default tip
Uploaded
author | dereeper |
---|---|
date | Thu, 30 May 2024 20:07:55 +0000 |
parents | e42d30da7a74 |
children |
line wrap: on
line source
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/GCskew.txt", "outputs/cactus_outdir/output_cactus.sv.gfa.gz", "outputs/cactus_outdir/full.gfa", "outputs/cactus_outdir/full.gfa.png", "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), expand("outputs/genomes/{sample}.nuc", sample=SAMPLES), genomes="outputs/genomes/genomes.txt", strains="outputs/genomes/strains.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 cactus: input: expand("outputs/genomes/{sample}.fasta", sample=SAMPLES), expand("outputs/genomes/{sample}.gb.rmdup.gff", sample=SAMPLES), params: reference=SAMPLES[0] output: gfasv="outputs/cactus_outdir/output_cactus.sv.gfa.gz", gfafull="outputs/cactus_outdir/full.gfa" shell: """ mkdir outputs/cactus_work cactus-pangenome outputs/jobstore outputs/genomes/seqfile --outDir outputs/cactus_outdir --outName output_cactus --reference {params.reference} --workDir outputs/cactus_work --vcf --mapCores 12 --gfa clip full --viz zcat outputs/cactus_outdir/output_cactus.full.gfa.gz >{output.gfafull} """ rule create_gene_paths: input: gff="outputs/genomes/{sample}.gb.rmdup.gff", gfa="outputs/cactus_outdir/full.gfa", output: basename="outputs/genomes/{sample}.gene_segments", gene_length="outputs/genomes/{sample}.gene_segments.gene_length.txt", bed="outputs/genomes/{sample}.gene_segments.bed", shell: """ perl $PANEX_PATH/Perl/CreateGenePathsFromGFA.pl {input.gfa} {input.gff} {output.basename} """ rule odgi: input: "outputs/cactus_outdir/full.gfa", output: og="outputs/cactus_outdir/full.gfa.og", png="outputs/cactus_outdir/full.gfa.png", shell: """ odgi build -g {input} -o {output.og} odgi viz -i {output.og} -o {output.png} """ rule bedtools_intersect: input: strains="outputs/genomes/strains.txt", bedfiles=expand("outputs/genomes/{sample}.gene_segments.bed", sample=SAMPLES) output: "outputs/pav_matrix.tsv", shell: """ perl $PANEX_PATH/Perl/GeneratePAVfromBed.pl {input.strains} outputs/genomes {output} """ 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} """