Question: How to filter/view reads mapped to multiple locations using samtools or pysam
gravatar for AlchemistMoz
3.1 years ago by
AlchemistMoz20 wrote:

I've looked through the documentation for samtools and pysam but haven't found what I'm looking for. Pysam can check if read.is_secondary which returns true if the read is not primary alignment. Whilst in samtools, I think the corresponding flag is:

256 0x100 secondary alignment

But what I'm looking for is a way to view/filter reads that are mapped to multiple regions, is there a way to do this with samtools or pysam?

pysam samtools python • 4.0k views
ADD COMMENTlink modified 8 days ago by jaeyoung.chun30 • written 3.1 years ago by AlchemistMoz20
gravatar for Tom_L
3.1 years ago by
Tom_L320 wrote:

According to the SAM format specification, I believe you are looking for the NH tag.

NH: Number of reported alignments that contains the query in the current record

Page 6:

To filter them:

samtools view -h myBAM.bam | grep -P "(NH:i:1|^@)" | samtools view -Sb - > myBAM_filtered.bam

This will basically only grep SAM header (@ lines given by the -h argument) and single mapping reads (NH:i:1). It will also generate a BAM file (-b argument) based on SAM input (-S argument).


ADD COMMENTlink written 3.1 years ago by Tom_L320

Note that this is an optional tag, so it may or may not be present, depending on the aligner and command. I suggest filtering reads based on the required field mapq; a read with mapq 3 or lower is not uniquely mapped, and a read with mapq >3 is uniquely mapped.

ADD REPLYlink written 3.1 years ago by Brian Bushnell17k

Hi Brian, i would like to extract only the reads with mapping score 0, but as far as i understand the parameter "-q" only filters out what is equal or smaller that the argument. But i would actually isolate and study on only mapq 0 reads. Can you help me with that?

ADD REPLYlink written 11 months ago by stefano.meucci0

If you have bamtools,

bamtools filter -in pre-filtered.bam -out post-filtered.bam -mapQuality "0"
ADD REPLYlink modified 8 days ago • written 8 days ago by jaeyoung.chun30

With your command,there may be some multi alignments containing such as NH:i:12 in NH tag.Maybe the pattern "(NH:i:1\b|^@)" is better.

ADD REPLYlink written 24 months ago by dh9084075430
gravatar for jaeyoung.chun
8 days ago by
jaeyoung.chun30 wrote:

If you have bamtools, you can do this in a little more human-friendly way:

bamtools filter -in pre-filtered.bam -out post-filtered.bam -tag "NH:>1"

Above will create a new bam file called post-filtered.bam. If you want to test by outputting the first 300 lines to stdout (instead of writing to a file):

bamtools filter -in pre-filtered.bam -tag "NH:>1" | bamtools convert -format sam | head -n 300
ADD COMMENTlink written 8 days ago by jaeyoung.chun30
Please log in to add an answer.


Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1163 users visited in the last hour