RNA-seq analysis in R

Differential Expression of RNA-seq data

Recap of pre-processing

The previous section walked-through the pre-processing and transformation of the count data. Here, for completeness, we list the minimal steps required to process the data prior to differential expression analysis.

# Read the sample information into a data frame
sampleinfo <- read_tsv("data/SampleInfo.txt")

# Read the data into R
seqdata <- read_tsv("data/GSE60450_Lactation.featureCounts", comment = "#")

# Transform the data to matrix of counts
countdata <- as.data.frame(seqdata) %>% 
    column_to_rownames("Geneid") %>% # turn the geneid column into rownames
    rename_all(str_remove, ".bam") %>% # remove the ".bam" from the column names
    select(sampleinfo$Sample) %>% # keep sample columns using sampleinfo
    as.matrix()

# filter the data to remove genes with few counts
keep <- rowSums(countdata) > 5
countdata <- countdata[keep,]

Load the data

Alternatively we can load the `objects with the RData file we created in the pre-processing tutorial.

# load the RData object we created in the previous session
load("Robjects/preprocessing.RData")
ls()

 [1] "col.cell"            "countdata"           "countVar"           
 [4] "ddsObj"              "design"              "highVar"            
 [7] "hmDat"               "keep"                "librarySizes"       
[10] "logcounts"           "logNormalizedCounts" "morecols"           
[13] "mypalette"           "normalizedCounts"    "pcDat"              
[16] "rlogcounts"          "sampleinfo"          "seqdata"            
[19] "statusCol"          

dim(countdata)

[1] 22013    12

sampleinfo

The model formula and design matrices

First load the packages we need.

library(tidyverse)
library(DESeq2)

Now that we are happy that that the quality of the data looks good, we can proceed to testing for differentially expressed genes. There are a number of packages to analyse RNA-Seq data. Most people use DESeq2 or edgeR. They are both equally applicable. There is an informative and honest blog post here by Mike Love, one of the authors of DESeq2, about deciding which to use.

We will use DESeq2 for the rest of this practical.

Create a DESeqDataSet object with the raw data

Creating the design model formula

First we need to create a design model formula for our analysis. DESeq2 will use this to generate the model matrix, as we have seen in the linear models lecture.

We have two variables: “status”" and “cell type”. We will fit two models under two assumptions: no interaction and interaction of these two factors.

Let’s start with the model with only main effects, that is no interaction. The main assumption here is that the effect of the status is the same in both type of cells.

# Use the standard R 'formula' syntax for an additive model
design <- as.formula(~ CellType + Status)

What does this look like as a model matrix?

modelMatrix <- model.matrix(design, data = sampleinfo)
modelMatrix

   (Intercept) CellTypeluminal Statuspregnant Statusvirgin
1            1               0              0            1
2            1               0              0            1
3            1               0              1            0
4            1               0              1            0
5            1               0              0            0
6            1               0              0            0
7            1               1              0            1
8            1               1              0            1
9            1               1              1            0
10           1               1              1            0
11           1               1              0            0
12           1               1              0            0
attr(,"assign")
[1] 0 1 2 2
attr(,"contrasts")
attr(,"contrasts")$CellType
[1] "contr.treatment"

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

It would be nice if virgin were the the base line/intercept. To get R to use virgin as the intercept we need to use a factor. Let’s set factor levels on Status to use virgin as the intercept.

sampleinfo$Status <- factor(sampleinfo$Status, 
                              levels = c("virgin", "pregnant", "lactate"))
modelMatrix <- model.matrix(design, data = sampleinfo)
modelMatrix

   (Intercept) CellTypeluminal Statuspregnant Statuslactate
1            1               0              0             0
2            1               0              0             0
3            1               0              1             0
4            1               0              1             0
5            1               0              0             1
6            1               0              0             1
7            1               1              0             0
8            1               1              0             0
9            1               1              1             0
10           1               1              1             0
11           1               1              0             1
12           1               1              0             1
attr(,"assign")
[1] 0 1 2 2
attr(,"contrasts")
attr(,"contrasts")$CellType
[1] "contr.treatment"

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

Build a DESeq2DataSet

We don’t actually need to pass DESeq2 the model matrix, instead we pass it the design formula and the sampleinfo it will build the matrix itself.

# create the DESeqDataSet object
ddsObj <- DESeqDataSetFromMatrix(countData = countdata,
                              colData = sampleinfo,
                              design = design)

some variables in design formula are characters, converting to factors

Data exploration

Let’s plot a PCA from vst transformed data. Can you anticipate if the interaction term will be important?

vstcounts <- vst(ddsObj, blind=TRUE)
plotPCA(vstcounts, intgroup=c("Status", "CellType"))

Differential expression analysis with DESeq2

The DESeq2 work flow

The main DESeq2 work flow is carried out in 3 steps:

First, Calculate the “median ratio” normalisation size factors…

ddsObj <- estimateSizeFactors(ddsObj)

… then estimate dispersion …

ddsObj <- estimateDispersions(ddsObj)

gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates

… finally, apply Negative Binomial GLM fitting and calculate Wald statistics

ddsObj <- nbinomWaldTest(ddsObj)

The DESeq command

In practice the 3 steps above can be performed in a single step using the DESeq wrapper function. Performing the three steps separately is useful if you wish to alter the default parameters of one or more steps, otherwise the DESeq function is fine.

# rebuild a clean DDS object
ddsObj <- DESeqDataSetFromMatrix(countData = countdata,
                              colData = sampleinfo,
                              design = design)

some variables in design formula are characters, converting to factors

# Run DESeq
ddsObj <- DESeq(ddsObj)

estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing

Generate a results table

We can generate a table of differential expression results from the DDS object using the results function of DESeq2.

res <- results(ddsObj, alpha=0.05)
head(res)

log2 fold change (MLE): Status lactate vs virgin 
Wald test p-value: Status lactate vs virgin 
DataFrame with 6 rows and 6 columns
                            baseMean     log2FoldChange             lfcSE
                           <numeric>          <numeric>         <numeric>
ENSMUSG00000051951  193.628316332322  0.685666258943773 0.756712220191988
ENSMUSG00000102331 0.552133378452342   1.82614307383859  2.73990068529035
ENSMUSG00000025900  2.09442070179461  -3.02782655402946  1.54150999536643
ENSMUSG00000025902  52.4208304161579 -0.703252517806911  0.41735658148811
ENSMUSG00000098104 0.684213614702782  0.201738946689915   1.6866157166625
ENSMUSG00000103922  27.5871076674379   1.63592009157737 0.669531803485605
                                stat             pvalue               padj
                           <numeric>          <numeric>          <numeric>
ENSMUSG00000051951 0.906112311454162  0.364876409573963  0.520197122672024
ENSMUSG00000102331 0.666499732505843  0.505091734597736                 NA
ENSMUSG00000025900 -1.96419521321996 0.0495074577247786   0.11736643706379
ENSMUSG00000025902 -1.68501600070477  0.091985484452551   0.18995625581263
ENSMUSG00000098104 0.119611684331461  0.904790763422636                 NA
ENSMUSG00000103922  2.44337921374416 0.0145504392990467 0.0446034399346747

Independent filtering

You will notice that some of the adjusted p-values (padj) are NA. Remember in Session 2 we said that there is no need to pre-filter the genes as DESeq2 will do this through a process it calls ‘independent filtering’. The genes with NA are the ones DESeq2 has filtered out.

From DESeq2 manual: “The results function of the DESeq2 package performs independent filtering by default using the mean of normalized counts as a filter statistic. A threshold on the filter statistic is found which optimizes the number of adjusted p values lower than a [specified] significance level”.

The default significance level for independent filtering is 0.1, however, you should set this to the FDR cut off you are planning to use. We will use 0.05 - this was the purpose of the alpha argument in the previous command.

The default contrast of results

The results function has returned the results for the contrast “lactate vs virgin”. Let’s have a look at the model matrix to understand why DESeq2 has given us this particular contrast.

modelMatrix

   (Intercept) CellTypeluminal Statuspregnant Statuslactate
1            1               0              0             0
2            1               0              0             0
3            1               0              1             0
4            1               0              1             0
5            1               0              0             1
6            1               0              0             1
7            1               1              0             0
8            1               1              0             0
9            1               1              1             0
10           1               1              1             0
11           1               1              0             1
12           1               1              0             1
attr(,"assign")
[1] 0 1 2 2
attr(,"contrasts")
attr(,"contrasts")$CellType
[1] "contr.treatment"

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

By default, results has returned the contrast encoded by the final column in the model matrix. DESeq2 has the command resultsNames that allows us to view the contrasts that are available directly from the model matrix.

resultsNames(ddsObj)

[1] "Intercept"                 "CellType_luminal_vs_basal"
[3] "Status_pregnant_vs_virgin" "Status_lactate_vs_virgin" 

Let’s just rename res so that we know which contrast results it contains.

resLvV <- res
rm(res)

If we want a different contrast we can just pass the results function the name of the design matrix column that encodes it. Let’s retrieve the results for pregant versus virgin

resPvV <- results(ddsObj, 
                  name="Status_pregnant_vs_virgin", 
                  alpha = 0.05)
resPvV

log2 fold change (MLE): Status pregnant vs virgin 
Wald test p-value: Status pregnant vs virgin 
DataFrame with 22013 rows and 6 columns
                            baseMean     log2FoldChange             lfcSE
                           <numeric>          <numeric>         <numeric>
ENSMUSG00000051951  193.628316332322 -0.954402462919324 0.758163669842148
ENSMUSG00000102331 0.552133378452342  0.793605030301095  2.76321615795618
ENSMUSG00000025900  2.09442070179461 -0.422312119522277  1.28960999953541
ENSMUSG00000025902  52.4208304161579 -0.372373337372372 0.399035417679429
ENSMUSG00000098104 0.684213614702782   -1.4747163554633  1.79948663916514
...                              ...                ...               ...
ENSMUSG00000064367  54055.8663105821  0.344588799853029  0.16434684684871
ENSMUSG00000064368  13454.3113895505  0.426945077829144 0.156109922740747
ENSMUSG00000064370  111885.216173771  0.415862003187522 0.225727285622426
ENSMUSG00000064372   396.00873428211  0.262372947574112 0.261696769419213
ENSMUSG00000095742  632.014515495881 -0.191806454656629 0.355891725136255
                                 stat              pvalue               padj
                            <numeric>           <numeric>          <numeric>
ENSMUSG00000051951  -1.25883433997574   0.208090174010475  0.411522854118238
ENSMUSG00000102331  0.287203383642663   0.773956595641228                 NA
ENSMUSG00000025900 -0.327472739568099   0.743310358007849                 NA
ENSMUSG00000025902 -0.933183674616883    0.35072514229176   0.56448494555843
ENSMUSG00000098104 -0.819520591799159   0.412489458422432                 NA
...                               ...                 ...                ...
ENSMUSG00000064367   2.09671683065658  0.0360186497109148  0.124182953552384
ENSMUSG00000064368    2.7349003210909 0.00623991686862116 0.0341430901519056
ENSMUSG00000064370   1.84232048881823  0.0654282845127896  0.190916096241363
ENSMUSG00000064372   1.00258382308807   0.316061704202673  0.530023878510884
ENSMUSG00000095742 -0.538946092616217   0.589924051629706  0.762012734103323

Let’s get the top 100 genes by adjusted p-value

topGenesPvV <- as.data.frame(resPvV) %>%
    rownames_to_column("GeneID") %>% 
    arrange(padj) %>% 
    head(100)
topGenesPvV

Challenge 1

Obtain results for luminal vs basal and find the top 200 genes. Call the new results object resBvL.

Contrasts

Suppose we want to find differentially expressed genes between pregnant and lactate. We don’t have a parameter that explicitly will allow us to test that hypothesis. We need to provide a contrast.

resultsNames(ddsObj)

[1] "Intercept"                 "CellType_luminal_vs_basal"
[3] "Status_pregnant_vs_virgin" "Status_lactate_vs_virgin" 

resPvL <- results(ddsObj,
                  contrast=c("Status", "pregnant", "lactate"), 
                  alpha = 0.05)
resPvL

log2 fold change (MLE): Status pregnant vs lactate 
Wald test p-value: Status pregnant vs lactate 
DataFrame with 22013 rows and 6 columns
                            baseMean      log2FoldChange             lfcSE
                           <numeric>           <numeric>         <numeric>
ENSMUSG00000051951  193.628316332322    -1.6400687218631 0.760113817008618
ENSMUSG00000102331 0.552133378452342    -1.0325380435375  2.70726895585512
ENSMUSG00000025900  2.09442070179461    2.60551443450719   1.5811323316226
ENSMUSG00000025902  52.4208304161579   0.330879180434539 0.424588010118588
ENSMUSG00000098104 0.684213614702782   -1.67645530215321  1.83315246550806
...                              ...                 ...               ...
ENSMUSG00000064367  54055.8663105821  -0.302824789042762 0.164348578199074
ENSMUSG00000064368  13454.3113895505  -0.251495130843604 0.156113680276774
ENSMUSG00000064370  111885.216173771  -0.196775657448966 0.225728145700262
ENSMUSG00000064372   396.00873428211 -0.0171568600400927 0.262491263102943
ENSMUSG00000095742  632.014515495881   0.280640913843395 0.357511558787336
                                  stat             pvalue              padj
                             <numeric>          <numeric>         <numeric>
ENSMUSG00000051951   -2.15766203055943 0.0309541180695016 0.106766273830346
ENSMUSG00000102331  -0.381394704543258  0.702910389783292                NA
ENSMUSG00000025900    1.64787879065969 0.0993775464915822                NA
ENSMUSG00000025902   0.779294686965191  0.435806143702786 0.636130155679643
ENSMUSG00000098104  -0.914520386981878  0.360443470414368                NA
...                                ...                ...               ...
ENSMUSG00000064367   -1.84257626297171 0.0653909022816961 0.183136793003594
ENSMUSG00000064368   -1.61097432587412  0.107185318448375 0.258163585174442
ENSMUSG00000064370  -0.871737358398624  0.383351674854833 0.588905143355555
ENSMUSG00000064372 -0.0653616422782202  0.947886063772367 0.973183515620879
ENSMUSG00000095742   0.784984168890418  0.432462850616363 0.633552881519567

Comparing two design models

Suppose we thought that maybe status were irrelevant and really the only differences might be between cell types. We could fit a simpler model, this would give us more degrees of freedom and therefore more power, but how would we know if it was a better model of not? We can compare the two models using the “log ratio test” (LRT).

designC <- as.formula(~ CellType )
# Compare the designs
ddsObjC <- DESeq(ddsObj, test="LRT", reduced=designC)

using pre-existing size factors
estimating dispersions
found already estimated dispersions, replacing these
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing

resCvCS <- results(ddsObjC)
resCvCS

log2 fold change (MLE): Status lactate vs virgin 
LRT p-value: '~ CellType + Status' vs '~ CellType' 
DataFrame with 22013 rows and 6 columns
                            baseMean     log2FoldChange             lfcSE
                           <numeric>          <numeric>         <numeric>
ENSMUSG00000051951  193.628316332322  0.685666258943773 0.756712220191988
ENSMUSG00000102331 0.552133378452342   1.82614307383859  2.73990068529035
ENSMUSG00000025900  2.09442070179461  -3.02782655402946  1.54150999536643
ENSMUSG00000025902  52.4208304161579 -0.703252517806911  0.41735658148811
ENSMUSG00000098104 0.684213614702782  0.201738946689915   1.6866157166625
...                              ...                ...               ...
ENSMUSG00000064367  54055.8663105821  0.647413588895792 0.164350605694756
ENSMUSG00000064368  13454.3113895505  0.678440208672748 0.156129934723006
ENSMUSG00000064370  111885.216173771  0.612637660636488 0.225729381045944
ENSMUSG00000064372   396.00873428211  0.279529807614205 0.262569324325006
ENSMUSG00000095742  632.014515495881 -0.472447368500024 0.357092689972106
                                stat               pvalue
                           <numeric>            <numeric>
ENSMUSG00000051951  4.40903494919449    0.110303736820533
ENSMUSG00000102331 0.397659679654588    0.819689359949812
ENSMUSG00000025900  4.68578354199757   0.0960494831341325
ENSMUSG00000025902  2.89099132419746    0.235629255115163
ENSMUSG00000098104   1.1310712566648    0.568055798475988
...                              ...                  ...
ENSMUSG00000064367  15.2088899461904 0.000498231876872594
ENSMUSG00000064368   18.663347272815 8.85738725187499e-05
ENSMUSG00000064370  7.28491634377798   0.0261878904431057
ENSMUSG00000064372  1.38080328700747    0.501374654747371
ENSMUSG00000095742  1.63512832056145    0.441505783095116
                                   padj
                              <numeric>
ENSMUSG00000051951    0.203011375932151
ENSMUSG00000102331                   NA
ENSMUSG00000025900    0.182314817582248
ENSMUSG00000025902    0.361374169977104
ENSMUSG00000098104                   NA
...                                 ...
ENSMUSG00000064367  0.00246460256305281
ENSMUSG00000064368 0.000552135624072482
ENSMUSG00000064370   0.0653923895104458
ENSMUSG00000064372    0.622036554157213
ENSMUSG00000095742    0.569776660112915

The null hypothesis is that there is no significant difference between the two models, i.e. the simpler model is sufficient to explain the variation in gene expression between the samples. If the the adjusted p-value for a gene passes a significance threshold (e.g. padj < 0.05) then we should consider using the more complex model for this gene. In practice we would usually choose one model or the other and apply it to all genes.

Challenge 2

When we looked at the PCA it did seem that an interaction model might be warranted. Let’s test that.
1.Fit a model with interaction. 2. Use the LRT to compare the two models.
3. Is the number of replicates good enough to include the interaction?
4. Is the interaction needed in the model?

Testing log2 fold change versus a threshold

A common practice when considering the results of a differential expression analysis is to filter out genes that are statistically significant but have a low fold change. Perhaps you are only interested in very strong response to the experimental conditions and decide to eliminate from consideration all genes with an absolute fold change lower than 2x.

Rather than do this by filtering after running the differential analysis, DESeq2 provides a means for testing the hypothesis that fold change is greater (or actually lesser than if you want) a given threshold.

resPvL2 <- results(ddsObj,
                  contrast=c("Status", "pregnant", "lactate"), 
                  alpha = 0.05,
                  lfcThreshold=0.5, 
                  altHypothesis="greaterAbs")
sum(resPvL2$padj<0.05, na.rm = T)

[1] 352

sum(resPvL$padj<0.05 & abs(resPvL$log2FoldChange)>=2^0.5, na.rm = T)

[1] 824

There are four possible values for altHypothesis:

• greaterAbs - |β|>x - tests are two-tailed
• lessAbs - |β|<x - p values are the maximum of the upper and lower tests
• greater - β>x
• less - β<x

Finally save the results in a new RData object

save(resLvV, ddsObj, sampleinfo, file="results/DE.RData")

---

References