Commit 81c86bc7 authored by khalid's avatar khalid
Browse files

Changes to deal with empty samples

parent c05d3708
...@@ -8,10 +8,7 @@ import csv ...@@ -8,10 +8,7 @@ import csv
############# #############
# Wildcards # # Wildcards #
############# #############
wildcard_constraints:
sample="[^\/]+",
id="[^\/]+"
STEPS = config["steps"] STEPS = config["steps"]
PREPARE_REPORT_OUTPUTS = config["prepare_report_outputs"] PREPARE_REPORT_OUTPUTS = config["prepare_report_outputs"]
PREPARE_REPORT_SCRIPTS = config["prepare_report_scripts"] PREPARE_REPORT_SCRIPTS = config["prepare_report_scripts"]
...@@ -19,16 +16,6 @@ OUTPUTS = config["outputs"] ...@@ -19,16 +16,6 @@ OUTPUTS = config["outputs"]
PARAMS_INFO = config["params_info"] PARAMS_INFO = config["params_info"]
config = config["params"] config = config["params"]
##########
# Inputs #
##########
# raw_inputs function call
raw_reads = raw_reads(config['results_dir'], config['sample_dir'], config['SeOrPe'])
config.update(raw_reads)
SAMPLES = raw_reads['samples']
def get_individus(): def get_individus():
if (config["demultiplexing"] == "null"): if (config["demultiplexing"] == "null"):
return SAMPLES return SAMPLES
...@@ -39,6 +26,15 @@ def get_individus(): ...@@ -39,6 +26,15 @@ def get_individus():
individus = get_individus() individus = get_individus()
##########
# Inputs #
##########
# raw_inputs function call
raw_reads = raw_reads(config['results_dir'], config['sample_dir'], config['SeOrPe'])
config.update(raw_reads)
SAMPLES = raw_reads['samples']
## get reads ## get reads
def reads(): def reads():
inputs = dict() inputs = dict()
...@@ -164,7 +160,7 @@ def step_outputs(step): ...@@ -164,7 +160,7 @@ def step_outputs(step):
else: else:
outputs = expand(rules.merge_process_radtags_SE.output,id=individus) outputs = expand(rules.merge_process_radtags_SE.output,id=individus)
elif (config[step] == "null"): elif (config[step] == "null"):
reads = raw_reads reads = raw_reads()
if(config["SeOrPe"] == "PE"): if(config["SeOrPe"] == "PE"):
outputs = list(expand(reads["read"],sample=individus)) outputs = list(expand(reads["read"],sample=individus))
outputs.extend(expand(reads["read2"],sample=individus)) outputs.extend(expand(reads["read2"],sample=individus))
...@@ -379,6 +375,7 @@ rule bowtie_PE: ...@@ -379,6 +375,7 @@ rule bowtie_PE:
best = "--best" if (config["bowtie_best"]) else "", best = "--best" if (config["bowtie_best"]) else "",
strata = "--strata" if (config["bowtie_strata"]) else "" strata = "--strata" if (config["bowtie_strata"]) else ""
shell: shell:
"if [ -s {input.read} ] || [ -s {input.read2} ]; then "
"bowtie " "bowtie "
"--threads {threads} " "--threads {threads} "
"{params.indexPrefix} " "{params.indexPrefix} "
...@@ -391,9 +388,14 @@ rule bowtie_PE: ...@@ -391,9 +388,14 @@ rule bowtie_PE:
"-S "#output in sam format "-S "#output in sam format
"-1 {input.read} " "-1 {input.read} "
"-2 {input.read2} 2> {log} " "-2 {input.read2} 2> {log} "
"| samtools view -b 2>> {log} " "| samtools view -b | samtools addreplacerg -r 'ID:{wildcards.sample}' -r 'SM:{wildcards.sample}' - 2>> {log} "
"| samtools sort -@ {threads} > {output.bam} 2>> {log} &&" "| samtools sort -@ {threads} > {output.bam} 2>> {log} &&"
"samtools index -@ {threads} {output.bam} 2>> {log}" "samtools index -@ {threads} {output.bam} 2>> {log}"
"else "
"touch {output.bam};"
"echo fastq file is empty > {log};"
"fi"
rule bowtie_SE: rule bowtie_SE:
input: input:
...@@ -411,9 +413,10 @@ rule bowtie_SE: ...@@ -411,9 +413,10 @@ rule bowtie_SE:
params: params:
indexPrefix = config["bowtie_index_path"]+"/index", indexPrefix = config["bowtie_index_path"]+"/index",
mult_align_limit = config["bowtie_mult_align_limit"], mult_align_limit = config["bowtie_mult_align_limit"],
best = "--best" if (config["bowtie_best"]) else "", best = config["bowtie_best"],
strata = "--strata" if (config["bowtie_strata"]) else "" strata = config["bowtie_strata"]
shell: shell:
"if [ -s {input.read} ]; then "
"bowtie " "bowtie "
"--threads {threads} " "--threads {threads} "
"{params.indexPrefix} " "{params.indexPrefix} "
...@@ -422,9 +425,13 @@ rule bowtie_SE: ...@@ -422,9 +425,13 @@ rule bowtie_SE:
"{params.strata} " "{params.strata} "
"-S "#output in sam format "-S "#output in sam format
"{input.read} 2> {log} " "{input.read} 2> {log} "
"| samtools view -b 2>> {log} " "| samtools view -b | samtools addreplacerg -r 'ID:{wildcards.sample}' -r 'SM:{wildcards.sample}' - 2>> {log} "
"| samtools sort -@ {threads} > {output.bam} 2>> {log} &&" "| samtools sort -@ {threads} > {output.bam} 2>> {log} &&"
"samtools index -@ {threads} {output.bam} 2>> {log}" "samtools index -@ {threads} {output.bam} 2>> {log}"
"else "
"touch {output.bam};"
"echo fastq file is empty > {log};"
"fi"
ruleorder: bowtie_PE > bowtie_SE ruleorder: bowtie_PE > bowtie_SE
...@@ -468,6 +475,7 @@ rule bwa_mem_PE: ...@@ -468,6 +475,7 @@ rule bwa_mem_PE:
indexPrefix = config["bwa_index_path"]+"/index", indexPrefix = config["bwa_index_path"]+"/index",
quality0_multimapping = """| awk '{{pos = index("XA:Z",$0); if (pos>0) $5=0;print}}'""" if (config["bwa_mem_quality0_multimapping"] == True) else "" quality0_multimapping = """| awk '{{pos = index("XA:Z",$0); if (pos>0) $5=0;print}}'""" if (config["bwa_mem_quality0_multimapping"] == True) else ""
shell: shell:
"if [ -s {input.read} ] || [ -s {input.read2} ]; then "
"bwa mem " "bwa mem "
"-t {threads} " "-t {threads} "
"{params.indexPrefix} " "{params.indexPrefix} "
...@@ -478,6 +486,10 @@ rule bwa_mem_PE: ...@@ -478,6 +486,10 @@ rule bwa_mem_PE:
"| samtools view -b 2>> {log} " "| samtools view -b 2>> {log} "
"| samtools sort -@ {threads} > {output.bam} 2>> {log} &&" "| samtools sort -@ {threads} > {output.bam} 2>> {log} &&"
"samtools index -@ {threads} {output.bam} 2>> {log}" "samtools index -@ {threads} {output.bam} 2>> {log}"
"else "
"touch {output.bam};"
"echo fastq file is empty > {log};"
"fi"
rule bwa_mem_SE: rule bwa_mem_SE:
input: input:
...@@ -499,6 +511,7 @@ rule bwa_mem_SE: ...@@ -499,6 +511,7 @@ rule bwa_mem_SE:
indexPrefix = config["bwa_index_path"]+"/index", indexPrefix = config["bwa_index_path"]+"/index",
quality0_multimapping = """| awk '{{pos = index("XA:Z",$0); if (pos>0) $5=0;print}}'""" if (config["bwa_mem_quality0_multimapping"] == True) else "" quality0_multimapping = """| awk '{{pos = index("XA:Z",$0); if (pos>0) $5=0;print}}'""" if (config["bwa_mem_quality0_multimapping"] == True) else ""
shell: shell:
"if [ -s {input.read} ] ; then "
"bwa mem " "bwa mem "
"-t {threads} " "-t {threads} "
"{params.indexPrefix} " "{params.indexPrefix} "
...@@ -508,6 +521,10 @@ rule bwa_mem_SE: ...@@ -508,6 +521,10 @@ rule bwa_mem_SE:
"| samtools view -b 2>> {log} " "| samtools view -b 2>> {log} "
"| samtools sort -@ {threads} > {output.bam} 2>> {log} &&" "| samtools sort -@ {threads} > {output.bam} 2>> {log} &&"
"samtools index -@ {threads} {output.bam} 2>> {log}" "samtools index -@ {threads} {output.bam} 2>> {log}"
"else "
"touch {output.bam};"
"echo fastq file is empty > {log};"
"fi"
ruleorder: bwa_mem_PE > bwa_mem_SE ruleorder: bwa_mem_PE > bwa_mem_SE
...@@ -545,18 +562,37 @@ rule gstacks: ...@@ -545,18 +562,37 @@ rule gstacks:
config["gstacks_threads"] config["gstacks_threads"]
log: log:
config["results_dir"]+"/logs/gstacks/gstacks_log.txt" config["results_dir"]+"/logs/gstacks/gstacks_log.txt"
shell: shell: """
"gstacks " nonemptySamples='';
"-t {threads} " emptycount=0
"-I {params.bam_dir} " for i in `ls {params.bam_dir}/*.bam`; do [ -s $i ] || emptycount=$((emptycount + 1)) ; done;
"-M {input.popmap} " if [[ ${{emptycount}} -gt 0 ]];
"--model {params.model} " then
"--var-alpha {params.var_alpha} " # here the bam file must include RG and SM
"--gt-alpha {params.gt_alpha} " for i in `ls {params.bam_dir}/*.bam`;
"--min-mapq {params.min_mapq} " do
"--max-clipped {params.max_clipped} " if [[ -s $i ]];
"-O {params.output_dir} " then
"|& tee {log}" countRG=`samtools view -H $i | grep -c @RG`
if [[ ${{countRG}} -eq 0 ]];
then
sample=`basename $i .bam`;
samtools addreplacerg -r ID:$sample -r SM:$sample -o tmp.bam $i
mv tmp.bam $i
fi
nonemptySamples=${{nonemptySamples}}" -B "$i ;
fi
done;
gstacks -t {threads} $nonemptySamples --model {params.model} \
--var-alpha {params.var_alpha} --gt-alpha {params.gt_alpha} --min-mapq {params.min_mapq} \
--max-clipped {params.max_clipped} -O {params.output_dir} |& tee {log}
else
gstacks -t {threads} -I {params.bam_dir} -M {input.popmap} --model {params.model} \
--var-alpha {params.var_alpha} --gt-alpha {params.gt_alpha} --min-mapq {params.min_mapq} \
--max-clipped {params.max_clipped} -O {params.output_dir} |& tee {log}
fi
"""
rule populations: rule populations:
...@@ -594,7 +630,7 @@ rule populations: ...@@ -594,7 +630,7 @@ rule populations:
"--plink " "--plink "
"--phylip " "--phylip "
"|& tee {log} " "|& tee {log} "
"&& gzip {params.output_dir}/populations.snps.vcf" "&& gzip -f {params.output_dir}/populations.snps.vcf"
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment