RNA-seq analysis in R

Linear Models

Original Authors: Belinda Phipson, Anna Trigos, Matt Ritchie, Maria Doyle, Harriet Dashnow, Charity Law

Based on the course RNAseq analysis in R delivered on May 11/12th 2016

Simulating linear models

We are going to simulate an experiment with 2 factors:

Treatment <- rep(c("Treatment A", "Control"), 3)
X <- cbind(c(1,0,1,0,1,0), c(0,1,0, 1, 0, 1))
Treatment <- factor(Treatment, levels=c("Treatment A", "Control"))
model.matrix(~Treatment)

  (Intercept) TreatmentControl
1           1                0
2           1                1
3           1                0
4           1                1
5           1                0
6           1                1
attr(,"assign")
[1] 0 1
attr(,"contrasts")
attr(,"contrasts")$Treatment
[1] "contr.treatment"

model.matrix(~Treatment -1)

  TreatmentTreatment A TreatmentControl
1                    1                0
2                    0                1
3                    1                0
4                    0                1
5                    1                0
6                    0                1
attr(,"assign")
[1] 1 1
attr(,"contrasts")
attr(,"contrasts")$Treatment
[1] "contr.treatment"

Treatment <- rep(c("Treatment A", "No Treatment"), 4)
ER <- rep(c("+", "-"), c(4,4))

We can obtain the design matrix for a model without and with interaction with the following commands:

model.matrix(~Treatment + ER)

  (Intercept) TreatmentTreatment A ER+
1           1                    1   1
2           1                    0   1
3           1                    1   1
4           1                    0   1
5           1                    1   0
6           1                    0   0
7           1                    1   0
8           1                    0   0
attr(,"assign")
[1] 0 1 2
attr(,"contrasts")
attr(,"contrasts")$Treatment
[1] "contr.treatment"

attr(,"contrasts")$ER
[1] "contr.treatment"

model.matrix(~Treatment * ER)

  (Intercept) TreatmentTreatment A ER+ TreatmentTreatment A:ER+
1           1                    1   1                        1
2           1                    0   1                        0
3           1                    1   1                        1
4           1                    0   1                        0
5           1                    1   0                        0
6           1                    0   0                        0
7           1                    1   0                        0
8           1                    0   0                        0
attr(,"assign")
[1] 0 1 2 3
attr(,"contrasts")
attr(,"contrasts")$Treatment
[1] "contr.treatment"

attr(,"contrasts")$ER
[1] "contr.treatment"

Now we will simulate a model without interaction. We need to provide a set of theoretical variable for our parameters, and a level of noise in our model:

params <- c(3, 2, -4)
names(params) <- c("Intercept", "TreatA", "ER+")
noise <- 1
X <- model.matrix(~Treatment + ER)
Y <- X %*% params
errors <- rnorm(length(Y), 0, 1)
Y <- Y + errors

Now we can compute the estimates of our parameters: We need to use specific functions for matrix algebra, like: Transpose: t(X) Multiply matrices: %*% Inverse of a matrix: solve(X)

beta.pars <- solve(t(X) %*% X) %*% t(X) %*% Y

However, we can use the linear model function in R, lm() and inspect the results of the fit, compute confidence intervals, etc.:

m <- lm(Y ~ Treatment + ER)
summary(m)

Call:
lm(formula = Y ~ Treatment + ER)

Residuals:
       1        2        3        4        5        6        7        8 
-1.09088  1.47036  0.03393 -0.41342  0.22952 -0.32015  0.82742 -0.73679 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)   
(Intercept)            3.2550     0.6066   5.366  0.00303 **
TreatmentTreatment A   1.8778     0.7005   2.681  0.04378 * 
ER+                   -3.8231     0.7005  -5.458  0.00281 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9906 on 5 degrees of freedom
Multiple R-squared:  0.8809,    Adjusted R-squared:  0.8332 
F-statistic: 18.49 on 2 and 5 DF,  p-value: 0.004897

confint(m)

                           2.5 %    97.5 %
(Intercept)           1.69558831  4.814381
TreatmentTreatment A  0.07715716  3.678429
ER+                  -5.62371468 -2.022443

Large Scale Hypothesis testing: FDR

When we are doing thousands of tests for differential expression, the overall significance level of a test is very difficult to control. Let’s see why: First, we simulate 40,000 genes not differentially expressed (with a mean of zero). We assume that we have 10 replicates of this experiment:

N <- 40000
R <- 10
X <- matrix(rnorm(N* R, 0, 1), nrow=N)

Now we assume that we run a t-test under the null hypothesis that the mean is zero for each of these genes, that is each row in the matrix:

t.test(X[1,])$p.value

[1] 0.5803333

pvals <- apply(X, 1, function(y) t.test(y)$p.value)

Because we have generated this data with mean zero, we know that none of these genes are differentially expressed, so we would like to be able to not reject any of the hypothesis. However, if you choose a significance level of 0.05 we get

sum(pvals<0.05)

[1] 2064

Too many rejections!!! In fact, if we look at the distributions of the p-values obtained we get:

hist(pvals)

That is, if the null hypothesis is true, the p-values will follow a uniform distribution. This is the key to all methods that aim to control the proportion of false positives amongs the genes that we call differentially expressed. Let’s add 1000 genes to our set that are really differentially expressed (mean of 1):

df <- 1000
Y <- matrix(rnorm(df* R, 1, 1), nrow=df)
Z <- rbind(X, Y)
pvals <- apply(Z, 1, function(y) t.test(y)$p.value)

Let’s look at the distribution of p-values now:

hist(pvals)

What would be the number of false positives now? How many would we expect if we reject p-values samller than our significance level, 0.05?

exp.sig<- (nrow(Z))*0.05
obs.sig <- sum(pvals<0.05)
FDR <- exp.sig / obs.sig
FDR

[1] 0.7167832

We can compare this with the Benjamini-Hochberg method:

pvals.adj <- p.adjust(pvals, method="BH")
plot(pvals, pvals.adj)
abline(v=0.05, col=2)