Mercurial > repos > dereeper > pangenome_explorer
diff Snakemake_files/Snakefile_wget_pggb_heatmap_upset_COG @ 3:e42d30da7a74 draft
Uploaded
author | dereeper |
---|---|
date | Thu, 30 May 2024 11:52:25 +0000 |
parents | |
children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Snakemake_files/Snakefile_wget_pggb_heatmap_upset_COG Thu May 30 11:52:25 2024 +0000 @@ -0,0 +1,197 @@ +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/pggb_out/all_genomes.fa.smooth.final.gfa", + "outputs/pav_matrix.tsv", + "outputs/heatmap.svg.gz", + "outputs/cog_output.txt", + "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", + shell: + """ + perl $PANEX_PATH/Perl/bp_genbank2gff3.pl -o outputs/genomes {input} + """ + +rule pggb: + input: + expand("outputs/genomes/{sample}.fasta", sample=SAMPLES), + expand("outputs/genomes/{sample}.gb.gff", sample=SAMPLES), + output: + gfa="outputs/pggb_out/all_genomes.fa.smooth.final.gfa", + png1="outputs/pggb_out/all_genomes.fa.lay.draw.png", + png2="outputs/pggb_out/all_genomes.fa.og.viz_multiqc.png", + vcf="outputs/all_genomes.vcf", + shell: + """ + samtools faidx outputs/genomes/all_genomes.fa + reference=$(head -1 outputs/genomes/strains.txt | awk '{{print $2}}') + pggb -i outputs/genomes/all_genomes.fa -o outputs/pggb_out -V $reference -m -M + mv outputs/pggb_out/all_genomes.*smooth.final.gfa {output.gfa} + mv outputs/pggb_out/all_genomes.*lay.draw.png {output.png1} + mv outputs/pggb_out/all_genomes.fa.*.og.viz_multiqc.png {output.png2} + mv outputs/pggb_out/all_genomes.*vcf {output.vcf} + """ + +rule odgi: + input: + gfa="outputs/pggb_out/all_genomes.fa.smooth.final.gfa", + output: + og="outputs/pggb_out/all_genomes.fa.smooth.final.gfa.og", + distance="outputs/pggb_out/all_genomes.fa.smooth.final.gfa.og.distance", + shell: + """ + odgi build -g {input} -o {output.og} + odgi similarity -i {output.og} -d >{output.distance} + """ + + +rule create_gene_paths: + input: + gff="outputs/genomes/{sample}.gb.gff", + gfa="outputs/pggb_out/all_genomes.fa.smooth.final.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 bedtools_intersect: + input: + strains="outputs/genomes/strains.txt", + bedfiles=expand("outputs/genomes/{sample}.gene_segments.bed", sample=SAMPLES) + params: + identity=30 + output: + "outputs/pav_matrix.tsv", + shell: + """ + perl $PANEX_PATH/Perl/GeneratePAVfromBed.pl {input.strains} outputs/genomes {output} {params.identity} + """ + + +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} + """