Skip to content
Snippets Groups Projects
Commit b09393be authored by james94's avatar james94
Browse files

merqury syntax, using stable (non-experimental) meryl

parent aa915446
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 24 additions and 17 deletions
rules/rulesNew.smk
+ 24
− 17
View file @ b09393be
...@@ -146,7 +146,8 @@ rule meryl_R1: ...@@ -146,7 +146,8 @@ rule meryl_R1:
read1= rules.combine_illuminaReads_R1.output read1= rules.combine_illuminaReads_R1.output
# check = os.path.join(config['Results'], "{sample}" + "/busco5/short_summary.specific." + config['busco5Lineage'] + "_odb10." + "{sample}" + ".txt") # check = os.path.join(config['Results'], "{sample}" + "/busco5/short_summary.specific." + config['busco5Lineage'] + "_odb10." + "{sample}" + ".txt")
params: params:
script = os.path.join(workflow.basedir, "programs/meryl/*/bin/meryl"), # script = os.path.join(workflow.basedir, "programs/meryl/build/bin/meryl"),
script = os.path.join(workflow.basedir, "programs/meryl-1.0/Linux-amd64/bin/meryl"),
kmer = 21, kmer = 21,
r1= "{sample}" + "_R1.21.meryl", r1= "{sample}" + "_R1.21.meryl",
path= os.path.join(config['Results'], "{sample}" +"/2_QVstats_merylAndMerqury/") path= os.path.join(config['Results'], "{sample}" +"/2_QVstats_merylAndMerqury/")
...@@ -170,7 +171,8 @@ rule meryl_R2: ...@@ -170,7 +171,8 @@ rule meryl_R2:
read2= rules.combine_illuminaReads_R2.output read2= rules.combine_illuminaReads_R2.output
# check = os.path.join(config['Results'], "{sample}" + "/busco5/short_summary.specific." + config['busco5Lineage'] + "_odb10." + "{sample}" + ".txt") # check = os.path.join(config['Results'], "{sample}" + "/busco5/short_summary.specific." + config['busco5Lineage'] + "_odb10." + "{sample}" + ".txt")
params: params:
script = os.path.join(workflow.basedir, "programs/meryl/*/bin/meryl"), # script = os.path.join(workflow.basedir, "programs/meryl/build/bin/meryl"),
script = os.path.join(workflow.basedir, "programs/meryl-1.0/Linux-amd64/bin/meryl"),
kmer = 21, kmer = 21,
r2= "{sample}" + "_R2.21.meryl", r2= "{sample}" + "_R2.21.meryl",
path= os.path.join(config['Results'], "{sample}" +"/2_QVstats_merylAndMerqury/") path= os.path.join(config['Results'], "{sample}" +"/2_QVstats_merylAndMerqury/")
...@@ -195,7 +197,8 @@ rule meryl: ...@@ -195,7 +197,8 @@ rule meryl:
read2=os.path.join(config['Results'], "{sample}" +"/2_QVstats_merylAndMerqury/" + "{sample}" + "_R2.21.meryl") read2=os.path.join(config['Results'], "{sample}" +"/2_QVstats_merylAndMerqury/" + "{sample}" + "_R2.21.meryl")
# check = os.path.join(config['Results'], "{sample}" + "/busco5/short_summary.specific." + config['busco5Lineage'] + "_odb10." + "{sample}" + ".txt") # check = os.path.join(config['Results'], "{sample}" + "/busco5/short_summary.specific." + config['busco5Lineage'] + "_odb10." + "{sample}" + ".txt")
params: params:
script = os.path.join(workflow.basedir, "programs/meryl/*/bin/meryl"), # script = os.path.join(workflow.basedir, "programs/meryl/build/bin/meryl"),
script = os.path.join(workflow.basedir, "programs/meryl-1.0/Linux-amd64/bin/meryl"),
kmer = 21, kmer = 21,
# r1= "{sample}" + "_R1.21.meryl", # r1= "{sample}" + "_R1.21.meryl",
# r2= "{sample}" + "_R2.21.meryl", # r2= "{sample}" + "_R2.21.meryl",
...@@ -226,17 +229,18 @@ rule merqury: ...@@ -226,17 +229,18 @@ rule merqury:
script = os.path.join(workflow.basedir, "programs/merqury-1.1/merqury.sh"), script = os.path.join(workflow.basedir, "programs/merqury-1.1/merqury.sh"),
outFile= "{sample}" + "_merq", outFile= "{sample}" + "_merq",
# outFile= os.path.join(config['Results'], "{sample}" +"/QVstats_merylAndMerqury/{sample}_merq"), # outFile= os.path.join(config['Results'], "{sample}" +"/QVstats_merylAndMerqury/{sample}_merq"),
merylPath = os.path.join(workflow.basedir, "programs/meryl/*/bin:$PATH"), # merylPath = os.path.join(workflow.basedir, "programs/meryl/build/bin:$PATH"),
merylPath = os.path.join(workflow.basedir, "programs/meryl-1.0/Linux-amd64/bin:$PATH"),
outPath= os.path.join(config['Results'],"{sample}" + "/2_QVstats_merylAndMerqury"), outPath= os.path.join(config['Results'],"{sample}" + "/2_QVstats_merylAndMerqury"),
export= os.path.join(workflow.basedir, "programs/merqury-1.1") export= os.path.join(workflow.basedir, "programs/merqury-1.1")
# symlink=os.path.join(config['Results'], "{sample}" + "/{sample}.fasta")
threads: threads:
workflow.cores * 0.5 workflow.cores * 0.5
output: output:
report(os.path.join(config['Results'],"{sample}" + "/2_QVstats_merylAndMerqury/" + "{sample}" + "_merq.qv"), caption="../report/merqury/qvScore.rst", category="QV Estimate", subcategory="{sample}"), report(os.path.join(config['Results'],"{sample}" + "/2_QVstats_merylAndMerqury/" + "{sample}" + "_merq.qv"), caption="../report/merqury/qvScore.rst", category="QV Estimate", subcategory="{sample}"),
report(os.path.join(config['Results'],"{sample}" + "/2_QVstats_merylAndMerqury/" + "{sample}" + "_merq.{sample}_assembly.spectra-cn.st.png"), caption="../report/merqury/spectra-cn_stacked.rst", category="Copy Number Spectra Plots", subcategory="{sample}"), report(os.path.join(config['Results'],"{sample}" + "/2_QVstats_merylAndMerqury/" + "{sample}" + "_merq.{sample}.spectra-cn.st.png"), caption="../report/merqury/spectra-cn_stacked.rst", category="Copy Number Spectra Plots", subcategory="{sample}"),
report(os.path.join(config['Results'],"{sample}" + "/2_QVstats_merylAndMerqury/" + "{sample}" + "_merq.{sample}_assembly.spectra-cn.fl.png"), caption="../report/merqury/spectra-cn_filled.rst", category="Copy Number Spectra Plots", subcategory="{sample}"), report(os.path.join(config['Results'],"{sample}" + "/2_QVstats_merylAndMerqury/" + "{sample}" + "_merq.{sample}.spectra-cn.fl.png"), caption="../report/merqury/spectra-cn_filled.rst", category="Copy Number Spectra Plots", subcategory="{sample}"),
os.path.join(config['Results'], "{sample}" +"/2_QVstats_merylAndMerqury/" + "{sample}" + "_dB.21.hist") os.path.join(config['Results'], "{sample}" +"/2_QVstats_merylAndMerqury/" + "{sample}" + "_dB.21.hist"),
symlink=os.path.join(config['Results'], "{sample}" + "/{sample}.fasta")
log: log:
os.path.join(config['Results'],"{sample}" + "/2_QVstats_merylAndMerqury/logs/{sample}_merqury.log") os.path.join(config['Results'],"{sample}" + "/2_QVstats_merylAndMerqury/logs/{sample}_merqury.log")
priority: priority:
...@@ -246,10 +250,11 @@ rule merqury: ...@@ -246,10 +250,11 @@ rule merqury:
shell: shell:
""" """
cd {params.outPath} cd {params.outPath}
ln -s {input.assembly} {output.symlink}
export PATH={params.merylPath} export PATH={params.merylPath}
export MERQURY={params.export} export MERQURY={params.export}
export OMP_NUM_THREADS={threads} export OMP_NUM_THREADS={threads}
(bash {params.script} {input.merylDB} {input.assembly} {params.outFile}) &> {log} (bash {params.script} {input.merylDB} {output.symlink} {params.outFile}) &> {log}
""" """
# bash {params.script} {input.merylDB} {params.symlink} {params.outFile} # bash {params.script} {input.merylDB} {params.symlink} {params.outFile}
...@@ -263,11 +268,12 @@ rule merqury: ...@@ -263,11 +268,12 @@ rule merqury:
rule genomescope2: rule genomescope2:
input: input:
hist=os.path.join(config['Results'], "{sample}" +"/2_QVstats_merylAndMerqury/" + "{sample}" + "_dB.21.hist"), hist=os.path.join(config['Results'], "{sample}" +"/2_QVstats_merylAndMerqury/" + "{sample}" + "_dB.21.hist"),
check1=os.path.join(config['Results'],"{sample}" + "/2_QVstats_merylAndMerqury/" + "{sample}" + "_merq.{sample}_assembly.spectra-cn.fl.png"), check1=os.path.join(config['Results'],"{sample}" + "/2_QVstats_merylAndMerqury/" + "{sample}" + "_merq.{sample}.spectra-cn.fl.png"),
check2=os.path.join(config['Results'],"{sample}" + "/2_QVstats_merylAndMerqury/" + "{sample}" + "_merq.{sample}_assembly.spectra-cn.st.png") check2=os.path.join(config['Results'],"{sample}" + "/2_QVstats_merylAndMerqury/" + "{sample}" + "_merq.{sample}.spectra-cn.st.png")
params: params:
outFolder=os.path.join(config['Results'], "{sample}" +"/3_genomescopeProfile/"), outFolder=os.path.join(config['Results'], "{sample}" +"/3_genomescopeProfile/"),
cpHist=os.path.join(config['Results'], "{sample}" +"/3_genomescopeProfile/" + "{sample}" + "_dB.21.histo") cpHist=os.path.join(config['Results'], "{sample}" +"/3_genomescopeProfile/" + "{sample}" + "_dB.21.histo"),
findGscope=os.path.join(workflow.basedir)
output: output:
summary=report(os.path.join(config['Results'], "{sample}" +"/3_genomescopeProfile/" + "{sample}" + "_summary.txt"), caption="../report/genomescope/summary.rst", category="Genomescope Profile", subcategory="{sample}"), summary=report(os.path.join(config['Results'], "{sample}" +"/3_genomescopeProfile/" + "{sample}" + "_summary.txt"), caption="../report/genomescope/summary.rst", category="Genomescope Profile", subcategory="{sample}"),
logPlot=report(os.path.join(config['Results'], "{sample}" +"/3_genomescopeProfile/" + "{sample}" + "_log_plot.png"), caption="../report/genomescope/logplot.rst", category="Genomescope Profile", subcategory="{sample}"), logPlot=report(os.path.join(config['Results'], "{sample}" +"/3_genomescopeProfile/" + "{sample}" + "_log_plot.png"), caption="../report/genomescope/logplot.rst", category="Genomescope Profile", subcategory="{sample}"),
...@@ -278,7 +284,8 @@ rule genomescope2: ...@@ -278,7 +284,8 @@ rule genomescope2:
shell: shell:
""" """
expand -t 1 {input.hist} > {params.cpHist} expand -t 1 {input.hist} > {params.cpHist}
Rscript genomescope2 -i {params.cpHist} -o {params.outFolder} -k 21 -n {wildcards.sample} export genomescope2script=$(find {params.findGscope} -name genomescope2)
Rscript $genomescope2script -i {params.cpHist} -o {params.outFolder} -k 21 -n {wildcards.sample}
grep "Genome Haploid Length" {output.summary} | awk {{'print $6'}} | sed 's/,//g'> {output.estimatedSize} grep "Genome Haploid Length" {output.summary} | awk {{'print $6'}} | sed 's/,//g'> {output.estimatedSize}
""" """
...@@ -342,8 +349,8 @@ rule saveConfiguration_and_getKeyValues: ...@@ -342,8 +349,8 @@ rule saveConfiguration_and_getKeyValues:
gscopeLin=os.path.join(config['Results'], "{sample}" +"/3_genomescopeProfile/" + "{sample}" + "_linear_plot.png"), gscopeLin=os.path.join(config['Results'], "{sample}" +"/3_genomescopeProfile/" + "{sample}" + "_linear_plot.png"),
busco=os.path.join(config['Results'], "{sample}" + "/1_busco5/short_summary.specific." + buscoDataBaseName + "_odb10." + "{sample}" + ".txt"), busco=os.path.join(config['Results'], "{sample}" + "/1_busco5/short_summary.specific." + buscoDataBaseName + "_odb10." + "{sample}" + ".txt"),
qv=os.path.join(config['Results'],"{sample}" + "/2_QVstats_merylAndMerqury/" + "{sample}" + "_merq.qv"), qv=os.path.join(config['Results'],"{sample}" + "/2_QVstats_merylAndMerqury/" + "{sample}" + "_merq.qv"),
spectraStacked=os.path.join(config['Results'],"{sample}" + "/2_QVstats_merylAndMerqury/" + "{sample}" + "_merq.{sample}_assembly.spectra-cn.st.png"), spectraStacked=os.path.join(config['Results'],"{sample}" + "/2_QVstats_merylAndMerqury/" + "{sample}" + "_merq.{sample}.spectra-cn.st.png"),
spectraFlat=os.path.join(config['Results'],"{sample}" + "/2_QVstats_merylAndMerqury/" + "{sample}" + "_merq.{sample}_assembly.spectra-cn.fl.png"), spectraFlat=os.path.join(config['Results'],"{sample}" + "/2_QVstats_merylAndMerqury/" + "{sample}" + "_merq.{sample}.spectra-cn.fl.png"),
# os.path.join(config['Results'],"{sample}" + "/QVstats_merylAndMerqury/logs/" + "{sample}" + "_merqury.log"), # os.path.join(config['Results'],"{sample}" + "/QVstats_merylAndMerqury/logs/" + "{sample}" + "_merqury.log"),
# os.path.join(config['Results'],"{sample}" + "/assemblyStats/{sample}_contig_stats.tsv"), # os.path.join(config['Results'],"{sample}" + "/assemblyStats/{sample}_contig_stats.tsv"),
scaffStats=os.path.join(config['Results'],"{sample}" + "/4_assemblyStats/{sample}_scaffold_stats.tsv"), scaffStats=os.path.join(config['Results'],"{sample}" + "/4_assemblyStats/{sample}_scaffold_stats.tsv"),
...@@ -367,8 +374,8 @@ rule saveConfiguration_and_getKeyValues: ...@@ -367,8 +374,8 @@ rule saveConfiguration_and_getKeyValues:
gscopeLin=os.path.join(config['Results'], "{sample}" +"/5_Key_Results/" + "{sample}" + "_gScope_linear_plot.png"), gscopeLin=os.path.join(config['Results'], "{sample}" +"/5_Key_Results/" + "{sample}" + "_gScope_linear_plot.png"),
busco=os.path.join(config['Results'], "{sample}" + "/5_Key_Results/short_summary.specific." + buscoDataBaseName + "_odb10." + "{sample}" + ".txt"), busco=os.path.join(config['Results'], "{sample}" + "/5_Key_Results/short_summary.specific." + buscoDataBaseName + "_odb10." + "{sample}" + ".txt"),
qv=os.path.join(config['Results'],"{sample}" + "/5_Key_Results/" + "{sample}" + "_merq.qv"), qv=os.path.join(config['Results'],"{sample}" + "/5_Key_Results/" + "{sample}" + "_merq.qv"),
spectraStacked=os.path.join(config['Results'],"{sample}" + "/5_Key_Results/" + "{sample}" + "_merq.{sample}_assembly.spectra-cn.st.png"), spectraStacked=os.path.join(config['Results'],"{sample}" + "/5_Key_Results/" + "{sample}" + "_merq.{sample}.spectra-cn.st.png"),
spectraFlat=os.path.join(config['Results'],"{sample}" + "/5_Key_Results/" + "{sample}" + "_merq.{sample}_assembly.spectra-cn.fl.png"), spectraFlat=os.path.join(config['Results'],"{sample}" + "/5_Key_Results/" + "{sample}" + "_merq.{sample}.spectra-cn.fl.png"),
# os.path.join(config['Results'],"{sample}" + "/QVstats_merylAndMerqury/logs/" + "{sample}" + "_merqury.log"), # os.path.join(config['Results'],"{sample}" + "/QVstats_merylAndMerqury/logs/" + "{sample}" + "_merqury.log"),
# os.path.join(config['Results'],"{sample}" + "/assemblyStats/{sample}_contig_stats.tsv"), # os.path.join(config['Results'],"{sample}" + "/assemblyStats/{sample}_contig_stats.tsv"),
scaffStats=os.path.join(config['Results'],"{sample}" + "/5_Key_Results/{sample}_scaffold_stats.tsv"), scaffStats=os.path.join(config['Results'],"{sample}" + "/5_Key_Results/{sample}_scaffold_stats.tsv"),
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment