Flipped ASE because of mapping issues

[weishwu]

Hi Ernest,

@eturro Sorry to repost this question. You replied to my other post and closed it. I replied there with more questions but the post was gone. I'm not sure if you are still able to see that post, so I'm reposting here:

I asked in that post:
As the Bowtie manual says, "--best --strata" does not apply for paired-end reads. What's the best solution to this? I really need to get only the best stratum.

You replied "This is already done internally by bam2hits, don't worry. Using --best --strata just helps get rid of some of the sub-optimal alignments beforehand."

And below is my follow-up questions:

But I'm seeing some obviously wrong results which was caused by this issue. For example, I have a transcript that contains only a single SNP. When I quantified allele-specific expression at this SNP, it clearly showed strong paternal-biased expression (more than 95% reads support paternal allele), however the results from MMSeq showed strong maternal-biased expression. After some detective work, I figured it out that the major (not 100%) problem is that bowtie could not distinguish between the paternal version and the maternal version of this transcript because it is only a single nucleotide mismatch (because the default allows 2 mismatches) so it reported alignment to both in equal amount. I tried redoing the alignment by requiring zero mismatch then the results showed correct paternal-biased expression. For other transcripts, this does not necessarily cause flipped direction but can reduce the amplitude of allele-specific expression. Only after I required zero mismatch the results started to make sense. But rather than requiring zero mismatch, the ideal way is to let bowtie report only the best stratum.

Let me show you an example:

Below is mmseq results for the transcripts of IGF2. IGF2 is a maternally-imprinted gene in my samples. Column 4 and 5 is log_mu of paternal and maternal versions respectively and column 5 is column4-column5 (log ratio):
ENST00000381392.5 ENSG00000167244.21 IGF2 1.56433 -2.22791 3.79224
ENST00000381395.5 ENSG00000167244.21 IGF2 9.766810000000001 4.06606 5.700750000000001
ENST00000381406.8 ENSG00000167244.21 IGF2 8.01555 -6.4431 14.458649999999999
ENST00000416167.7 ENSG00000167244.21 IGF2 7.79913 -4.31979 12.11892
ENST00000434045.6 ENSG00000167244.21 IGF2 -7.343839999999999 4.61276 -11.956599999999998

As you can see, while all other transcripts are maternally imprinted, ENST00000434045.6 is paternally imprinted. And the first two transcripts show lower degree of imprinting than the next two. When I looked at SNP level ASE, all the SNPs in these transcripts have strong maternal imprinting.

I tried rerunning it by requiring zero mismatch in bowtie mapping, and got these results:
ENST00000381392.5 ENSG00000167244.21 IGF2 7.73735 -4.58045 12.3178
ENST00000381395.5 ENSG00000167244.21 IGF2 9.3704 2.71391 6.65649
ENST00000381406.8 ENSG00000167244.21 IGF2 7.44733 -5.74605 13.193380000000001
ENST00000416167.7 ENSG00000167244.21 IGF2 7.612660000000001 -2.5921 10.20476
ENST00000434045.6 ENSG00000167244.21 IGF2 4.867319999999999 -9.99843 14.86575

As you can see, now all transcripts show strong maternal imprinting, consistent with the SNP level results. Well, I'd expect ENST00000381395.5 to be similar with others, but at least there is no flipping any more.

The diploid transcriptome was created using the phased SNPs, so the transcript level ASE should be consistent with the SNP level ASE when all SNPs in this gene are strongly maternally imprinted.