SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Differential Expression: is it better to map reads to genome or transcriptome? dalesan Bioinformatics 26 04-22-2014 04:02 AM
Tophat RNASeq mapping - right reads map and left reads map 50% less airad22 Bioinformatics 2 08-14-2013 08:26 AM
Reads that can align to multiple places ashwatha Bioinformatics 2 09-21-2011 08:21 PM
How to map short reads to a distant genome? ynwh Illumina/Solexa 5 08-03-2011 06:56 AM
How to map 454 reads/contigs to a mitochondrial genome? fruktimport Bioinformatics 2 03-28-2011 11:35 AM

Reply
 
Thread Tools
Old 08-17-2017, 10:11 AM   #1
ammoon
Junior Member
 
Location: pa usa

Join Date: Aug 2017
Posts: 2
Question how to eliminate reads that map to multiple places in genome

I am trying to compare expression of histone genes in an rnaseq dataset using igv; many of the reads map to multiple places in the genome,,not surprising given the homology in this protein family but, i need to get accurate expression data . How do i fix the alignments to get rid of these kinds of reads?
ammoon is offline   Reply With Quote
Old 08-17-2017, 11:14 AM   #2
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

Typically, ambiguously-mapped reads are either given a single random mapping location, or mapped to all locations via secondary alignments. Expression quantification tools understand this.

You can get rid of multimapping reads by filtering out everything with a MAPQ of 3 or less, but that will cause underrepresentation of the expression of homologous genes.
Brian Bushnell is offline   Reply With Quote
Old 08-17-2017, 11:21 AM   #3
ammoon
Junior Member
 
Location: pa usa

Join Date: Aug 2017
Posts: 2
Default

thanks for the response

i tried setting the igv threshold for map quality to 50 and i am still getting reads that map to multiple places based on the "blat read sequence" readout table

some reads start in one histone gene and end in another one 90Kb away...?

how do i tell if my reads were given a single random location or mapped to all via secondary alignments?
thank you
ammoon is offline   Reply With Quote
Old 08-17-2017, 11:48 AM   #4
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 7,080
Default

It depends on what aligner you used and if it followed the consensus of giving a low MAPQ value to the multi-mapped reads. Brian was referring to default behavior of his BBMap aligner. You can either re-map the data (if you are not sure what your original aligner did) or use one of the suggestions in this Biostars thread to filter your BAM files.
GenoMax is offline   Reply With Quote
Old 08-17-2017, 12:38 PM   #5
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

You can also look at the NH tag:

Code:
NH i Number of reported alignments that contains the query in the current record
That's an optional field, but it should tell you how many alignments there are for a read, in the sam file. There's also

Code:
H0 i Number of perfect hits
...but that is much less common.
Brian Bushnell is offline   Reply With Quote
Old 08-22-2017, 01:14 PM   #6
aprice67
Member
 
Location: New York

Join Date: Nov 2012
Posts: 49
Default

Check out figure 12 and the area around it in this paper.

http://journals.plos.org/plosone/art...l.pone.0180904

It might not be such a good idea to remove them.
aprice67 is offline   Reply With Quote
Old 08-23-2017, 04:12 AM   #7
Markiyan
Senior Member
 
Location: Cambridge

Join Date: Sep 2010
Posts: 116
Lightbulb Some commandline examples for mapping quality filtering.

Mapping score distribution depends on mapper program/version used.

First check the actual MAPQ score histogram for your bam file (first 1M of reads):

samtools view $1 | head - -n 1000000 | cut -f5 | sort | uniq -c

Than filter with desired threshold using samtools and -q parameter (assuming the input is test.bam):

#!/bin/sh
MAP_QUAL=30
STRAIN=test
samtools view -q $MAP_QUAL -b ${STRAIN}.bam -o ${STRAIN}.Q$MAP_QUAL.bam -U ${STRAIN}.Q${MAP_QUAL}U.bam


The passing MAPQ filter reads will end up in test.Q30.bam

And failing ones would end up in test.Q30U.bam file.
Markiyan is offline   Reply With Quote
Old 08-23-2017, 06:43 AM   #8
aprice67
Member
 
Location: New York

Join Date: Nov 2012
Posts: 49
Default

Be very careful when dealing with MAPQ scores. It seems every aligner has a different way to determine them and they aren't standardized at all. See here:

http://www.acgt.me/blog/2015/3/17/mo...-documentation
aprice67 is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 07:57 AM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2021, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO