differential gene expression complex design no replicates

Question

We have an experimental design as seen below

Where we administered drug at 0 min for each mouse genotype and took them down at given below intervals. wt and ko mouse models were administered only with the vehicle which would be our negative control.

Performed rnaseq (total rna) truseq stranded protocol would like to see differences between each time point compared to the vehicle, used star aligner to align to mm9 , htseq for gene counts .

sample  type    time
WT  Vehicle     30 min
KO  Vehicle     30 min
WT  Drug    30min
KO  Drug    30min
WT  Drug    1 hour
KO  Drug    1 hour
WT  Drug    2 hour
KO  Drug    2hour
WT  Drug    3 hour
KO  Drug    3hour

We do not have any replicates and hence cannot use any of the differential expression statistical tools (asked this question at bioconductor where Dr. Love himself replied link to question asked) where he suggested a linear model.

we were wondering if someone could help us with any such design (papers/previously performed analysis) what we would like to see is comparison/differences (significant) between these groups (veh vs 30min, veh vs 1 hour, beh vs 2 hours, veh vs 3 hours)

Answer 1

Although it is not recommended to use no replicates, in the edgeR manual they give some advice on how to go on with no replicates design. See page 21 of user guide. You can e.g. estimate a BCV value. Of course it is still a trick and not sound statistics.

Answer 2

Mike Love is right. If the response you are looking for is linear in terms of change per minute, the most productive approach is likely to be fitting a linear model. You might get something out of this because the difference between successive time points represents multiple measurements of the rate of change over time. The biggest problem is that the time points aren't independent. You would need to fit a model that estimated that rate for both WT and KO and then test the difference in rates.

I see three ways forward, in all cases you want to test the last term in the model ~time + genotype + time:genotype

1. Ideally you would like to fit a negative binomial linear model to the read counts. But that is going to be difficult as you have no replicates from which to assess the dispersion. Thus you could generate read counts using your favorite read count pipeline (STAR+featureCounts/Salmon/Kallisto) and take the approach in the edgeR guide noted by b.nota and then use edgeR to test the fitted model.

2. You could use something like salmon or kallisto to estimate the TPMs assume that log TPMs are normally distributed and model them with plain lm in R, using something like broom to model each gene separately. A tutorial by David Robinson, the author of broom on doing this sort of thing can be found here

3. You could use the Salmon or Kallisto bootstraps to estimate technical variation (but not biological) and then I suspect you could use Sleuth to fit the model.

In any case these results should be treated with scepticism andmostly treated as hypothesis generating. You might find genes that you want to go back and confirm and you might see some patterns of interest if you look at some enrichments etc. which you can then design experiments to test.