# differential gene expression complex design no replicates

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)

• What is your design, exactly? What are we supposed to understand from that vague table? What are the times? The time you administered the drug? The length of time the drug had to act before administering an antagonist? The time you spent observing it? How is this related to bioinformatics? Please edit your question and explain what you are doing and what you expect from a good answer. – terdon Jul 25 '17 at 17:45
• It's unlikely that you'll find any papers doing this because it'd be hard to get one accepted if it were largely based on such an experiment. I have to wonder why you expended time/resources on this experiment when you had to know that this was going to be a major problem. – Devon Ryan Jul 25 '17 at 18:01
• @DevonRyan we have experimental (wet lab results supporting) this, realized it will be difficult to go ahead with rna seq only after i posted it on deseq forum (bioconductor), we are not sure if this is a dead end or we can at least see some basic differences using very basic stats – novicebioinforesearcher Jul 25 '17 at 18:53
• @terdon thank you for feedback i have edited the question, please let me know if further information is required – novicebioinforesearcher Jul 25 '17 at 18:56
• So you administered the drug and then killed the mice after the specified time. Then what? Did you perform rna seq? Of what? Whole transcriptome? What are you trying to show? What data do you have? Are you really trying to draw conclusions based on a single individual in each condition? – terdon Jul 25 '17 at 19:23

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.

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.

• thank you for the information is there a tutorial for the lm and broom methods that you have mentioned would help me understand it better, as R is something I have just started using appreciate it. – novicebioinforesearcher Jul 26 '17 at 13:53
• I have added a link to a broom tutorial. – Ian Sudbery Jul 26 '17 at 14:29
• it is really an awesome resource many thanks!! (I would not have found this), and keep everyone posted hopefully i can pick up on using R and work this out within a couple of days. – novicebioinforesearcher Jul 26 '17 at 14:53