search

eturro / mmseq

Haplotype, isoform and gene level expression analysis using multi-mapping RNA-seq reads
GNU General Public License v2.0
67 stars 20 forks source link

Quantified expression does not go along with read counts #45

Open weishwu opened 4 years ago

weishwu commented 4 years ago

Hi,

I've run bowtie, bam2hits and MMSeq to quantify allele specific expression, and observed some discrepancy between the alignment counts and the final log_mu.

One problematic transcript is ENST00000381392.5 which is an IGF2 transcript. I have ENST00000381392.5_A and ENST00000381392.5_B in the reference representing paternal and maternal transcripts respectively.

In the Bowtie alignment, I counted the alignment counts in them and these are the results: ENST00000381392.5_A: 725,664 ENST00000381392.5_B: 721,570

ENST00000381392.5_A has more counts than ENST00000381392.5_B, so the expression of the former should be higher than the latter. However the log_mu from MMSeq is: ENST00000381392.5_A: -4.56312 ENST00000381392.5_B: 4.48249

How could this happen?

Thanks.

eturro commented 4 years ago

Presumably because reads that align to _A also align to some other transcript but this doesn't happen with reads that align to _B?

You can look at this using hitstools inspect

On 13 Dec 2019, at 16:33, weishwu notifications@github.com wrote:

Hi,

I've run bowtie, bam2hits and MMSeq to quantify allele specific expression, and observed some discrepancy between the alignment counts and the final log_mu.

One problematic transcript is ENST00000381392.5 which is a IGF2 transcript. I have ENST00000381392.5_A and ENST00000381392.5_B in the reference representing paternal and maternal transcripts respectively.

In the Bowtie alignment, I counted the alignment counts in them and these are the results: ENST00000381392.5_A: 725,664 ENST00000381392.5_B: 721,570

ENST00000381392.5_A has more counts than ENST00000381392.5_B, so the expression of the former should be higher than the latter. However the log_mu from MMSeq is: ENST00000381392.5_A: -4.56312 ENST00000381392.5_B: 4.48249

How could this happen?

Thanks.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub, or unsubscribe.

weishwu commented 4 years ago

Yes, it's true. Because IGF2 is a maternally-imprinted gene in our sample, so all the paternal IGF2 transcripts have much more alignments than their maternal counterparts. And because multi-mapping is allowed in the alignment, many reads can be aligned to multiple IGF2 paternal transcripts.

Most of the other IGF2 transcripts have the log2_mu in the right direction, that is, much higher in _A than in _B, however ENST00000381392.5 is just opposite.

Is this an expected behavior of MMSeq for imprinted genes that have a bunch of isoforms? It seems to distort the true biology (i.e. all IGF2 transcripts should be maternally imprinted). Is there a way to correct it?

Since our reads were from a targeted capture of a list of imprinted genes in the studied tissue, we have a lot of cases that are like IGF2.