[chadn737]

Not to open a new tread for my question:

I was wondering what is the optimal number of samples per lane on GAIIx if we want to do an RNA-seq of unsequenced plant?

That is the wrong question. Unless you are doing a massive number of samples, you can and should multiplex all your samples and spread it across however many lanes of you need to achieve the desired depth.

The more important question is how many samples do you need and how many reads do you need per sample. That will be determined by the goal of the experiment, not the sequencer.

[vebaev]

Not to open a new tread for my question:

I was wondering what is the optimal number of samples per lane on GAIIx if we want to do an RNA-seq of unsequenced plant?

[DerSeb]

This would be my estimates, looking at human, mouse or rat:

1. for a run of the mill (not looking for anything fancy, microarray replacement type of RNA-seq) 60M
2. a run looking for alternate splicing 150M
3. looking for an elusive gene 100M
4. nothing is known about this transcriptome ?

About the PCR duplicates: yes, I agree that you will sequence more and more duplicates. But with bigger depth I'm not interested in finding new reads. Let's say I capture one exon with 5-10 reads, these reads are present in my library and it will still give me much better extimates of expression if I sequence several replicates of them. So to do a comparison between samples I still consider the information of how many fragments are present in my library as important.
Therefor I do not discard duplicate reads but I use them for quantification. For DNA sequencing I understand why one would discard duplicate reads, but not for quantitative analyses.

[sdriscoll]

The fact that different programs cuffdiff/bayseq/DEseq/DEXseq give different lists of genes makes for yet another interesting discussion also in another thread.

indeed.

[severin]

start with the question

Hi,

I just got my first RNA-seq dataset alignment results today, and I am wondering what would be considered a "good run".
This is a pilot experiment, whose objective is to determine how many samples we can multiplex in one lane of Illumina HiSeq to get the information we want and the best trade off between cost, specificity and sensitivity.

I have some basic questions I ask everyone that asks me this.

1) What is the biological question you are trying to answer? (depending on your answer to this will determine if you want gene level expression differences or differential splice variation)
2) What is your organism and how different are your biological replicates from each other at a genomic level (clones, siblings, population).
-Organism matters because different species have different evolutions and polyploidy like events(more important for plants than animals).
-large variation (overdispersion) between biological replicates is common and more so depending on how you controlled the environment.
3) Are you hoping to find differential expression of a particular gene or is it ok to find differential expression of gene families (multiple mapping can affect gene expression estimates)
4) How good is your annotation? Do you anticipate new genic information not contained in the current gff file?

This just scratches the surface. I had another thread discussing what people are using in the counting (paired concordant, too-long, scramble, inverse, half-mapping, translocations, unpaired reads that map individual but not together and then normalization for paired vs unpaired data)

The fact that different programs cuffdiff/bayseq/DEseq/DEXseq give different lists of genes makes for yet another interesting discussion also in another thread.

[GenoMax]

So before we go to page 2 in this interesting discussion .. which makes it harder for someone coming along afterwards (or searching) to find the conclusive information

Shall we start collecting read levels people feel comfortable recommending (in "Mil" reads for a genome of "X" size).

Possible scenarios (feel free to add more):

1. for a run of the mill (not looking for anything fancy, microarray replacement type of RNA-seq) ---
2. a run looking for alternate splicing ---
3. looking for an elusive gene ---
4. nothing is known about this transcriptome --

[sdriscoll]

Not sure why this bothers you. I think they are trying to address the question: How few reads am I likely to be able to get away with while retaining most of the information I need. So you can look at the curves they generate and dial in a number that might work for the genes you are interested in.

I don't think their results were intended to surprise.

--
Phillip

I'm not completely sure why it bothers me but it does non the less. I may eventually be able to put it into words. Maybe because within those plots you're drawn into comparing the different subsets of the complete data set to one another. So in a 100M read set you might be looking at what 50% of those reads quantifies as to 60% of the reads. I'd wager that is not as informative as looking at actual technical replicates sequenced to different depths. In the subset approach 83% of the reads in the 60% subset were quantified previously in the 50% subset so very little new information is added into the quantification and the chance that the quantification is different isn't nearly as great as a 20% quantification compared to a 10% quantification where only 50% of the data was in the previous subset. If you think of it in the other direction, starting with 100% of your data, how much of an impact would you expect to see to the quantification by removing 10% of the data. What about 70%?

I guess what I'm getting at is I think those plots show a property of sub-setting a fixed size sample at different percentages. It's not a property of RNA-Seq data and I'm afraid that the analysis shows common results independent of the actual number of reads collected. For example if you've got 40 million reads the analysis might show that once you're past 20 or 25 million reads you've got a pretty good picture of what's happening at 40 million reads. If you had 900 million reads it might show you that around 400 or 500 million reads you've got a good picture of what's happening at 900 million reads. I suspect that in all cases when you have N total reads then this analysis will show you that you have a pretty good picture of the quantification of N reads at N/2 or 2N/3 reads. It's just a property of percentages. Take 10% of your data and that's 100% different than no data. At 20%, 50% of the data came from the last subset. At 30%, 66% of the data was quantified in the 20% subset. By 60%, as I mentioned before, you quantified 83% of the that sample in the 50% subset. At each additional subset the the chance of the quantification looking different from the previous subset gets smaller and smaller independent of the numbers of reads involved in each subset.

Of course depth is important for other things - like guaranteeing good coverage of specific length transcripts down to some RPKM level and for splicing analysis. But to say that the gene expression values at 100 million reads are more correct than the values at 50 million seems to me irrelevant and it's certainly not demonstrated by comparing 50% of your reads to 100% of your reads within a single sample. The gene expressions of a single sample sequenced to 500 million reads are still just from a single sample. The distribution of expressions from 10 biological replicates sequenced at 50 million reads each would be much more reliable...which is obvious, right? It's just like any other kind of experiment.

[pmiguel]

I think the deeper the better. No dataset has yet reached the maximum coverage needed to answer all questions.

Past a certain point you might begin to see limits imposed by the library itself. Your major bottleneck is how many molecules you have of your library prior to enrichment PCR. That is, you converted your RNA (probably after ribo-depleting and fragmenting) to DNA and ligated adapters onto it. At that point how many "good" library molecules do you have?

We actually have started doing qPCR at this step. The results are quite variable. But if you are seeing that your concentration is 10 pM -- that is ~6,000,000 "good" molecules/ul. For 20 ul of library you are looking at 120 million molecules. If you want to sequence 30 million reads -- no problem. But at 10x deeper you are guaranteed to be wasting more than 1/2 your reads on PCR duplicates.

--
Phillip

[DerSeb]

I think the deeper the better. No dataset has yet reached the maximum coverage needed to answer all questions.

If you want to only look at Gene expression, probably ~ 60 Million reads would be ok (it's my feeling).

But in my experience, we often hit the limits in our datasets (we have ~140 mill 100 bp PE reads, several biological replicates which is very important!)

What do you want to do after the global gene expression analysis? As soon as you look into specific candidates you will be happy for every read.

Often the expression of your gene of interest might be low (like a FPKM of 3) and you are interested in an alternative first exon. Then you just might end up having to compare 5 reads versus 1

The problem is that expression levels are so different and small regions like exons have much lower read counts than entire genes.

[pmiguel]

you'll see is the values tend to be closer to the mean of all the numbers as you use larger subsets of the data. so...duh?

Not sure why this bothers you. I think they are trying to address the question: How few reads am I likely to be able to get away with while retaining most of the information I need. So you can look at the curves they generate and dial in a number that might work for the genes you are interested in.

I don't think their results were intended to surprise.

--
Phillip

[Geneus]

that seems to be the state-of-the-art.

LOL...now that is funny. But I thought RNAseq has been evolving as a replacement to the former "state of the art"...expression microarrays. Data I've seen in the past from one very large service provider compared a given number of read to the genes detected/quantified compared to the number of genes detected/quantified using a microarray. Based on their work 20M reads was more than sufficient to outperform the microarray. But maybe that was too simplistic an analysis.

[sdriscoll]

so there's one thing that bothers me with these "depth" investigations when they collected N million reads and then randomly subsetted the reads, did alignments and quantified expressions at each subset. so they find some variation across the subsets that decreases as they approach the full alignment.

wouldn't that be obvious though? i mean - as you're using more and more of the total, amount of the data wouldn't there naturally be less and less variation between adjacet subset sizes just due to the fact that you're including more of the data? so the differences between two random subsets of 10% of the data would naturally have more variation than two random subsets of 75% of the data. there is potential for much more overlap between the two 75% subsets than the to 10% subsets.

i guess what i mean is it doesn't really make sense to compare the data to itself over and over again - that's why there aren't statistical tests designed for such comparisons and of course 75% of the data is going to look more like 100% of the data than 10% of the data. duh...right?

try this in R:

# 100 random values between 0 and 10
sset <- runif(100,max=10)
# compute mean of different random subsets at different percentages of the
# total set of numbers at 10% increments
mmean <- 0
for(i in 1:10) { mmean[i] <- mean(sample(sset(,i*10)) }
plot(mmean,type='b')

you'll see is the values tend to be closer to the mean of all the numbers as you use larger subsets of the data. so...duh?

[billstevens]

Hey, I just wanted to say this was a very helpful discussion. As the only person in our lab who is just starting RNA-Seq, it can seem very intimidating, and in many things, there are set standards. I wanted to point out that sometimes when we first start, we don't realize how complex it all is. At the beginning, you can sometimes think,
"hey, what's a good depth for my genome? I'll figure out what all the things that are interesting later."

[sdriscoll]

"give it a shot and see what you get".

that seems to be the state-of-the-art.

[pmiguel]

fair enough!

so then can you tell me how many reads I should collect if i'm interested in genes expressed below the 5 FPKM level? is it worth it to spread a single sample across 8 lanes of a HiSeq 2000 run? those are the exact questions researchers want answers to. for example when you're proposing a grant for a sequencing based project and you have to think in terms of cost.

I was going to write "I am in exactly the same position here in the Purdue genomics core." But I have never had anyone pose the question you just did. If they know what 5 FPKM means, they probably will already know how many reads they want. But I get the more general "how many reads/sample is enough." To tell you the truth I mainly just tell them how many reads/sample the head of our bioinformatics core recommends. If they don't like that answer I can expand somewhat. But unless you have some idea what the shape of your transcriptome is, there is not much to add other than "give it a shot and see what you get".

--
Phillip