SEQanswers

SEQanswers (http://seqanswers.com/forums/index.php)
-   Bioinformatics (http://seqanswers.com/forums/forumdisplay.php?f=18)
-   -   RNA-seq and normalization numbers (http://seqanswers.com/forums/showthread.php?t=586)

zee 09-22-2008 01:15 AM

RNA-seq and normalization numbers
 
In many of the experiments our lab is doing with Illumina reads, we always seem to end up with the task of normalizing data.
If I have 3 experimental conditions, I've sequenced a lane for each and there needs to be a way to compare counts of my known RNA in mirbase to my sequence reads mapped to the genome (with MAQ / novoalign).

I've read about people doing counts as reads per million and log transforming these values to fit Poisson distribution, but it's sprung multiple ideas in my mind. Would this be as simple as dividing my counts for each experiment by

1) 1 Million
2) the total number of reads sequenced
3) the total number of uniquely mapped reads


I'm inclined to option (3) because that represents the amount of usable sequence data.

I'm just wondering if anybody has a more intelligent way of tackling this problem with nextgen data or perhaps there's some software to help out.

I have :

a) Alignment locations of all reads on a ref. genome for each experiment
b) Location of my reference RNAs on the same genome

I am already able to count the number of overlapping locations with each reference RNA in each experiment, and that gives me raw counts.

I have about 4 experiments, but this varies from study to study.

Chema 09-22-2008 03:11 AM

What we are doing is to apply kind of quantile normalization. After to obtain a score for each miRNA , since we are expecting that most of the miRNA have similar expression in both tissues at study, we normalize the score values to have similar distribution.

The main reason to apply this kind of normalization is that we were very worry about the different coverage of the different samples at study.

ECO 09-22-2008 07:46 AM

http://www.ncbi.nlm.nih.gov/pubmed/18516045

This paper introduces a concept RPKM for quantifying tags in RNA-seq. Might be worth a look if you haven't seen it already.

Chema 09-23-2008 05:26 AM

Hi,

The point with RPKM that I do not like, it is that I do not feel that it can handle different coverages. Perhaps I can explain it better through an example.

Let say that we are working with a genome with three genes A, B and C with the same length (I know not very realist, but just an example), and we want to study their expression in two conditions 1 and 2.

The real expression of the genes is:
Condition 1 Condition 2
Gene A 1 1
Gene B 1 1
Gene C 1 0

We run a RNA-seq experiment and we get the next number of reads
Condition 1 Condition 2
Gene A 3*10^5 4.5*10^5
Gene B 3*10^5 4.5*10^5
Gene C 3*10^5 0
Total 9*10^5 9*10^5

Translate to RPKM , since they have the same length, it should be something like:
Condition 1 Condition 2
Gene A 333333 5*10^5
Gene B 333333 5*10^5
Gene C 333333 0

As you can see, it seems that Gene A and B are also differentially express. This is because, since the expression of gene C is lower in condition 2 than 1, we have more reads that will improve of the coverage of the other genes.

Anyway, I think that always it is nice to normalize the data in some way. Mainly, when you are working with so low number of replicates.

vruotti 10-07-2008 09:46 AM

RPKM concern?
 
Hi Chema,

We have used Wold's (RPMK) method and got very good results. Are you saying that if a particular gene has less coverage (coverage=fewer number of mapped reads) this gene will contribute to a change regarding the accuracy of detecting the expression of the other genes? I think the this is the whole point of normalizing. The difference is that the overall expression should not be too different. Also, can you please check your numbers for us? Are they really 333,333 or 333,000,000 for genes a, b and c in condition 1 after converting to RPKM?

Maybe I'm doing this wrong. The formula I see in their paper is:

RPKM = 10^9 x C / NL, which is really just simply C/N

C= the number of mappable reads that felt onto the gene's exons
N= total number of mappable reads in the experiment
L= the sum of the exons in base pairs.

So, let's plug in your numbers.
Condition 1 Condition 2
Gene A 3*10^5 4.5*10^5
Gene B 3*10^5 4.5*10^5
Gene C 3*10^5 0
Total 9*10^5 9*10^5

Translate to RPKM , since they have the same length, it should be something like:
Condition 1 Condition 2
Gene A 333,000,000 500,000,000
Gene B 333,000,000 500,000,000
Gene C 333,000,000 500,000,000

If you look at these numbers you could argue whether the two expression values are differentially expressed. They are not that far apart. Sorry, did I miss your point? Can you explain your concern again?

Victor

Chema 10-10-2008 07:45 AM

Hi Victor,

What I wanted to say is that if you have a group of genes that are less expressed in condition 2 than 1, then, you will have less reads that represent these genes in condition 2 than 1, and therefore the genes that are not in this group will be represented for more reads in condition 2 than in condition 1 (not necessarily, because they are differentially expressed). But I guess this problem will be only important when you have a high number of genes differentially expressed.

Yes, you are right the numbers are as you said. Only:

Translate to RPKM , since they have the same length, it should be something like:
Condition 1 Condition 2
Gene A 333,000,000 500,000,000
Gene B 333,000,000 500,000,000
Gene C 333,000,000 0

About if the difference in gene A and B are significant. Well, I didn’t generate any noise for these data, this is a theorical example, they are not estimations if not the real values of the variable. So, if they are different is because the method will declare it as differentially expressed.

Any way, if you want to see a bigger difference, just run an example with more genes that are downregulated in condition 2. You will see that the diference for gene A and B increase.

vruotti 10-11-2008 08:07 AM

normalization
 
Hi Chema,
Yeah I see what you mean. I guess most of the conditions we ran so far do not have the a huge amount of genes downregulated in comparison with the rest of the genes studied. We are doing whole genome RNA-seq studies and I think the house keeping genes do help avoid this problem.

In case having a huge number of genes that are indeed different, then maybe we have to come up with a slight modification of the RPKM formula/method to cope with this elevated expression or lack of. Maybe include a coverage based on the estimated level of expression variable.

Does anybody else used the RPKM method? Can you share with us the different normalization methods you might use for doing RNA-seq outside of RPKM?

Victor

fabio25 11-12-2008 09:04 AM

Dear everybody,
I would be very interested in asking how you fix the miRNA targets.
just use mirBASE or have you a much nicer way?
Thanks a lot,
Fabio

doxologist 02-27-2009 06:01 AM

Hi Victor:

It seems that from your problem, looking at percentage of the tags a certain gene occupies would solve the problem... and ALSO like suggested, quantile normalization (assuming that the differentially regulated genes do not make up a large percentage of your genes).

steffenp 01-15-2010 03:55 AM

RPKM for Tag-Seq data
 
Hello everyone!
My question is also related to the RPKM normalization. Is this measure suited for Tag-Seq. data, where we have only reads at the 3'end and not all over the whole transcript. My suggestion is no. What would be an adequate normalization for this kind of data?

RockChalkJayhawk 04-02-2010 07:42 PM

Quote:

The point with RPKM that I do not like, it is that I do not feel that it can handle different coverages.
I totally agree with you Chema. RPKM is biased for testing differential expression for longer genes. See the following paper:

HTML Code:

http://www.biology-direct.com/content/4/1/14
and also other papers by Oshlack and Wakefield.

Simon Anders 04-03-2010 01:55 AM

If you want to test for differential expression, it is a good idea to stay on the level of raw, integer counts, and not use RPKM or related data that is normalized by transcript length. This is because significance depends on the number of actual reads that you count. If you have low count you need to see a high fold-change to call significance.

See this thread for more details: http://seqanswers.com/forums/showthread.php?t=4349 (especially from post #6 onwards)

If you work with count data, your testing procedure needs to be aware of the ratios of sequencing depths of the libraries. This functionality is offered by several tools, namely edgeR, DESeq, and cuffdiff. I recommend DESeq, of course, as it is our work. ;-)

Note that this does not alleviate the bias towards longer genes: If two genes have the same expressions (same number of transcript molecules per volume) in two samples and hence the same fold change, the longer one may be called significant and the shorter one not, because the longer one produces more fragments.

If this bothers you, you have a couple of options:
- use Tag-Seq instead of RNA-Seq
- additionally filter with with a rather large threshold on log fold change
- for GO enrichment test and the like, use the test by Young et al. (2010), which takes the length bias into account: http://genomebiology.com/2010/11/2/R14

Here is a figure, that shows, how the log fold change required for significance (red dots: genes with significant DE; black dots: other genes) depends on the counts when using DESeq for testing:

http://www.embl.de/~anders/misc_data/DE_fly_sm.png

For more information, see our paper: http://dx.doi.org/10.1038/npre.2010.4282.1

Simon Anders 04-03-2010 02:05 AM

To estimate the library size, simply taking the total number of (mapped or unmapped) reads is, in our experience, not a good idea.

Sometimes, a few very strongly expressed genes are differentially expressed, and as they make up a good part of the total counts, they skew this number. After you divide by total counts, these few strongly expressed genes become equal, and the whole rest looks differentially expressed.

The following simple alternative works much better:

- Construct a "reference sample" by taking, for each gene, the geometric mean of the counts in all samples.

- To get the sequencing depth of a sample relative to the reference, calculate for each gene the quotient of the counts in your sample divided by the counts of the reference sample. Now you have, for each gene, an estimate of the depth ratio.

- Simply take the median of all the quotients to get the relative depth of the library.

This is what the 'estimateSizeFactors' function of our DESeq package does.

Siva 04-04-2010 09:52 PM

Re: Bowtie-Tophat SAM output for read count assembly
 
Quote:

Originally Posted by Simon Anders (Post 16467)
If you want to test for differential expression, it is a good idea to stay on the level of raw, integer counts, and not use RPKM or related data that is normalized by transcript length. This is because significance depends on the number of actual reads that you count. If you have low count you need to see a high fold-change to call significance.

See this thread for more details: http://seqanswers.com/forums/showthread.php?t=4349 (especially from post #6 onwards)

If you work with count data, your testing procedure needs to be aware of the ratios of sequencing depths of the libraries. This functionality is offered by several tools, namely edgeR, DESeq, and cuffdiff. I recommend DESeq, of course, as it is our work. ;-)

Hi Simon
I have used Bowtie, Tophat and Cufflinks to align and assemble maize RNA-seq data. Cufflinks reports FPKM values which our statistics collaborator is not able to use as it has already been normalized. Can I use the Bowtie, Tophat generated 'SAM' output in some other program to assemble the data in way that I will have absolute read counts per gene or transcript? If so what other programs would you recommend?

thanks
Siva

lpachter 04-05-2010 08:50 AM

Dear Siva,

I'd like to understand why your statistics collaborator cannot use the Cufflinks FPKM values together with their confidence intervals?

Regarding absolute counts, the whole point of Cufflinks is that it is not possible to obtain absolute read counts per transcript, because for many reads there is ambiguity as to which transcript they belong to. Cufflinks is probabilistically assigning reads to transcripts and thereby able to estimate expression of individual transcripts.

RockChalkJayhawk 04-05-2010 09:23 AM

Quote:

Originally Posted by lpachter (Post 16518)
Regarding absolute counts, the whole point of Cufflinks is that it is not possible to obtain absolute read counts per transcript, because for many reads there is ambiguity as to which transcript they belong to. Cufflinks is probabilistically assigning reads to transcripts and thereby able to estimate expression of individual transcripts.

Why not just take the FPKM values categorized by cufflinks for each individual transcripts and multiply them by transcript length and number of aligned sequences? That would give you the tag count per transcript.

lpachter 04-05-2010 09:58 AM

Thats correct- the procedure RCJ suggests will give you an estimate of the actual tag count for each transcript.

Siva 04-05-2010 08:32 PM

Quote:

Originally Posted by lpachter (Post 16518)
Dear Siva,

I'd like to understand why your statistics collaborator cannot use the Cufflinks FPKM values together with their confidence intervals?

Regarding absolute counts, the whole point of Cufflinks is that it is not possible to obtain absolute read counts per transcript, because for many reads there is ambiguity as to which transcript they belong to. Cufflinks is probabilistically assigning reads to transcripts and thereby able to estimate expression of individual transcripts.

Dear L Pachter
Thank you very much for your reply thanks also to RCJ. I understand your point about probabilistically assigning reads to transcripts. I will get back to our statistical group about this. However the procedure suggested by RCJ and seconded by you has me a bit confused. I used to think that the algorithm that Cufflinks uses assigns reads based on probability and it will be different for each transcript and also vary according to genome location. So how can a mere multiplication of FPKM by transcript length and aligned sequences give us tag count per transcript? Am I missing something too obvious here?

thanks
Siva

RockChalkJayhawk 04-06-2010 05:46 AM

To answer your question, see this page <a href="http://cufflinks.cbcb.umd.edu/howitworks.html"> here</a>
Quote:

To estimate isoform-level abundances, one must assign fragments to individual transcripts, which may be difficult because a read may align to multiple isoforms of the same gene. Cufflinks uses a statistical model of paired-end sequencing experiments to derive a likelihood for the abundances of a set of transcripts given a set of fragments.
Therefore, cufflinks is using probabilities to assign individual fragments to individual transcripts. Since the units of the counts are in fragments per KB transcript per million reads, you can convert them back to raw integers by the multiplying the length of a given transcript and number of reads in millions.

RockChalkJayhawk 04-06-2010 06:08 AM

Unfortunately I just noticed that when cufflinks calls transcripts, it doesn't report thier length. Only in the tmap reference files. All the other files only show the genomic coordinates, which isn't as helpful.

I will e-mail the author to see if this functional;ity can be easily added or worked around.


All times are GMT -8. The time now is 03:17 PM.

Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2020, vBulletin Solutions, Inc.