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
