filtering for uniquely mappable regions

[dphansti]

I have paired end 101 bp reads for some yeast strains. I am currently using bwa to map reads and the GATK to call and filter SNPs and indels. For various reasons I would like to only deal with uniquely mappable reads. I would also like to know which regions of the genome are uniquely mappable.

Are there simple ways to do this?

Can I simply do the following...

samtools -bS -q 1 input.sam > output.bam

when doing my sam to bam conersions? Is there a more appropriate value than 0 for the "-q" argument?

In order to find non-uniquely mappable regions i found every 101 bp sequence in the reference genome and made a fastq of them. I then aligned these and filtered them for MAPQ scores of greater than 0. using the command shown above.

Is this an appropriate method? Is there a different standard method I should be applying?

Thanks so much.

[gringer]

bowtie has a -m <x> option that filters out non-unique reads (i.e. a read that maps to more than x genomic positions). I don't think the -q option of samtools view will work, because you're still likely to get perfect matches that map to multiple positions.

[gringer]

FWIW, when I was mapping the yeast data from the short read archive SRP000637, here are the bowtie commands I ran:

# (#) filter out rRNA sequences using bowtie, seed length of 21
bowtie -p 10 ~/bowtie/scer_rRNA $(printf '%s,' polyAclipped_* | perl -pe 's/,$//') -l 21 -t --un rRNAfiltered_SRP000637_all.fastq 2>bowtie_summary_all_rRNA.txt > /dev/null
# (#) find cDNA matches with <10 hits, direct misses to separate file
bowtie --un non-cDNA_rRNAfiltered_SRP000637_all.fastq --max high-repeat_rRNAfiltered_SRP000637_all.fastq -l 21 -t -a -m 10 --best --strata -S ~/bowtie/Saccharomyces_cerevisiae.EF3.64.cdna.all rRNAfiltered_SRP000637_all.fastq 2>bowtie_summary_cDNA_all.txt > bowtie_cDNA_rRNAfiltered_SRP000637_all.sam

I've used -m 10 there because I wanted to pick up matches to different isoforms, but -m 1 is also possible. The --max option makes sure that these high-map sequences are put in a separate file for further analysis, if necessary.

[cedance]

If you use BWA, the output bam or sam will have an optional flag = XT:A:U set in case of unique reads. So, you could just write a script to pull those lines containing XT:A:U.
Similarly, for tophat there's an optional NH:i:1 (1 for unique hit).

[swbarnes2]

If you use BWA, the output bam or sam will have an optional flag = XT:A:U set in case of unique reads. So, you could just write a script to pull those lines containing XT:A:U.
Similarly, for tophat there's an optional NH:i:1 (1 for unique hit).

True, but you also want the M reads. Those are the ones where the read wasn't unique, but the mate was, so that made the mapping of the repetative bit sure. And I've unfortunately seen R reads where the mate mapped uniquely. I think those are cases where it should have been changed to M, but it wasn't.

So I don't think that pulling from the flags alone will do it. You'll have to pull down pairs, and check to see if either end maps uniquely.

[cedance]

swbarnes2, Thank you. In my case, if both ends don't map uniquely, I throw them. If one of the pair maps uniquely (that would have XT:A:U and its mate with XT:A:M), I consider the unique hit as a single-end read and remove the other one alone.

couple questions, 1) what do you mean here? "so that made the mapping of the repetative bit sure". 2) why do you want the other pair when its not uniquely mapped? (assuming you remove non-unique reads otherwise).

[swbarnes2]

If you have one read that could map at positions 3000, 5000, and 10,000, and the other end best maps at 10,200 uniquely, then the first end must have come from position 10,000, and not from the other two places. So why throw that away?

[cedance]

makes total sense. thank you once again!

[dphansti]

Thanks cedance. Very helpful. The XT:A: tag works perfectly. For others viewing this thread, those tags at the end of the line are all added by BWA so I needed to look here for the information about what they mean. I had been looking at the description of the sam format instead.

But I do have a follow up question. One of the lines has an XT:A:U tag but it also has XA:Z:chr5,+111884617,36M,1;. I realize that this indicates an alternate match. My 'real' alignment also has a CIGAR string of 36M but an edit distance of 0 (NM:i:0). I assume that this means my 'real' alignment is slightly better and therefore I should accept it. But what exactly does edit distance mean? Thanks so much.

[cedance]

X = AGCT
Y = AACA
Edit distance = 2. It denotes the number of operations (depending on the software, it allows for substitution alone or insertions, deletions and substitution) that are required to start at X and end at Y. Here, from G->A and T->A requires 2 changes.

So, in your case, you might have an alternate match (I haven't looked at the optional column for alternate match before..) where the edit distance was not 0 (which means that hit had a mismatch; either ins, del or sub). So, when reads are non-unique, I guess it should mean there are more than 1 read that maps uniquely having at least the same edit distance, I guess. Good to know.