Colleagues,
Imagine I have RNAseq data in the form of reads count measured in 2 environments, with 3 replicates each, in a single sequencing experiment. Like this:
Env reads
A 9
A 10
A 11
B 4
B 5
B 6
Of course I can treat each reads count as an independent variable and do a t test and I will have very low power. Instead, I want to be able to treat each read as an observation and then notice that in a single given sequencing effort I got 15 reads in samples from environment B, but 30 reads in samples from environment A. This is highly unlikely to be observed by chance and my power of detecting differences between environments A and B is way, way higher and will increase with the sequencing effort, unlike that of a t-test. (In a typical RAseq experiment number of reads is much higher, for some genes often in hundreds per sample; 2fold count difference will be very significant).
I want to be able to do this contingency test with a more complex model, perhaps with 3 different independent variables (and replicate as a block effect).
The easy solution is then this: I create an additional categorical variable, lets call it "Model" with values Observed and Expected, with "Observed" corresponding to the lines above and "Expected" in equal number of lines, but with the column :reads containing the expected count if no effect is significant, i.e., in, this case, average read count (assuming equal effort in each replicate).
And then run this:
Fit Model( Freq( :reads ), Y( :Model ), Effects( :Env, :rep ), Personality( Nominal Logistic ), Run( Likelihood Ratio Tests( 1 ), Wald Tests( 0 ) ) )
...or this (which gives me identical results, I don't quite understand why, but I am happy with that):
Fit Model( Freq( :reads ), Y( :Model ), Effects( :Env, :rep ), Personality( Generalized Linear Model ), GLM Distribution( Binomial ), Link Function( Logit ), Overdispersion Tests and Intervals( 0 ), Name( "Firth Bias-adjusted Estimates" )(0), Run )
Am I doing something ridiculously wrong or does this make sense as a valid test?
If it makes sense, is there perhaps a way to achieve the same without the need to create additional lines with expected counts?
Thanks in advance!!
L.Yampolsky
Hi @l_yampolsky,
To be clear, for each environment you have three samples (replicates). But instead of combining the replicates for each environment and do a t-test like analysis, you want to treat each sample as though the count represents the number of observations of that particular transcript (transcript A) in the population of transcripts from that sample. The next sample from that same environment is like drawing a new set of observations from the same population of transcripts and observing how many are of the same transcript A. Did I understand that correctly?
Question: Is each sample from the same individual (i.e. a sample of the population of cells within that individual, thus technical replicates) or are they from different individuals in the same environment (i.e. a sample of each individual and thus representing a sample different populations of cells that have experienced the same environment)?
If the samples are from different individuals, I don't think this is valid approach. The counts are the response to some change in environment. It is like counting the number of photons of a particular wavelength from a light source in a test environment and then pooling it with the count of photons of the same wavelength from another, yet different, device of similar light source and also in the same environment. They are not the same light source (or device), but representatives of the same light source and must be treated independently of each other. This would represent device to device variability and must be tracked separately since each device is not an exact copy of the other. If they were, then multiple devices would not be needed for the test since all photons would come from a single source and each sample would be a representative sample of that source.
RNAseq is similar in that it is an effort of quantification of entity, or the abundance of a transcript, in a sample. The next sample is not typically a sample from the same source population. It is usually a sample from a different individual and thus must be treated that way. Individual to individual variability. Even samples from the same tissue from the same individual have variability (cells in different states of division, stratification of functionality, etc.) and should not be pooled as though they came from the same population of cells.
Ultimately, this all really depends on what the scientific question that is being asked and the experimental design used to discover the answer to that question. You might be able to pool them and do what you want or maybe not. It depends....
Replicates are from different individuals; biological replicates. And you are absolutely right, pooling them would not be a good idea. However, I can include replicate as a random independent variable into my model, nested within environment treatments, essentially treating this design as a block design - the same experiment repeated 3 times.
I am more worrying about the general approach - creating dummy observations with expected "counts" under null hypothesis (i.e., equal counts of reads everywhere) and using a contingency table /nominal logistics model to see if the actual counts differ.
I will, as another user suggested, post an example with the actual data.
