startWorkflow.nf

#! /usr/bin/env nextflow


switch(params.platform) {
	case "i":
	case "Illumina":
		print("Paired end Illumina mode")
		platform = "Illumina"
		isPairedEnd = true
		break;
	case "s":
		print("Single read Illumina mode")
		platform = "Illumina"
		isPairedEnd = false
		break;
	case "n":
	case "ONT":
	case "nanopore":
		print("Nanopore mode")
		platform = "ONT"
		isPairedEnd = false
		break;
	case "p":
	case "pb":
	case "PB":
	case "pacbio":
	case "PacBio":
		print("PacBio mode")
		platform = "PacBio"
		isPairedEnd = false
		break;
	case "h":
		printUsage()
		return 0
	default:
		print("ERROR: unknown sequencing modality")
		printUsage()
		return 1
}


// Convert relative path to absolute path
if (params.primerBedFile[0]=='/' || params.primerBedFile[0]=='~')
	primerBedFile=params.primerBedFile
else
	primerBedFile="$launchDir/$params.primerBedFile"


def printUsage () {
	println("Example usage: ")
	println("\t./startWorkflow.nf --platform i --primers path/to/bed --in path/to/fastq/ --out path/to/outputDir")
	println()
	println("Run the C-WAP workflow that will import the fastq files, trim adaptors, apply quality trimming, call variants and generate a html and pdf report and quit. All fastq files in the provided directory will be processed.")
	println()
	
	println("Platform options:")
	println("\t-h\t\t: Show this help message and exit.")
	println("\t-i\t\t: Use parameters optimised for the Illumina platform, paired-end mode")
	println("\t-s\t\t: Use parameters optimised for the Illumina platform, single read mode")
	println("\t-n\t\t: Use parameters optimised for the ONT platform")
	println("\t-p\t\t: Use parameters optimised for PacBio platform")
}


// Given a file name, extracts a human readable sample name to be used in the output.
allSampleNames = []
def getSampleName(filename) {
	// Typical for most users using Illumina output as is
	// Ex: /path/to/dir/something_S1_L2_R1.fastq -> something
	sampleName = filename.name.split("/")[-1].split("\\.")[0].split("_")[0]
    	
	// For complicated file names involving underscores that cannot be eliminated
	// Ex: /path/to/dir/some_thing_R1.fastq -> some-thing	
	// sampleName = filename.name.split("/")[-1].split("\\.")[0].split("_R")[0].replace('_','-')
	
	// If there are special fixed substrings available within all file names.
	// Ex: /path/to/dir/some_thing_1art_out_R1.fastq -> some-thing-1
	// sampleName = filename.name.split("/")[-1].split("\\.")[0].split("art_out")[0].replace('_','-')
    
    // PacBio Lima output pattern
    // some-name.bc1234.fastq.gz
    // sampleName = filename.name.split("/")[-1].split("\\.")[1].split("--")[0]

	
	// Check if there is name collision and immediately abort the execution if this is the case.
	if (allSampleNames.contains(sampleName)) {
		throw new IllegalArgumentException ('The sample name is not unique, please adjust the getSampleName function.')
	}
	else {
		allSampleNames += sampleName
		return sampleName
	}
}


// Import the list of files to process
if (isPairedEnd) {
	FQs = Channel
	    .fromFilePairs("$params.in/*_R{1,2}*.fastq*", checkIfExists: true, flat:true)
		.map{ tuple(getSampleName(it[1]), it[1], it[2]) }
}
else {
	FQs = Channel
		.fromPath( "$params.in/*.fastq*", checkIfExists: true)
		.map{ tuple(getSampleName(it), it, null) }
}


FQs
	.view()
	.into{ input_fq_a; input_fq_b; input_fq_c; input_fq_d; input_fq_e }



//////////////////////////////////////////////////////////////////////////
// Read processing using common tools, QC etc.
//////////////////////////////////////////////////////////////////////////

// Align the reads to the reference sequence to obtain a sorted bam file
process referenceAlignment {
	label 'high_cpu'
	
	input:
		tuple val(sampleName), file('R1.fastq.gz'), file('R2.fastq.gz') from input_fq_a
	
	output:
		tuple val(sampleName), file('aligned.sam') into aligned_sam

	if (platform == "Illumina")
		conda "$projectDir/conda/env-bowtie2"
	else
		conda "$projectDir/conda/env-minimap2"
	
	shell:
	refSeqBasename = params.referenceSequence.replaceAll('.fa$', '')
	"""
		if [[ -n \$SLURM_CPUS_ON_NODE ]]; then
			numThreads=\$SLURM_CPUS_ON_NODE
		else
			numThreads=\$(nproc)
		fi
		
		case $platform in
			Illumina)
				if $isPairedEnd; then
					bowtie2 --no-unal --threads \$numThreads -x $refSeqBasename -1 R1.fastq.gz -2 R2.fastq.gz \
						-S aligned.sam
				else
					bowtie2 --no-unal --threads \$numThreads -x $refSeqBasename -U R1.fastq.gz -S aligned.sam
				fi
				;;
			ONT)
				minimap2 -a --sam-hit-only -2 -x map-ont ${refSeqBasename}.mmi R1.fastq.gz \
						-t \$numThreads -o aligned.sam
				;;
			PacBio)
				minimap2 -a --sam-hit-only -2 -x map-hifi ${refSeqBasename}.mmi R1.fastq.gz \
						-t \$numThreads -o aligned.sam
				;;
		esac
	"""
}


// Align the reads to the reference sequence to obtain a sorted bam file
process trimming {
	label 'high_cpu'
	
	input:
		tuple val(sampleName), file('aligned.sam') from aligned_sam
	
	output:
		tuple val(sampleName), env(numReads), path('resorted.bam') into resorted_bam_a, resorted_bam_b,  resorted_bam_c,  resorted_bam_d
		tuple val(sampleName), path('sorted.stats'), path('resorted.stats'), path('primer_hit_counts.tsv') into samtools_stats
	
	// An older version of samtools is automatically provided by ivar's dependency
    // Does one need to provide samtools=1.15 here? 
	conda "$projectDir/conda/env-ivar"

	shell:
	"""
		if [[ -n \$SLURM_CPUS_ON_NODE ]]; then
			numThreads=\$SLURM_CPUS_ON_NODE
		else
			numThreads=\$(nproc)
		fi
		
		# Generation of a sorted bam file from the alignment output
		samtools sort aligned.sam -o sorted.bam -@ \$numThreads
		
		# Nanopore has a much lower read quality, so the quality trimming should be much more lax.
		if [[ $platform == ONT ]]; then
			ivar trim -e -b $primerBedFile -p trimmed -i sorted.bam -q 1 | tee ivar.stdout
		else
			ivar trim -e -b $primerBedFile -p trimmed -i sorted.bam | tee ivar.stdout
		fi
		
        # Generate a tsv file tabulating the number of reads vs trimmer primer name in the bed file
        cat ivar.stdout | grep -A 10000 "Primer Name" | head -n -5 > primer_hit_counts.tsv
        
		samtools sort trimmed.bam -o resorted.bam -@ \$numThreads
		
		# Evaluate read statistics
		samtools stats sorted.bam | grep ^SN | cut -f 2- > sorted.stats
		samtools stats resorted.bam | grep ^SN | cut -f 2- > resorted.stats
		numReads=\$(cat resorted.stats | grep "raw total sequences" | awk '{ print \$4 }')
        
   		rm sorted.bam trimmed.bam
	"""
}



process trimmedBam2Fastq {
	input:
		tuple val(sampleName), env(numReads), path('resorted.bam') from resorted_bam_a
	
	output:
		tuple val(sampleName), env(numReads), path('resorted.fastq.gz') into resorted_fastq_gz_a, resorted_fastq_gz_b, resorted_fastq_gz_c
	
	conda "$projectDir/conda/env-samtools"
	
	shell:
	"""
		if [[ \$numReads -gt 1000000 ]]; then
			# A maximum of 1 000 000 reads are kept to limit the computation time of variant calling processes.
			echo Subsampling the dataset...
			retentionRatio=\$(echo 3k 1000000 \$numReads /p | dc)
			samtools view --threads 2 --subsample \$retentionRatio resorted.bam -o resorted.sam
		else
			samtools view --threads 2 resorted.bam -o resorted.sam
		fi
		
		$projectDir/sam2fastq.py resorted.sam resorted.fastq
		gzip resorted.fastq
		rm resorted.sam
	"""
}


process generatePileup {
	input:
		tuple val(sampleName), env(numReads), path('resorted.bam') from resorted_bam_b
	
	output:
		tuple val(sampleName), path('pile.up') into pile_up_a, pile_up_b
	
	conda "$projectDir/conda/env-samtools"
	
	shell:
	"""
		samtools mpileup -aa -A -d 10000 -B -Q 0 --reference $params.referenceSequence -o pile.up resorted.bam
	"""
}



process variantCalling {
	input:
		tuple val(sampleName), path('pile.up') from pile_up_a
	
	output:
		tuple val(sampleName), path('rawVarCalls.tsv') into ivar_out
	
	conda "$projectDir/conda/env-ivar"
	
	shell:
	"""
		cat pile.up | ivar variants -p rawVarCalls -g $projectDir/covidRefSequences/covidGenomeAnnotation-NCBI.gff \
		-r $params.referenceSequence -m 10
	"""
}


// Kraken2-based taxonomic classification
process k2stdDB {
	memory '70 GB'
	label 'high_IO'
	
	input:
		tuple val(sampleName), file('R1.fastq.gz'), file('R2.fastq.gz') from input_fq_b
	
	output:
		tuple val(sampleName), path('k2-std.out') into k2_std_out
	
	conda "$projectDir/conda/env-kraken2"
	
	shell:
	"""	
		if $isPairedEnd; then
			kraken2 --paired R1.fastq.gz R2.fastq.gz --db \$K2_STD_DB_PATH --threads 2 --report k2-std.out > /dev/null
		else
			kraken2 R1.fastq.gz --db \$K2_STD_DB_PATH --threads 2 --report k2-std.out > /dev/null
		fi	
	"""
}


// The pileup file is parsed to calculate the positionwise quality and depth parameters.
// The result is stored as png files that are added to the html report
process QCplots {
	input:
		tuple val(sampleName), path('pile.up') from pile_up_b
	
	output:
		tuple val(sampleName), path('pos-coverage-quality.tsv'), path('coverage.png'), path('depthHistogram.png'), path('quality.png'), path('qualityHistogram.png'), path('discontinuitySignal.png'), path('genesVSuncovered_abscounts.png'), path('genesVSuncovered_scaled.png'), path('breadthVSdepth.png') into QChists
	
	conda "$projectDir/conda/env-python"
	
	shell:
	"""
		python3 $projectDir/plotQC.py pile.up $primerBedFile
	"""
}


// Draw a histogram of all read lengths
process readLengthHist {
	input:
		tuple val(sampleName), file('R1.fastq.gz'), file('R2.fastq.gz') from input_fq_d
	
	output:
		tuple val(sampleName), path('readLengthHist.png') into readLengthHist_png
        tuple val(sampleName), path('readLengthHist.png'), path('timeVSreadcounts.png') into readLengthHist_out


	conda "$projectDir/conda/env-python"
	label 'high_IO'

	shell:
	"""
		gzip -dc R1.fastq.gz > allreads.fastq
		if $isPairedEnd; then
			gzip -dc R2.fastq.gz >> allreads.fastq
		fi
		
		# Only up to 1 million reads will be considered for the length histogram
		head -n 4000000 allreads.fastq | awk 'NR%4==2' | awk "{print length}" | python3 $projectDir/plotLengthHist.py
        
        # If this is an ONT run with detailed headers, also plot read count w.r.t time to show if the data collection
        # was essentially complete.
        # An example ONT fastq entry header with a timestamp looks like this:
        # @b6a37669-02b7-4d48-bcff-ad7e2eb4fa06 runid=52014692ae58b6d24cb1dcd29fd35d118e5f6a42 read=16 ch=304 start_time=2022-06-29T15:57:00.583046-04:00 flow_cell_id=FAT09511 protocol_group_id=VSS_spikein_June22toJune28_062922 sample_id=no_sample barcode=barcode49 barcode_alias=barcode49 parent_read_id=b6a37669-02b7-4d48-bcff-ad7e2eb4fa06 basecall_model_version_id=2021-05-17_dna_r9.4.1_minion_96_29d8704b
   		touch timeVSreadcounts.png
        if [[ $platform == ONT ]]; then
            cat allreads.fastq | grep start_time | cut -d ' ' -f 5 | awk -F 'start_time=' '{print \$2}' > timestamps
            if [[ -s timestamps ]]; then
                $projectDir/plotTimeVSreadcounts.py ./timestamps timeVSreadcounts.png
            fi
        fi
        
   		rm allreads.fastq
	"""
}



// ////////////////////////////////////////////
// VARIANT CALLING
// ////////////////////////////////////////////

process krakenVariantCaller {
	input:
		tuple val(sampleName), env(numReads), path('resorted.fastq.gz') from resorted_fastq_gz_a
	
	output:
		tuple val(sampleName), path('k2-allCovid_bracken*.out'), path('k2-majorCovid_bracken*.out'), path('k2-allCovid.out'), path('k2-majorCovid.out') into k2_covid_out
	
	conda "$projectDir/conda/env-kraken2"
    
	shell:
	"""
		if [[ -n \$SLURM_CPUS_ON_NODE ]]; then
			numThreads=\$SLURM_CPUS_ON_NODE
		else
			numThreads=\$(nproc)
		fi

		# Check the number of reads. Ignore if there are too few reads
		kraken2 resorted.fastq.gz --db $projectDir/customDBs/allCovidDB --threads \$numThreads --report k2-allCovid.out > /dev/null
        num_phylum_hits=\$(cat k2-allCovid.out | grep -w P | cut -f 2 | head -n 1)
        if [[ -z \${num_phylum_hits} || \${num_phylum_hits} -lt 10 ]]; then
			# There is a bug in bracken that throws a "no reads found" error, if there
            # are fewer hits in all target level taxons. We check that there are a minimum of 10 reads.
			echo 100.00\$'\t'0\$'\t'0\$'\t'R\$'\t'1\$'\t'root > k2-allCovid_bracken.out
		else
			bracken -d $projectDir/customDBs/allCovidDB -i k2-allCovid.out -o allCovid.bracken -l P
		fi

		kraken2 resorted.fastq.gz --db $projectDir/customDBs/majorCovidDB --threads \$numThreads --report k2-majorCovid.out > /dev/null
        num_class_hits=\$(cat k2-majorCovid.out | grep -w C | cut -f 2 | head -n 1)
        if [[ -z \${num_class_hits} || \${num_class_hits} -lt 10 ]]; then
			# There is a bug in bracken that throws a "no reads found" error, if there
            # are fewer hits in all target level taxons. We check that there are a minimum of 10 reads.
            echo 100.00\$'\t'0\$'\t'0\$'\t'R\$'\t'1\$'\t'root > k2-majorCovid_bracken.out
		else
			bracken -d $projectDir/customDBs/majorCovidDB -i k2-majorCovid.out -o majorCovid.bracken -l C
        fi
	"""
}


// Calculation of the consensus sequence using bcftools
process consensusSequence {
	input:
		tuple val(sampleName), env(numReads), path('resorted.bam') from resorted_bam_c
	
	output:
		tuple val(sampleName), path('consensus.fa') into consensus_fa
	
	conda "$projectDir/conda/env-bcftools"
	
	shell:
	"""		
		bcftools mpileup -d 10000 -Ou -f $params.referenceSequence resorted.bam | bcftools call --ploidy 1 -mv -Oz -o calls.vcf.gz
		bcftools index calls.vcf.gz
		cat $params.referenceSequence | bcftools consensus calls.vcf.gz > consensus.fa
	"""
}



// Characterisation of the consensus sequence based on Pangolin output
// Calculation of the consensus sequence is used to determine the predominant lineage.
process pangolinVariantCaller {
	input:
		tuple val(sampleName), path('consensus.fa') from consensus_fa
	
	output:
		tuple val(sampleName), env(consensusLineage), path('lineage_report.csv'), path('consensus.fa') into pangolin_out
	
	conda "$projectDir/conda/env-pangolin"
	
	shell:
	"""
		pangolin --alignment consensus.fa --threads 2 --outdir ./
		
		# Check the pangolin result to ensure the failure cases are handled properly
		# In non-convergence case, it either throws multiple options or throws "None"
		consensusLineage=\$(tail -n 1 lineage_report.csv | awk -F "," '{ print \$3 }')
		if [[ \$(echo \$consensusLineage | wc -w) -lt 1 || \$consensusLineage == "None" ]]; then
			consensusLineage=Unknown
		fi
	"""
}



process linearDeconVariantCaller {
	input:
		tuple val(sampleName), path('rawVarCalls.tsv') from ivar_out
	
	output:
		tuple val(sampleName), path('linearDeconvolution_abundance.csv'), path('mutationTable.html'), path('VOC-VOIsupportTable.html'), env(mostAbundantVariantPct), env(mostAbundantVariantName), env(linRegressionR2) into linearDeconvolution_out

	conda "$projectDir/conda/env-python"
		
	shell:
	"""
		deconvolutionOutput=\$($projectDir/deconvolveVariants.py rawVarCalls.tsv ./ $params.variantDBfile)
		
		mostAbundantVariantPct=\$(echo \$deconvolutionOutput | awk '{ print \$1 }')
		mostAbundantVariantName=\$(echo \$deconvolutionOutput |awk '{ print \$2 }')
		linRegressionR2=\$(echo \$deconvolutionOutput | awk '{ print \$3 }')
	"""
}


process kallistoVariantCaller {
	input:
		tuple val(sampleName), env(numReads), path('resorted.fastq.gz') from resorted_fastq_gz_b

	output:
		tuple val(sampleName), path('abundance.tsv') into kallisto_out
		
	conda "$projectDir/conda/env-kallisto"
		
	shell:
	"""
		# Check the number of reads. Ignore if there are too few reads
		if [[ $task.attempt -lt 2 ]] && [[ \$numReads -gt 10 ]]; then
			kallisto quant --index $projectDir/customDBs/variants.kalIdx --output-dir ./ \
					--plaintext --threads 2 --single -l 300 -s 50 resorted.fastq.gz
		else
			echo target_id\$'\t'length\$'\t'eff_length\$'\t'est_counts tpm > abundance.tsv
			echo Error\$'\t'29903\$'\t'29903\$'\t'100\$'\t'100 >> abundance.tsv
		fi
	"""
}


// https://github.com/rvalieris/LCS
process LCSvariantCaller {
	input:
		tuple val(sampleName), env(numReads), path('resorted.fastq.gz') from resorted_fastq_gz_c 
	
	output:
		tuple val(sampleName), path('LCS/outputs/decompose/lcs.out') into lcs_out
		
	conda "$projectDir/conda/env-LCS"
	time = '5 min'
	
	shell:
	"""
		if [[ $task.attempt -lt 2 ]] && [[ \$numReads -gt 10 ]]; then
			echo Fetching the LCS repository...
			# git clone https://github.com/rvalieris/LCS.git
			# rm .git -rf
			cp -r $projectDir/LCS ./
			cd LCS
			
			echo Preparing the DB...
			mkdir -p outputs/variants_table
			zcat data/pre-generated-marker-tables/pango-designation-markers-v1.2.9.tsv.gz > outputs/variants_table/pango-markers-table.tsv
			
			echo Preparing the sample dataset...
			mkdir data/fastq
			mv ../resorted.fastq.gz data/fastq/resorted.fastq.gz
			echo "resorted" > data/tags_pool_lcs
			
			echo Executing LCS...
			snakemake --config markers=pango dataset=lcs --cores 2
		else
			echo Not enough covid reads for LCS, skipped.
			mkdir -p LCS/outputs/decompose
			echo sample\$'\t'variant_group\$'\t'proportion\$'\t'mean\$'\t'std_error > LCS/outputs/decompose/lcs.out
			echo ERROR\$'\t'Error\$'\t'1\$'\t'1\$'\t'1 >> LCS/outputs/decompose/lcs.out
		fi
	"""
}


process freyjaVariantCaller {
	label 'high_cpu'

	input:
		tuple val(sampleName), env(numReads), path('resorted.bam') from resorted_bam_d
		
	output:
		tuple val(sampleName), path('freyja.demix'), path('freyja_boot_lineages.csv'), path('freyja_bootstrap.png') into freyja_out

	conda "$projectDir/conda/env-freyja"
	
	shell:
	"""
        if [[ \$numReads -lt 100 ]]; then
            echo INSUFFICIENT DATA > freyja.demix
            echo summarized\$'\t'"[('Undetermined', 1.00)]" >> freyja.demix
            echo lineages\$'\t'Undetermined >> freyja.demix
            echo abundances\$'\t'1.00 >> freyja.demix
            echo resid\$'\t'-1 >> freyja.demix
            echo coverage\$'\t'-1 >> freyja.demix
            
            echo "Undetermined" > freyja_boot_lineages.csv 
            touch freyja_bootstrap.png
        else 
            if [[ $task.attempt -lt 4 ]]; then
                echo Pileup generation for Freyja...
                freyja variants resorted.bam --variants freyja.variants.tsv --depths freyja.depths.tsv --ref $params.referenceSequence
                
                echo Demixing variants by Freyja and bootstrapping
                freyja demix freyja.variants.tsv freyja.depths.tsv --output freyja.demix --confirmedonly &
                freyja boot freyja.variants.tsv freyja.depths.tsv --nt \$(nproc) --nb 1000 --output_base freyja_boot
                wait
                
                echo Parsing bootstrapping output...
                export PYTHONHASHSEED=0
                python3 $projectDir/parseFreyjaBootstraps.py freyja.demix freyja_boot_lineages.csv freyja_bootstrap.png
            else
                # Due to a potential bug, some big fastqs result in a pandas error.
                # Generate an empty file to circumvent such failure cases
                echo FATAL ERROR > freyja.demix
                echo summarized\$'\t'"[('Error', 1.00)]" >> freyja.demix
                echo lineages\$'\t'Error >> freyja.demix
                echo abundances\$'\t'1.00 >> freyja.demix
                echo resid\$'\t'-1 >> freyja.demix
                echo coverage\$'\t'-1 >> freyja.demix
                
                echo "ERROR" > freyja_boot_lineages.csv 
                touch freyja_bootstrap.png
            fi
        fi
	"""
}


// A metadata fetch attemp from NCBI via Entrez-Direct
// Only works if file names explicitly carries an SRR number.
process getNCBImetadata {
	// Allow access to this section only 1 thread at a time to avoid network congestion.
	executor = 'local'
	// queueSize = 1
	maxForks 1
	submitRateLimit = '3/1min'
	
	// NCBI bandwidth limit might cause lookup failures. If so, the next attempt should start with a time delay.
	// Wait some random time so that threads go out of sync.
	errorStrategy { sleep(1000 * Math.random() as long); return 'retry' }
	maxRetries = 105
	
	input:
		tuple val(sampleName), file('R1.fastq.gz'), file('R2.fastq.gz') from input_fq_c
		
	output:
		tuple val(sampleName), env(libraryProtocol), env(seqInstrument), env(isolate), env(collectionDate), env(collectedBy), env(sequencedBy), env(sampleLatitude), env(sampleLongitude), env(sampleLocation) into metadata
	
	conda "$projectDir/conda/env-entrez-direct"
	
	shell:
	"""
	srrNumber=$sampleName
	if [[ $task.attempt -lt 100 ]] && [[ \${srrNumber:0:3} == 'SRR' ]]; then
		# The tool returns error: too many requests, bypassing by redirection of error		
		sraQueryResult=\$(esearch -db sra -query \$srrNumber 2>/dev/null)
		sleep 1
		
		if echo \$sraQueryResult | grep -q "<Count>1</Count>"; then
			# Get runinfo from SRA
			echo Downloading metadata for \$srrNumber...	
			SRRmetadata=\$(echo "\$sraQueryResult" | efetch -format runinfo 2>/dev/null | grep \$srrNumber)
			
			echo Parsing...
			libraryProtocol=\$(echo \$SRRmetadata | awk -F ',' '{print \$13}')
			seqInstrument=\$(echo \$SRRmetadata | awk -F ',' '{print \$20}')
			isolate=\$(echo \$SRRmetadata  | awk -F ',' '{print \$30}')	
			
			# Get metadata out of biosample db
			echo Fetching biosample data...
			SAMN=\$(echo \$SRRmetadata | awk -F ',' '{print \$26}')
			SAMNmetadata=\$(efetch -db biosample -id \$SAMN 2>/dev/null)
			
			echo Parsing...
			collectionDate=\$(echo "\$SAMNmetadata" | grep "collection date" | awk -F '"' '{print \$2}')
			collectedBy=\$(echo "\$SAMNmetadata" | grep "collected by" | awk -F '"' '{print \$2}')
			sequencedBy=\$(echo "\$SAMNmetadata" | grep SEQUENCED_BY | awk '{ \$1=""; print \$0 }')
			sampleLatitude=\$(echo "\$SAMNmetadata" | grep "latitude and longitude" | awk -F '"' '{print \$2}'\
								| awk '{ print \$1\$2 }')
			sampleLongitude=\$(echo "\$SAMNmetadata" | grep "latitude and longitude" | awk -F '"' '{print \$2}'\
								| awk '{ print \$3\$4 }')
			sampleLocation=\$(echo "\$SAMNmetadata" | grep "geographic location" | awk -F '"' '{print \$2}')
		fi
	fi
	
	if [[ -z \$SAMN ]]; then
		SAMN=Missing
	fi

	if [[ -z \$libraryProtocol ]]; then
		libraryProtocol=Missing
	fi

	if [[ -z \$seqInstrument ]]; then
		seqInstrument=Missing
	fi

	if [[ -z \$isolate ]]; then
		isolate=Missing
	fi

	if [[ -z \$collectionDate ]]; then
		collectionDate=Missing
	fi

	if [[ -z \$collectedBy ]]; then
		collectedBy=Missing
	fi

	if [[ -z \$sequencedBy ]]; then
		sequencedBy=Missing
	fi

	if [[ -z \$sampleLatitude ]]; then
		sampleLatitude="?"
	fi

	if [[ -z \$sampleLongitude ]]; then
		sampleLongitude="?"
	fi

	if [[ -z \$sampleLocation ]]; then
		sampleLocation=Missing
	fi
	"""
}



// Computation is now mostly over. All threads need to synchronise here.
// We will group based on the sample name and pass everything to the report
// generation steps.
reportInputCh = metadata.join(samtools_stats).join(k2_std_out).join(QChists)
                        .join(readLengthHist_out).join(linearDeconvolution_out)
                        .join(k2_covid_out).join(pangolin_out).join(kallisto_out)
                        .join(freyja_out).join(lcs_out)



///////////////////////////////////////////////
// Report generation and final output
///////////////////////////////////////////////

// Generates a report based on the computation results generated by executeAnalysis.sh
// One separate html report per each sample (i.e. per fastq/fastq pair)
process generateReport {
	input:
		tuple val(sampleName), env(libraryProtocol), env(seqInstrument), env(isolate), env(collectionDate), env(collectedBy), env(sequencedBy), env(sampleLatitude), env(sampleLongitude), env(sampleLocation),
		path('sorted.stats'), path('resorted.stats'), path('primer_hit_counts.tsv'),
		path('k2-std.out'),
		path('pos-coverage-quality.tsv'), path('coverage.png'), path('depthHistogram.png'), path('quality.png'), path('qualityHistogram.png'), path('discontinuitySignal.png'), path('genesVSuncovered_abscounts.png'), path('genesVSuncovered_scaled.png'), path('breadthVSdepth.png'),
		path('readLengthHist.png'), path('timeVSreadcounts.png'),
		path('linearDeconvolution_abundance.csv'), path('mutationTable.html'), path('VOC-VOIsupportTable.html'), env(mostAbundantVariantPct), env(mostAbundantVariantName), env(linRegressionR2),
		path('k2-allCovid_bracken.out'), path('k2-majorCovid_bracken.out'), path('k2-allCovid.out'), path('k2-majorCovid.out'),
		env(consensusLineage), path('lineage_report.csv'), path('consensus.fa'),
		path('kallisto_abundance.tsv'),
		path('freyja.demix'), path('freyja_boot_lineages.csv'), path('freyja_bootstrap.png'),
		path('lcs.out') from reportInputCh
	
	output:
		file "outfolder" into reportCh
		
	conda "$projectDir/conda/env-python"
	
	shell:
	"""
		echo Making pie charts...
		export PYTHONHASHSEED=0
		$projectDir/plotPieChartsforAbundance.py ./ $params.variantDBfile linearDeconvolution_abundance.csv \
				kallisto_abundance.tsv k2-allCovid_bracken.out k2-majorCovid_bracken.out freyja.demix lcs.out
		
		export kallistoTopName=\$(cat kallisto.out | sort -k 2 -n | tail -n 1 | awk '{ print \$1 }')
		
		echo generating report.html...
		$projectDir/generateReport.sh $sampleName $projectDir/htmlHeader.html $isPairedEnd $primerBedFile $projectDir
	"""
}



// The below process runs once per folder and generates a concise summary of all samples after all other
// executions are over.
process summaryPage {
	input:
		file 'report' from reportCh.collect()
	
	output:
		file "analysisResults" into results_with_summary

	conda "$projectDir/conda/env-python"
	
	shell:
	"""
		$projectDir/generateSummary.sh $projectDir/htmlHeader.html $params.variantDBfile $projectDir
	"""
}



// OPTIONAL: convert html reports into pdf and then generate a combined pdf of all results.
// Useful if need to share the result with external collaborators
process html2pdf {
	input:
		file "analysisResults" from results_with_summary
	
	output:
		file "analysisResults" into analysisResults

	conda "$projectDir/conda/env-gs-wkhtmltopdf"
	label 'high_cpu'
	publishDir "$params.out", mode: 'copy', overwrite: true	
	
	shell:
	if (params.make_pdfs) {
	"""
		echo Generating report.pdf...
		cd analysisResults
		for sampleName in \$(ls */ -d | tr -d '/'); do
			awk '1; /Detected mutations/{exit}' \$sampleName/\${sampleName}_report/report.html > \$sampleName/\${sampleName}_report/temp.html
			echo "Excluded from this pdf version due to file size limitations." >> \$sampleName/\${sampleName}_report/temp.html
			echo "<br>"\$'\n'"</body>"\$'\n'"</html>" >> \$sampleName/\${sampleName}_report/temp.html
			
			wkhtmltopdf --enable-local-file-access --page-size Letter --margin-top 10mm --margin-bottom 0 --margin-left 0 \
				--margin-right 0 --print-media-type --title "Wastewater report" \$sampleName/\${sampleName}_report/temp.html \
				\$sampleName/\${sampleName}_report/report.pdf &
		done
		
		echo Generating summary.pdf...
		wkhtmltopdf --enable-local-file-access --page-size Letter --margin-top 10mm --margin-bottom 0 \
			--margin-left 0 --margin-right 0 --print-media-type --title "Wastewater report" summary.html summary.pdf &
		
		echo Waiting for the conversion processes to complete
		wait
		rm ./*/*/temp.html
		
		echo Merging PDFs...
		gs -dNOPAUSE -dQUIET -dBATCH -sDEVICE=pdfwrite -dPreserveAnnots=false -sOUTPUTFILE=./consolidated.pdf ./summary.pdf ./*/*report/report.pdf
	"""
	}
	else {
	"""
		echo make_pdfs was set to false, so skipped the pdf generation.
	"""
	}
}


analysisResults.subscribe onComplete: {
	println('Pipeline execution complete. Thank you for choosing C-WAP')
}