This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Extract entire reads, rather than aligned segments, from a SAM/BAM file

Hello,

I am trying to extract long reads from the Blasr .sam output. However, I want to collect the entire raw reads that produced alignments, and not just the aligned segments.

Is there an easy way to do this?

Here is the blasr command I am using:

blasr <reads.fasta> <reference.fasta> -noSplitSubreads -clipping none -sam -out output.sam
sam bam blasr
• 2,609 views • link
written by gbdias gbdias

Grab the read ID's from your SAM file and then pull them out of the original fastq by using filterbyname.sh from BBMap.

• 0 views • link
written by GenoMax GenoMax

1 answer

I'd assume that the .bam does not have that information. You'd have to get the read ids and get them from the fastq. Something like the following could work. That's what I used to do to extract paired-end reads (in that case the bam needs to be sorted by name).

grep -A3 -F -f <(samtools view <bamfile> | cut -f 1) <source-fastq> | grep -v "\-\-" > <extracted-reads-fastq>
• 0 views • link
written by cschu181 cschu181

Thank you. This strategy worked.

• 0 views • link
written by gbdias gbdias

An alternative I found is to use seqtk. It takes a list of sequence names in a file and extract the sequences from FASTA/Q.

seqtk subseq in.fq name.lst > out.fq
• 0 views • link
written by gbdias gbdias

Log in to answer this question.

More posts like this