Allow genome size to be configurable in downsample_amplicons script (…

…#51)

bin/downsample_amplicons.py

@@ -229,7 +229,7 @@ def main(args):
 
     required_depth_achieved = [False] * len(genome_checkpoints)
 
-    depths = [0] * 30000
+    depths = [0] * (args.genome_size + 100)
 
     infile = pysam.AlignmentFile(args.bam, "rb")
     bam_header = infile.header.copy().to_dict()
@@ -303,6 +303,7 @@ def main(args):
     parser.add_argument('--min-depth', type=int, default=200, help='Subsample to n coverage')
     parser.add_argument('--mapping-quality', type=int, default=20, help='Minimum mapping quality to include read in output')
     parser.add_argument('--amplicon-subdivisions', type=int, default=3, help='How many times to divide amplicons to detect coverage')
+    parser.add_argument('--genome-size', type=int, default=29903, help='')
     parser.add_argument('--verbose', action='store_true', help='Debug mode')
     args = parser.parse_args()
     main(args)

modules/illumina.nf

@@ -228,7 +228,7 @@ process callConsensusFreebayes {
         bcftools consensus -f ${ref} -I ${sampleName}.ambiguous.norm.vcf.gz > ${sampleName}.ambiguous.fa
 
         # apply remaninng variants, including indels
-        bcftools consensus -f ${sampleName}.ambiguous.fa -m ${sampleName}.mask.txt ${sampleName}.fixed.norm.vcf.gz | sed s/MN908947.3/${sampleName}/ > ${sampleName}.consensus.fa
+        bcftools consensus -f ${sampleName}.ambiguous.fa -m ${sampleName}.mask.txt ${sampleName}.fixed.norm.vcf.gz | sed s/${params.viral_contig_name}/${sampleName}/ > ${sampleName}.consensus.fa
         """
 }
 

modules/utils.nf

@@ -27,14 +27,21 @@ process downsampleAmplicons {
     publishDir "${params.outdir}/${task.process.replaceAll(":","_")}", pattern: "${sampleId}.mapped.primertrimmed.downsampled.sorted{.bam,.bam.bai}", mode: 'copy'
 
     input:
-        tuple val(sampleId), path(trimmed_bam), path(trimmed_bam_index), path(bedfile)
+        tuple val(sampleId), path(trimmed_bam), path(trimmed_bam_index), path(bedfile), path(ref)
     output:
         tuple val(sampleId), path("${sampleId}.mapped.primertrimmed.downsampled.sorted.bam"), path("${sampleId}.mapped.primertrimmed.downsampled.sorted.bam.bai"), emit: alignment
 	path("downsampling_summary.csv"), emit: summary
 
     script:
         """
-        downsample_amplicons.py --bed ${bedfile} --min-depth ${params.downsampleMinDepth} --mapping-quality ${params.downsampleMappingQuality} --amplicon-subdivisions ${params.downsampleAmpliconSubdivisions}  ${trimmed_bam} 2> downsampling_summary_no_sample_id.csv | \
+        samtools faidx ${ref}
+        downsample_amplicons.py \
+          --bed ${bedfile} \
+          --min-depth ${params.downsampleMinDepth} \
+          --mapping-quality ${params.downsampleMappingQuality} \
+          --amplicon-subdivisions ${params.downsampleAmpliconSubdivisions} \
+          --genome-size \$(cut -d \$'\\t' -f 2 ${ref}.fai) \
+          ${trimmed_bam} 2> downsampling_summary_no_sample_id.csv | \
             samtools sort - -o ${sampleId}.mapped.primertrimmed.downsampled.sorted.bam
         samtools index ${sampleId}.mapped.primertrimmed.downsampled.sorted.bam
         paste -d ',' <(echo "sample_id" && echo "${sampleId}") downsampling_summary_no_sample_id.csv > downsampling_summary.csv

workflows/illuminaNcov.nf

@@ -108,7 +108,7 @@ workflow sequenceAnalysis {
 
       trimPrimerSequences(readMapping.out.combine(ch_bedFile))
 
-      downsampleAmplicons(trimPrimerSequences.out.ptrim.combine(ch_bedFile))
+      downsampleAmplicons(trimPrimerSequences.out.ptrim.combine(ch_bedFile).combine(Channel.fromPath(params.ref)))
 
       downsampledBamToFastq(downsampleAmplicons.out.alignment)