Statistical analysis for experimental design with 3 factors from RNA-seq data

Question

[alfonseca]

Hi,

I'm analyzing RNA-seq data trying to find genes which RNA levels change depending on 3 factors of an experimental design and I couldn't find information about which could be the best method to do that (I have read that all these programs like edgeR, DEseq and limma just work with paired comparisons and I have not found good results with 3-way ANOVA)

In my experimental design I have 2 treatment conditions (treated and control), 4 times of treatment and 2 cellular fractions. I'm interested to find genes which transcripts have differences in their RNA levels between the cellular fractions during the time after treatment.

I will thank all kind of help

Thanks in advance

[Benn]

Try limma (voom or trend) for 3 factorial analysis. In the limma manual they describe how to do multi level experiments (9.7), in the same line you can make your design. If you need to test interaction, see how they do 2 x 2 factorial analysis in 9.5.

[alfonseca]

Thanks, I don't really understand how to use voom or trend. I have read the limma manual and we did the design table, but we don't understand how to test the interactions, according to what I understand it still continue doing pair analysis, according to the contrasts that I write. Probably I'm not understanding well. Do you know where I could read more about it? Thanks

[Benn]

If you can't find it in the manual, you can ask it to the authors of limma at bioconductor support site. They will answer you mostly within a day.

Answer 1

[Asaf]

You can use DESeq2 for the analysis although it will not treat the time points as linear on time, just four different conditions. I think that analyzing all the data together might be too complicated, I would take all the data and normalize it and then split it in different ways: take one cellular fraction and test the time:treatment interaction, take the treatment only and test time:fraction interaction, take time 0 and test treatment:fraction interaction etc. There are packages that deal with the time factor better but I'm not experienced with them.

[alfonseca]

Thanks, I'm not analyzing the time as a continous component. My problem with using DESeq2 is that I finish with a lot of lists, where each tells me a different story. I would like to do a more global analysis first and then start to analyze for time or fraction.

[Asaf]

I would start with plotting all the data with PCA. You can definitely analyze everything together and ask the one question you're interested in as an interaction which will probably be time:treatment:fraction. I use this code to plot PCA, it's in their manual:

dds <- DESeqDataSetFromMatrix(countData = counts,
                          colData = colData,
                          design = ~ Treatment + Time + Fraction + Treatment:Time:Fraction)
dds <- estimateSizeFactors(dds)
dds <- estimateDispersions(dds)

vsd = vst(dds, blind=FALSE)
pcaData <- plotPCA(vsd, intgroup=c("Fraction", "Treatment", "Time"), returnData=TRUE)
percentVar <- round(100 * attr(pcaData, "percentVar"))
gplt = ggplot(pcaData, aes(PC1, PC2, color=Time, shape=Treatment:Fraction)) +
  geom_point(size=3) + 
  xlab(paste0("PC1: ",percentVar[1],"% variance")) +
  ylab(paste0("PC2: ",percentVar[2],"% variance")) + 
  ggtitle("All data PCA") +
  coord_fixed()
print(gplt)

[alfonseca]

Thanks Asaf. I've already checked the PCA and I can see that the main components are cellular fraction and treatment (in that order). But I don't understand how the PCA can give me information about the interaction, and neither about the genes that change their levels depending on the 3 factors.

[Asaf]

It can't, it can just give you an overview