question about identifying difference from sequencing data statistically

[gene_x]

Folks,
I'd like to hear your advice, suggestion and practice on how to statistically identify difference between two sequencing datasets. This applies to many analysis I encountered.

Suppose I have only one control and one experimental sample. I only have one sequencing data on each sample and the sequencing data can be RNA-Seq, ChIP-seq, FAIRE-seq etc, I want to compare the two sequencing data and identify either those most differentially expressed genes for RNA-Seq, or those regions that show the most differentially transcription factor binding between two samples, or those most differentially methylated regions between two samples etc... given I only have one sequencing data on each sample, what would you do to get those differentially "changed" sets of genes/regions?

I know some people calculate rpkm difference for the same gene/region and take the top difference as the most differentially "changed" sets of genes/regions. But this doesn't have statistical significance I guess. So how do you go about doing this in a statistically sound way?

I hope I stated my questions clearly. Thanks!

[dpryan]

I know some people calculate rpkm difference for the same gene/region and take the top difference as the most differentially "changed" sets of genes/regions. But this doesn't have statistical significance I guess. So how do you go about doing this in a statistically sound way?

You don't go about doing things in a statistically sound way, because you can't. The best you can do is make a ranked list of candidates for further validation. Maybe you'll get lucky and someone else will have published a similar enough dataset that you can take a stab at estimating biological variability, but this is generally a terrible way to do things.

[gene_x]

You don't go about doing things in a statistically sound way, because you can't. The best you can do is make a ranked list of candidates for further validation. Maybe you'll get lucky and someone else will have published a similar enough dataset that you can take a stab at estimating biological variability, but this is generally a terrible way to do things.

Can you explain in more detail about how you do a ranked list of candidates? Do you do anything to "control" or "normalize" things?

[dpryan]

Can you explain in more detail about how you do a ranked list of candidates? Do you do anything to "control" or "normalize" things?

That would depend a bit on the experiment type. I'll speak to RNAseq, since that's most of what I do. There, one would generally do the library size normalization in DESeq or edgeR or whatever you're using. For RNAseq data, you can then calculate dispersion with method="blind", followed by the usual tests to get a list of candidates (they'll have p-values, but those are only really useful for sorting the list). Alternatively, one could instead sort by fold-change, after applying a signal strength (i.e. read number) threshold so that the resulting list isn't dominated by genes whose large fold-change is due to having low read counts. In general, it's better to just run more samples.

In general, a ranked list would just be done by fold-change between the two groups. You would obviously need to normalize by library size since otherwise that alone could throw things off (I wouldn't just scale things blindly, often people will use some variant of quantile normalization, google around for examples). RNAseq complicates things a bit since there's a relationship between expression level and dispersion. For that reason, people (who made the poor decision to not use biological replicates) would normally use the more complicated method outlined above.

Anyway, you're really better off sequencing some replicates of things before spending too much time with the analysis.

[gene_x]

For RNAseq data, you can then calculate dispersion with method="blind", followed by the usual tests to get a list of candidates (they'll have p-values, but those are only really useful for sorting the list).

How are these p-values for a list of candidates calculated?

Alternatively, one could instead sort by fold-change, after applying a signal strength (i.e. read number) threshold so that the resulting list isn't dominated by genes whose large fold-change is due to having low read counts. In general, it's better to just run more samples.

Besides fold-change, do you use difference between signals (normalized signal, like rpkm?) instead?

In general, a ranked list would just be done by fold-change between the two groups. You would obviously need to normalize by library size since otherwise that alone could throw things off (I wouldn't just scale things blindly, often people will use some variant of quantile normalization, google around for examples). RNAseq complicates things a bit since there's a relationship between expression level and dispersion. For that reason, people (who made the poor decision to not use biological replicates) would normally use the more complicated method outlined above.

What is the dispersion you mentioned here? Dispersion of reads?

Anyway, you're really better off sequencing some replicates of things before spending too much time with the analysis.

I guess when you have replicates for each group, you'll get statistics on each candidate comparison? sorry for my limited understanding of statistics...

[Jeremy]

n=1 is not a sufficient sample size for statistics or biology.

[swbarnes2]

Right; one control and one experimental can't tell you anything statistically significant. You have no idea how much your control values will naturally vary, so you can't know if the experimental value is significantly different.

One compromise that we use in our lab is to do multiple bio replicates, and then they are pooled, so that only one RNA-seq library prep is done. An average value is better than a single replicate. And any hit is checked with another technology.

Really, you need to tell the people doing this experiment that they are pretty much wasting your time, and the resources of whoever sequenced their samples. They need to do biological replicates, to get a handle on the natural variance for all those genes.

[chadn737]

Besides the obvious problem of no biological replicates, you cannot simply make generalizations about all these different types of data. Analyzing RNA-seq data is very different than DNA Methylation data.

[Jeremy]

Besides the obvious problem of no biological replicates, you cannot simply make generalizations about all these different types of data. Analyzing RNA-seq data is very different than DNA Methylation data.

It's not that different for the topic we are talking about, methylated sites have natural variation in the population too, without biological replicates you have no way to know if the difference you see is expected based on the degree of natural variation or not.
It is simple statistics, I remember hearing it in first year statistics and in first year biology, n=1 does not equal a pattern.

[chadn737]

It's not that different for the topic we are talking about, methylated sites have natural variation in the population too, without biological replicates you have no way to know if the difference you see is expected based on the degree of natural variation or not.
It is simple statistics, I remember hearing it in first year statistics and in first year biology, n=1 does not equal a pattern.

I never said you didn't need biological replicates. I said " Besides the obvious problem of no biological replicates" which basically means "in addition too this problem".

I am saying that the methods of analyzing and detecting methylation differential methylation are different from those of detecting differential expression. How one normalizes the data, the statistical models used....these all differ depending on what sort of analysis one is doing.

[Heisman]

Not to hijack this thread, but I'll be doing some type of statistical analysis like this over the summer and am currently ignorant in regards to how to do this in R. I basically need to teach myself by the summer. Could somebody provide links to tutorials or classical papers that describe how to do this type of analysis (assuming there are biological replicates)? I think that is somewhat what the OP is getting at anyway (without realizing he needed replicates). Any help would be very greatly appreciated. I'll be working with methylation data.

[chadn737]

Look up the vignettes for DESeq and EdgeR. These are two of the best written vignettes I have ever seen. They not only guide the user through the use of their methods, but also do an excellent job of explaining the underlying rationale for the statistics used. You will learn a lot if just starting with these two if nothing else.

[gene_x]

Look up the vignettes for DESeq and EdgeR. These are two of the best written vignettes I have ever seen. They not only guide the user through the use of their methods, but also do an excellent job of explaining the underlying rationale for the statistics used. You will learn a lot if just starting with these two if nothing else.

I typically find it hard to understand those vignettes but thanks for pointing to those two.

[Heisman]

Look up the vignettes for DESeq and EdgeR. These are two of the best written vignettes I have ever seen. They not only guide the user through the use of their methods, but also do an excellent job of explaining the underlying rationale for the statistics used. You will learn a lot if just starting with these two if nothing else.

Excellent, I will start there. Thank you!

[danwiththeplan]

If you already have the sequencing done you obviously can't do this, but adding spike-in controls would go a long way to letting you compare expression levels across datasets. Also, biological reps.