Wednesday, March 29, 2017

RNA-seq: How many replicates?

Do I need replicates for RNA-seq?
Short answer: yes.

Longer answer: Depends on  the goal of your analysis. For simplicity, let's break RNA-seq into two types of studies: Qualitative and Quantitative. An example of a qualitative analysis is a transcriptome assembly. For this kind of study sequencing depth is your friend, and paired-end sequencing a must. Replicates in this case would be helpful but optional. By contrast in a quantitative analysis, like a differential expression experiment, replicates are essential. These days many experimenters, with good reason, try and do both at the same time so the long answer is still 'yes'.

The number of replicates you decide to sequence is quite literally a cost-benefit analysis since the major limiter these days seems to be the price of sequencing. Two important things to consider are:

    1: How noisy is your data.
    2: What you’re looking for.

Let’s look (hear?) at noise. Below are two example sets both from the starlet sea anemone Nematostella vectensis. The top two plots are RNAseq data from embryonic samples and the bottom two from regeneration. In the Embryonic case, we see that the Biological coefficient of variation (which is the square root of the dispersion) is fairly low with the common BCV being ~0.18. In the principle component plot on the right, each replicate of the three samples cluster tightly together. We can conclude that our inter-replicate variation is low for these data. In the case of Regeneration (lower plots) the common BCV is higher ~0.21 and the replicates don't exhibit the same tightness in their clustering. This is because regeneration is not as synchronized or stereotypical as embryogenesis and this makes the data a bit more noisy. For this reason we opted for triplicates in our sampling which will increase our ability to detect differentially expressed genes later on.


The second thing to consider when deciding the number of replicates is the goals of the study. More replicates will yield more gene discovery in a simple differential expression experiment, although with diminishing returns (Liu et al. 2014; Schruch et al. 2016). Consider this figure from Liu et al.  As you can see the number of discovered DE genes increases with replicates while sequencing depth gives diminishing returns after 10M reads.




So if your goal is to find  subtle changes in gene expression in inherently noisy data, the more replicates the better. If however you just want a list of the most highly DE genes from samples that are fairly stereotypical, two full biological replicates may suffice.

One last quick note about transcriptome assembly. One does not need replicates per se, and a greater depth of sequencing will give you more detectable transcripts. However you will most likely want to quantify the results later on so might as well include those samples (and their replicates) in your experimental design.

Oh yeah and here's the code for the BCV plot and the PCA:

library(edgeR)
library(ggplot2)

#EdgeR and BCV plot 
m <- data.matrix(rawcounts)#raw counts is a count matrix with columns as samples and rows as genes
m <- m[rowSums(m > 5) >=2,] #drop low counts. for large datasets I set >=(0.5*sample number)

counts <-  m
group <- c(1,1,2,2,3,3) #designate replicates
BCV.true <- 0.1
y <- DGEList(counts=counts, group=group)
y <- estimateCommonDisp(y)
y <- estimateTrendedDisp(y)
y <- estimateTagwiseDisp(y)
plotBCV(y, main="Embryo", xlim=c(-5,15), ylim=c(0.0,0.8)) 
#I set the limits so the two plots look the same

#PCA

z <- cpm(y, normalized.lib.size=TRUE)
log_data = log(z+1,2) #log transforming is optional
condition <- factor(c(1,1,2,2,3,3))
replicate <- factor(c(1,2,1,2,1,2))
transpose <- t(log_data) #this puts samples as rows, genes as columns
transpose_df <- as.data.frame(transpose)
pca.data <- prcomp(transpose_df)
summary(pca.data)
scores = as.data.frame(pca.data$x) #get scores for plotting
p <- ggplot(data = scores, aes(x = PC1, y = PC2, label = rownames(scores), colour=factor(group), shape=factor(grouprep))) + 

geom_point(size=6) + scale_fill_hue(l=40) + 
coord_fixed(ratio=1, xlim=c(-100, 100), ylim=c(-100, 100)) +
labs(title= "Regeneration",color = "Sample", shape= "Replicate")
p

  1. Liu Y, Zhou J, White KP. RNA-seq differential expression studies: more
    sequence or more replication? Bioinformatics. 2014 Feb 1;30(3):301-4. doi:
    10.1093/bioinformatics/btt688. Epub 2013 Dec 6. PubMed PMID: 24319002; PubMed
    Central PMCID: PMC3904521.
  2. Schurch NJ, Schofield P, Gierliński M, et al. How many biological replicates are needed in an RNA-seq experiment and which differential expression tool should you use? RNA. 2016;22(6):839-851. doi:10.1261/rna.053959.115.


No comments:

Post a Comment

Clustering RNAseq data using K-means: how many clusters?

Clustering RNAseq data using K-means: how many clusters? Note this is part 2 of a ser...