Any tool to output unmapped or mapped_with_clips reads from a BAM to two fastq files for PE data?

Question

[chen]

Hi,

I am processing pair-end sequencing data. I want to output the reads that are unmapped to reference genome or mapped with soft/hard clips, to two fastq files (say R1.fq, R2.fq).

Is there any existing tool available to do this? If not, I will have to develop it by myself.

Answer 1

[Devon Ryan]

Unmapped reads can be found with samtools view -f 4. If you need those with soft/hard clipping then use pysam and parse the cigartuples item.

[chen]

But this only output the single reads, not the pair.

For pair-end data, if any read of a pair is unmapped or mapped with clips, I'd like to output both reads of this pair.

[ATpoint]

Use -f 5 then. It will give you unmapped but paired reads. You can play around with the flags, using this tool.

[Devon Ryan]

You'll need to do two passes, where you make a hash of the read names in the first pass and then output them in the second. You can do that in python with pysam, C with htslib, or java with htsjdk. I would suggest using python, it'll be quicker to write.

Answer 2

[Pierre Lindenbaum]

select unmappeds read or clipped reads with , say samjdk http://lindenb.github.io/jvarkit/SamJdk.html (many other solutions are available, search biostars.org)

java -jar dist/samjdk.jar -o out.bam -e 'return record.getReadUnmappedFlag() || record.getCigar().getCigarElements().stream().anyMatch(C->C.getOperator().isClipping());'  in.bam

then use picard SamToFastq to convert back to fastq

http://broadinstitute.github.io/picard/command-line-overview.html#SamToFastq

[Pierre Lindenbaum]

I'd like to output both reads of this pair.

more complicated: You need to get the names of the reads after filtering (with samtools view + cut) , and then use http://broadinstitute.github.io/picard/command-line-overview.html#FilterSamReads + READ_LIST_FILE to get the read with the same names.