Fastq from unmapped reads

[maxibor]

To answer #187
Comes with a python script (using pysam) to recreate fastq files from the unmapped reads.

However this is using Python3 and the Python in the conda env is currently 2.7 :(

[apeltzer]

But the current code already does this?

eager/main.nf

Lines 846 to 898 in bc55df3

	process samtools_filter {
	tag "$prefix"
	publishDir "${params.outdir}/samtools/filter", mode: 'copy',
	saveAs: {filename ->
	if (filename.indexOf(".fq.gz") > 0) "unmapped/$filename"
	else if (filename.indexOf(".unmapped.bam") > 0) "unmapped/$filename"
	else if (filename.indexOf(".filtered.bam")) filename
	else null
	}
	input:
	file bam from ch_mapped_reads_filter.mix(ch_mapped_reads_filter_cm,ch_bwamem_mapped_reads_filter)
	output:
	file "*filtered.bam" into ch_bam_filtered_qualimap, ch_bam_filtered_dedup, ch_bam_filtered_markdup, ch_bam_filtered_pmdtools, ch_bam_filtered_angsd, ch_bam_filtered_gatk
	file "*.fastq.gz" optional true
	file "*.unmapped.bam" optional true
	file "*.{bai,csi}"
	script:
	prefix="$bam" - ~/(\.bam)?/
	size = "${params.large_ref}" ? '-c' : ''
	if("${params.bam_discard_unmapped}" && "${params.bam_unmapped_type}" == "discard"){
	"""
	samtools view -h -b $bam -@ ${task.cpus} -F4 -q ${params.bam_mapping_quality_threshold} -o ${prefix}.filtered.bam
	samtools index "${size}" ${prefix}.filtered.bam
	"""
	} else if("${params.bam_discard_unmapped}" && "${params.bam_unmapped_type}" == "bam"){
	"""
	samtools view -h $bam | tee >(samtools view - -@ ${task.cpus} -f4 -q ${params.bam_mapping_quality_threshold} -o ${prefix}.unmapped.bam) >(samtools view - -@ ${task.cpus} -F4 -q ${params.bam_mapping_quality_threshold} -o ${prefix}.filtered.bam)
	samtools index "${size}" ${prefix}.filtered.bam
	"""
	} else if("${params.bam_discard_unmapped}" && "${params.bam_unmapped_type}" == "fastq"){
	"""
	samtools view -h $bam | tee >(samtools view - -@ ${task.cpus} -f4 -q ${params.bam_mapping_quality_threshold} -o ${prefix}.unmapped.bam) >(samtools view - -@ ${task.cpus} -F4 -q ${params.bam_mapping_quality_threshold} -o ${prefix}.filtered.bam)
	samtools index "${size}" ${prefix}.filtered.bam
	samtools fastq -tn ${prefix}.unmapped.bam | pigz -p ${task.cpus} > ${prefix}.unmapped.fastq.gz
	rm ${prefix}.unmapped.bam
	"""
	} else if("${params.bam_discard_unmapped}" && "${params.bam_unmapped_type}" == "both"){
	"""
	samtools view -h $bam | tee >(samtools view - -@ ${task.cpus} -f4 -q ${params.bam_mapping_quality_threshold} -o ${prefix}.unmapped.bam) >(samtools view - -@ ${task.cpus} -F4 -q ${params.bam_mapping_quality_threshold} -o ${prefix}.filtered.bam)
	samtools index "${size}" ${prefix}.filtered.bam
	samtools fastq -tn ${prefix}.unmapped.bam | pigz -p ${task.cpus} > ${prefix}.unmapped.fastq.gz
	"""
	} else { //Only apply quality filtering, default
	"""
	samtools view -h -b $bam -@ ${task.cpus} -q ${params.bam_mapping_quality_threshold} -o ${prefix}.filtered.bam
	samtools index "${size}" ${prefix}.filtered.bam
	"""
	}
	}

It even extracts the unmapped data to either BAM, FastQ depending on the users choice. I think we need something else :-(

[jfy133]

@maxibor are these fastq reads pre Trimming and merging, but without human reads?

[maxibor]

@maxibor are these fastq reads pre Trimming and merging, but without human reads?

in this PR, yes

[maxibor]

It even extracts the unmapped data to either BAM, FastQ depending on the users choice. I think we need something else :-(

But this step works on post AR fastq files

[jfy133]

Sorry - that close was courtesy of Maia

[jfy133]

I think you are good to go for testing once python version fixed.

Minor thing: the help message/description is unspecific of what actually is being output. Thus maybe Alex' confusion

[maxibor]

Updated Conda env so Pysam should now come. But the test should fail because the Docker container isn't rebuilt as of now with the new Pysam dependancy

[jfy133]

(Official review)

Looking more in depth (please correct me if I misunderstand anything).

Major

I can't comment on the python section. However,

1. The BAM file as input into this tool should be immediately after BWA (sorted.bam) as it should include all possible mapped reads, not just ones that mapped exactly. i.e. if you use the reads post-samtools filter -q 37, the discarded mapped reads would not be filtered by the new module, because they are removed from the mapped BAM file after the samtools filter.

2. Is L921 flipped? Shouldn't it be if (params.singleEnd) { - as you only indicate a single output file in that conditional block?

Minor

For the process name (L906), flag itself (L363, L218) and help message (L108), I suggest the following extra precision instead of -unmap:

--strip_input_fastq                       Create pre-Adapter Removal FASTQ files without reads that mapped to reference (e.g. for public upload of privacy sensitive non-host data)

publishDir should renamed to something like /samtools/stripped_fastq rather than unmapped_fastq, as the latter file already exists with the --extract_unmapped functionality.

Equally, the output FASTQ names should be e.g. stripped.fwd.fq.gz and stripped.rev.fq.gz or something. I think R1/R2 would be dangerous for novices re-analyzing the data with the same pipeline, if they forget to add underscores (i.e. there are two R1/R2s in the name). This may then lead to funky input regex for a EAGER2 re-run and subsequent errors.

The actual names/flags can be discussed here - you don't have to go exactly with my 'stripped' suggestion.

[maxibor]

Followed the suggestions of @jfy133
Also user can now choose the strip mode: Either stripping completely the mapped reads (--strip_mode strip) or just replace the sequence of the mapped reads by N (--strip_mode replace)

[jfy133]

@ivelsko owes you a beer ;)

[maxibor]

@apeltzer Can you have a look and merge if ok ?
The tests are working with the updated docker image using the new conda env. Travis should be fixed after the docker image is rebuilt.

[apeltzer]

Looking good 👍

[apeltzer]

Just failing because of missing python3.6 and pysam - merging to get the dev image updated 👍