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.
•
written
by
chen
2 answers
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.
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.
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
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.
Log in to answer this question.
-
How to run SPAdes assembler with bam files?
written by marongiu.luigi
-
Unmapped reads to fastq (hisat2 / samtools / bedtools) not matching pair-end reads
written by juan.crescente •
-
Tutorial: use the tool slicer to split big files (like FASTQ) to smaller files for parallel process…
written by chen
-
Removing shorter reads in PE mapping
written by biolab
-
How to align non overlapping fastq files
written by Rongen •
-
Strand specific pair end sequencing
written by Bioinfonext
-
How to get percentage of reads mapped to each Reference Genome
written by bandanaschapagain •
-
How to obtain insert sequence from paired-end mapping result (SAM/BAM)
written by gundalav
-
Tool: AfterQC: Automatic Filtering, Trimming, Error Removing and Quality Control for fastq data
written by chen
-
paired end illumina reads
written by midox
