Analysis of RNA-seq data in R

Learning Objectives

• Which statistical tests are appropriate for RNA-seq data
• Using the DESeq2 package to detect differential expression
• Basic visualisation of RNA-seq counts
• Using annotation databases to map between gene identifers

Differential expression with DESeq2

Now that we are happy that we have normalised the data and that the quality looks good, we can continue to testing for differentially expressed genes. There are a number of packages to analyse RNA-Seq data. Most people use DESeq2 or edgeR. We will use DESeq2 for the rest of this practical.

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.

Note that although we spent some time looking at the quality of our data , these steps are not required prior to performing differential expression so are not shown here. Remember, DESeq2 requires raw counts so the vst transformation is not shown as part of this basic protocol.

library(tximport)
library(DESeq2)
library(readr)
dirs <- list.files("salmon_quant/")
quant_files <- list.files("salmon_quant/",
                          pattern="quant.sf.gz",
                          recursive = TRUE,
                          full.names = TRUE)
names(quant_files) <- dirs

tx2gene <- read_csv("tx2gene.csv",col_names = FALSE)

txi <- tximport(quant_files,
                type="salmon",
                tx2gene = tx2gene)

sampleinfo <- read_tsv("meta_data/sampleInfo_corrected.txt")


dds <- DESeqDataSetFromTximport(txi, 
                                colData = sampleinfo,
                                design = ~Treated)
dds$condition <- as.factor(dds$condition)

It would be a good idea to save the results of the pre-processing so we don’t have to repeat it every time.

dir.create("Robjects/",showWarnings = FALSE)
saveRDS(dds, file="Robjects/dds.rds")

We will be using these raw counts throughout the workshop and transforming them using methods in the DESeq2 package. If you want to know about alternative methods for count normalisation they are covered on this page.

If you have problems running these steps, you can re-load a pre-processed object from the course materials

dds <- readRDS("Robjects/dds_BACKUP.rds")

The DESeq workflow in brief

We have previously defined the test condition using the design argument when we created the object. This can be checked using the design function.

Typically we decide the design for the analysis when we create the DESeq2 objects, but it can be modified prior to the differential expression analysis. The design tells DESeq2 which sample groups to compare in the differential analysis. The name specified must correspond to a column in the sample information.

colData(dds)

DataFrame with 9 rows and 5 columns
                    Run condition        Name Replicate
            <character>  <factor> <character> <numeric>
1_CTR_BC_2   1_CTR_BC_2       CTR       CTR_1         1
2_TGF_BC_4   2_TGF_BC_4       TGF       TGF_1         1
3_IR_BC_5     3_IR_BC_5       IR         IR_1         1
4_CTR_BC_6   4_CTR_BC_6       CTR       CTR_2         2
5_TGF_BC_7   5_TGF_BC_7       TGF       TGF_2         2
6_IR_BC_12   6_IR_BC_12       IR         IR_2         2
7_CTR_BC_13 7_CTR_BC_13       CTR       CTR_3         3
8_TGF_BC_14 8_TGF_BC_14       TGF       TGF_3         3
9_IR_BC_15   9_IR_BC_15       IR         IR_3         3
             Treated
            <factor>
1_CTR_BC_2         N
2_TGF_BC_4         Y
3_IR_BC_5          Y
4_CTR_BC_6         N
5_TGF_BC_7         Y
6_IR_BC_12         Y
7_CTR_BC_13        N
8_TGF_BC_14        Y
9_IR_BC_15         Y

design(dds) <- ~Treated

The counts that we have obtained via sequencing are subject to random sources of variation. The purpose of differential expression is to determine if potential sources of biological variation (e.g. counts observed from different sample groups) are greater than random noise.

The DESeq function runs a couple of processing steps automatically to adjust for different library size and gene-wise variability, which you can read about in the DESeq2 vignette and run some example code at the end of this session.

de_treated<- DESeq(dds)
de_treated

class: DESeqDataSet 
dim: 57914 9 
metadata(1): version
assays(6): counts avgTxLength ... H cooks
rownames(57914): ENSG00000000003 ENSG00000000005 ...
  ENSG00000284747 ENSG00000284748
rowData names(22): baseMean baseVar ... deviance
  maxCooks
colnames(9): 1_CTR_BC_2 2_TGF_BC_4 ... 8_TGF_BC_14
  9_IR_BC_15
colData names(5): Run condition Name Replicate Treated

The results of the analysis are not immediately accessible, but can be obtained using the results function. Each row is a particular gene measured in the study (i.e. all genes in the organism being studied) and each column reports some aspect of the differential expression analysis for that gene. Note that all genes are reported. At this stage the gene identifiers are not very informative, something we will fix in the next section. Furthermore, the results function displays results in a format which is not compatible with standard data manipulation tools (i.e. tidyverse), so we will have to convert.

results(de_treated)

log2 fold change (MLE): Treated Y vs N 
Wald test p-value: Treated Y vs N 
DataFrame with 57914 rows and 6 columns
                   baseMean log2FoldChange     lfcSE      stat
                  <numeric>      <numeric> <numeric> <numeric>
ENSG00000000003 1447.265802     -0.2051187  0.121418 -1.689365
ENSG00000000005    0.113575      0.5069982  3.533803  0.143471
ENSG00000000419 1814.108471      0.0773555  0.189163  0.408936
ENSG00000000457  645.735472     -0.1603703  0.124985 -1.283112
ENSG00000000460  221.065743     -0.4662569  0.423195 -1.101753
...                     ...            ...       ...       ...
ENSG00000284744    8.377622       0.381924  0.672894  0.567584
ENSG00000284745    0.000000             NA        NA        NA
ENSG00000284746    0.102422      -0.935671  3.533803 -0.264777
ENSG00000284747   29.068925      -0.333544  0.374367 -0.890956
ENSG00000284748    0.548454       2.181462  3.495590  0.624061
                   pvalue      padj
                <numeric> <numeric>
ENSG00000000003 0.0911494  0.574408
ENSG00000000005 0.8859182        NA
ENSG00000000419 0.6825867  0.950327
ENSG00000000457 0.1994530  0.740864
ENSG00000000460 0.2705689  0.800604
...                   ...       ...
ENSG00000284744  0.570318  0.924953
ENSG00000284745        NA        NA
ENSG00000284746  0.791181        NA
ENSG00000284747  0.372953  0.861365
ENSG00000284748  0.532587        NA

Processing the DE results using tidyverse

The output can be converted into a data frame and manipulated in the usual manner. It is recommended to use dplyr to manipulate the data frames with the standard set of operations detailed on the dplyr cheatsheet

• select to pick which columns to display
• filter to restrict the rows
• mutate to add new variables to the data frame
• arrange to order the data frame according to values of a column

The %>% symbol refers to the piping operation in R, which is a way of chaining operations together.

library(dplyr)
 results(de_treated, tidy=TRUE)

We can sort the rows by adjusted p-value and then print the first 10 rows.

results(de_treated,tidy=TRUE) %>%
  arrange(padj) %>%  
  head(n=10)

Or we can sort the rows and then write the resulting data frame to a file.

dir.create("de_analysis",showWarnings = FALSE)
 results(de_treated,tidy=TRUE) %>%
  arrange(padj) %>% 
   write_csv("de_analysis/treated_Y_vs_N_DESeq_all.csv")

Filtering to the differentially-expressed genes can be achieved using the filter function from dplyr.

 results(de_treated,tidy=TRUE) %>%
  filter(padj < 0.05) %>% 
  write.csv("de_analysis/treated_Y_vs_N_DESeq_sig.csv")

It is also a good idea to save the results object itself so we can re-use later.

saveRDS(de_treated, file="Robjects/de_treated.rds")

We can discover how many differentially-expressed genes (at a particular p-value cut-off) using the count function

results(de_treated,tidy=TRUE) %>%
count(padj < 0.05)

Another overview of the results is to use the plotMA function. Each point on this plot represents and individual gene with the x- and y-axes being the overall expression level and magnitude of difference respectively. Significant genes are automatically highlighted. The fanning effect at low expression levels is often seen due to high relative fold-change at low expression levels.

plotMA(de_treated)

It is also instructive to perform a “sanity” check and plot the sample-level counts for genes with high significance. This could highlight any other technical factors that we are not currently taking into account. The plot is not particularly attractive, but is a good quick diagnostic.

plotCounts(dds,"ENSG00000158258",intgroup = "Treated")

plotCounts(dds,"ENSG00000136999",intgroup = "Treated")

If your study involves knocking-out a particular gene, or you have some positive controls that are known in advance, it would be a good idea to visualise their expression level with plotCounts.

Exercise

• Re-run the analysis to find differentially-expressed genes between the TGF treated samples and CTRL
• Write a csv file that contains results for the genes that have an adjusted p-value less than 0.05 and a log2 fold change more than 1, or less than -1 in the contrast of TGF vs CTRL.
  • HINT: So that we don’t overwrite our results so far, it may be convenient to create a new DESeqDataSet object for the new differential expression analysis. Check the colData to see which analyses can be made
• Use the plotCounts function to visually-inspect the most statistically-significant gene identified

dds_condition<- dds
colData(dds_condition)
## You will need to change this line to choose the correct comparison
design(dds_condition) <- ~...

Changing the direction of the contrast

In this initial analysis DESeq2 has automatically decided which member of our sample groups to use as our baseline (CTR in this case). If the log2 fold changes has a positive value this means higher expression in Treated sample. We can alter this by changing the contrast argument in the results function

## This should give the same as the table above
results(de_treated, contrast=c("Treated","Y","N"))
## Changing the direction of the contrast
results(de_treated, contrast=c("Treated","N","Y"))

If we want to perform differential expression analysis on the condition variable then there are various contrasts that can be made; IR vs CTR, TGF vs CTR etc. When the results function is run with no contrast argument specified, the table that is displayed is for the contrast TGF vs CTR. The resultsNames function can tell us which other contrasts we can access.

dds_condition<- dds
design(dds_condition) <- ~condition
de_condition <- DESeq(dds_condition)
results(de_condition, contrast = c("condition","IR","CTR"),tidy=TRUE)  %>%
  arrange(padj) 

results(de_condition, contrast = c("condition","IR","TGF"),tidy=TRUE) %>%
  arrange(padj) 

More complex designs

The examples we have used so far have performed a differential expression analysis using a named column in the colData object. The DESeq2 package is capable of performing more complex analyses that can take multiple factors into consideration at the same time; so-called “multi-factor designs”

• Multi-factor designs in DESeq2

The use of such a design could be motivated by discovering sources of technical variation in our data that might obscure the biological differences we would like to compare. e.g.

In the example image above the main source of variation is the batch in which the samples were sequenced. A multi-factor analysis to compare the various conditions, but “correct” for differences in batch, would be as follows.

### Don't run this. It's just a code example
design(MY_DATA) <- ~ batch + condition

Likewise, if we have different treatments applied to difference cell-lines, but the main source of variation is the cell line the following could be used.

### Don't run this. It's just a code example
design(MY_DATA) <- ~cell_line + treatment

Adding annotation to the DESeq2 results

We would love to share these results with our collaborators, or search for our favourite gene in the results. However, the results are not very useful in there current form as each row is named according to an Ensembl identifier. Whilst gene symbols are problematic and can change over time, they are the names that are most recognisable and make the results easier to navigate.

There are a number of ways to add annotation, but we will demonstrate how to do this using the org.Hs.eg.db package. This package is one of several organism-level packages in Bioconductor that are re-built every 6 months. These packages are listed on the annotation section of the Bioconductor, and are installed in the same way as regular Bioconductor packages. An alternative approach is to use biomaRt, an interface to the BioMart resource. BioMart is much more comprehensive, but the organism packages do not require online access once downloaded.

### Only execute when you need to install the package
install.packages("BiocManager")
BiocManager::install("org.Hs.eg.db")
# For Human
BiocManager::install("org.Hs.eg.db")

The packages are larger in size that Bioconductor software packages, but essentially they are databases that can be used to make offline queries. An alternatve (biomaRt) that connects to the ensembl biomart resource will be discussed later.

library(org.Hs.eg.db)

First we need to decide what information we want. In order to see what we can extract we can run the columns function on the annotation database.

columns(org.Hs.eg.db)

 [1] "ACCNUM"       "ALIAS"        "ENSEMBL"     
 [4] "ENSEMBLPROT"  "ENSEMBLTRANS" "ENTREZID"    
 [7] "ENZYME"       "EVIDENCE"     "EVIDENCEALL" 
[10] "GENENAME"     "GENETYPE"     "GO"          
[13] "GOALL"        "IPI"          "MAP"         
[16] "OMIM"         "ONTOLOGY"     "ONTOLOGYALL" 
[19] "PATH"         "PFAM"         "PMID"        
[22] "PROSITE"      "REFSEQ"       "SYMBOL"      
[25] "UCSCKG"       "UNIPROT"     

We are going to filter the database by a key or set of keys in order to extract the information we want. Valid names for the key can be retrieved with the keytypes function.

keytypes(org.Hs.eg.db)

 [1] "ACCNUM"       "ALIAS"        "ENSEMBL"     
 [4] "ENSEMBLPROT"  "ENSEMBLTRANS" "ENTREZID"    
 [7] "ENZYME"       "EVIDENCE"     "EVIDENCEALL" 
[10] "GENENAME"     "GENETYPE"     "GO"          
[13] "GOALL"        "IPI"          "MAP"         
[16] "OMIM"         "ONTOLOGY"     "ONTOLOGYALL" 
[19] "PATH"         "PFAM"         "PMID"        
[22] "PROSITE"      "REFSEQ"       "SYMBOL"      
[25] "UCSCKG"       "UNIPROT"     

We should see ENSEMBL, which is the type of key we are going to use in this case. If we are unsure what values are acceptable for the key, we can check what keys are valid with keys

keys(org.Hs.eg.db, keytype="ENSEMBL")[1:10]

 [1] "ENSG00000121410" "ENSG00000175899" "ENSG00000291190"
 [4] "ENSG00000171428" "ENSG00000156006" "ENSG00000196136"
 [7] "ENSG00000114771" "ENSG00000127837" "ENSG00000129673"
[10] "ENSG00000090861"

For the top gene in our analysis the call to the function would be:-

select(org.Hs.eg.db, keys="ENSG00000158258",
       keytype = "ENSEMBL",columns=c("SYMBOL","GENENAME")
)

Unfortunately, the authors of dplyr and AnnotationDbi have both decided to use the name select in their packages. To avoid confusion and errors, the following code is sometimes used:-

AnnotationDbi::select(org.Hs.eg.db, keys="ENSG00000158258",keytype = "ENSEMBL",columns=c("SYMBOL","GENENAME"))

To annotate our results, we definitely want gene symbols and perhaps the full gene name. Let’s build up our annotation information into a new data frame using the select function.

anno <- AnnotationDbi::select(org.Hs.eg.db,keys=rownames(dds),
              columns=c("SYMBOL","GENENAME"),
              keytype="ENSEMBL")
# Have a look at the annotation
head(anno)

However, we have a problem that the resulting data frame has more rows than our results table. This is due to the one-to-many relationships that often occur when mapping between various identifiers.

dim(anno)

[1] 58859     3

dim(dds)

[1] 57914     9

Such duplicated entries can be identified using the duplicated function. Fortunately, there are not too many so hopefully we won’t lose too much information if we discard the entries that are duplicated. The first occurrence of the duplicated ID will still be included in the table.

anno <- AnnotationDbi::select(org.Hs.eg.db,keys=rownames(dds),
            columns=c("ENSEMBL","SYMBOL","GENENAME","ENTREZID"),
            keytype="ENSEMBL") %>% 
filter(!duplicated(ENSEMBL))
dim(anno)

[1] 57914     4

We can bind in the annotation information to the results data frame.

results_annotated <- results(de_treated,tidy=TRUE) %>% 
  left_join(anno, by=c("row"="ENSEMBL"))

head(results_annotated)

We can save the results table using the write.csv function, which writes the results out to a csv file that you can open in excel.

write.csv(results_annotated,file="de_analysis/treatment_Y_vs_N_DESeq_annotated.csv",row.names=FALSE)
saveRDS(results_annotated, file="Robjects/treatment_Y_vs_N_DESeq_annotated.rds")

Exercise

• Join the annotation table to your results from the DESeq analysis of TGF vs CTR. Save the resulting data frame as a csv file. e.g. Robjects/results_TGF_vs_CTR_annotated.csv
• The publication gives examples of COL1A1, COL1A2 and COL3A1 as genes that are up-regulated in TGF-treated samples vs controls (Figure 6C). Use your data to verify this by
  • 1. extracting their p-values
  • 2. plotting the counts for these genes

Exporting normalized counts

The DESeq workflow applies median of ratios normalization that accounts for differences in sequencing depth between samples. The user does not usually need to run this step. However, if you want a matrix of counts for some application outside of Bioconductor the values can be extracted from the dds object.

dds <- estimateSizeFactors(dds) 
countMatrix <-counts(dds, normalized=TRUE) 
head(countMatrix)

                 1_CTR_BC_2  2_TGF_BC_4  3_IR_BC_5  4_CTR_BC_6
ENSG00000000003 1494.338259 1288.200211 1509.16643 1510.967861
ENSG00000000005    0.000000    0.000000    0.00000    0.000000
ENSG00000000419 1595.714814 1513.050580 2052.08122 1709.346140
ENSG00000000457  654.585928  536.353587  610.56771  680.178955
ENSG00000000460  247.768049  257.534062   82.94258  238.025591
ENSG00000000938    7.914116    3.633936   12.66761    4.767613
                 5_TGF_BC_7  6_IR_BC_12 7_CTR_BC_13
ENSG00000000003 1447.006807 1460.331190  1757.32486
ENSG00000000005    0.000000    0.000000     0.00000
ENSG00000000419 1882.326459 2161.853195  1943.82932
ENSG00000000457  672.369676  657.377968   747.96824
ENSG00000000460  300.401172   84.589046   326.58254
ENSG00000000938    2.013035    2.009128    25.63447
                8_TGF_BC_14  9_IR_BC_15
ENSG00000000003 1315.399310 1242.657289
ENSG00000000005    0.000000    1.022172
ENSG00000000419 1342.397617 2126.376893
ENSG00000000457  599.670748  652.546430
ENSG00000000460  290.916429  160.832221
ENSG00000000938    3.196483    6.971344

write.csv(countMatrix,file="normalized_counts.csv")

Full DESeq workflow

The median of ratios normalisation method is employed in DESeq2 to account for sequencing depth and RNA composition. Let’s go through a short worked example (courtesy of https://hbctraining.github.io/DGE_workshop/lessons/02_DGE_count_normalization.html) to explain the process.

## create a small example matrix of "counts"
test_data <- matrix(c(1489,22,793,76,521,906,13,410,42,1196),nrow=5)
rownames(test_data) <- c("EF2A","ABCD1","MEFV","BAG1","MOV10")
colnames(test_data) <- c("SampleA","SampleB")
test_data

      SampleA SampleB
EF2A     1489     906
ABCD1      22      13
MEFV      793     410
BAG1       76      42
MOV10     521    1196

Firstly, an “average” or reference sample is created that represents the counts on a typical sample in the dataset. The geometric mean is used rather than the arithmetic mean. In other words the individual counts are multiplied rather than summed and the measure should be more robust to outliers.

psuedo_ref <- sqrt(rowProds(test_data))
psuedo_ref

      EF2A      ABCD1       MEFV       BAG1      MOV10 
1161.47923   16.91153  570.20172   56.49779  789.37697 

A ratios of sample to “psuedo reference” are then calculated for each gene. We are assuming that most genes are not changing dramatically, so this ratio should be somewhere around 1.

test_data/psuedo_ref

        SampleA   SampleB
EF2A  1.2819859 0.7800398
ABCD1 1.3008873 0.7687061
MEFV  1.3907359 0.7190438
BAG1  1.3451854 0.7433919
MOV10 0.6600142 1.5151189

DESeq2 defines size factors as being the median of these ratios for each sample (median is used so any outlier genes will not affect the normalisation).

norm_factors <- colMedians(test_data/psuedo_ref)
norm_factors

  SampleA   SampleB 
1.3008873 0.7687061 

Individual samples can then normalised by dividing the count for each gene by the corresponding normalization factor.

test_data[,1] / norm_factors[1]

      EF2A      ABCD1       MEFV       BAG1      MOV10 
1144.60340   16.91153  609.58395   58.42166  400.49589 

and for the second sample…

test_data[,2] / norm_factors[2]

      EF2A      ABCD1       MEFV       BAG1      MOV10 
1178.60387   16.91153  533.36378   54.63727 1555.86118 

The size factors for each sample in our dataset can be calculated using the estimateSizeFactorsForMatrix function.

sf <- estimateSizeFactorsForMatrix(assay(dds))
sf

 1_CTR_BC_2  2_TGF_BC_4   3_IR_BC_5  4_CTR_BC_6  5_TGF_BC_7 
  1.0549852   1.1812947   0.8928435   1.0990482   1.0518069 
 6_IR_BC_12 7_CTR_BC_13 8_TGF_BC_14  9_IR_BC_15 
  0.8626010   0.9491286   1.0643407   0.9783918 

The estimation of these factors can also take gene-lengths into account, and this is implemented in the estimateSizeFactors function. Extra normalization factor data is added to the dds object.

dds <- estimateSizeFactors(dds)
dds

In preparation for differential expression DESeq2 also need a reliable estimate of the variability of each gene; which it calls dispersion.

dds <- estimateDispersions(dds)
dds

A statistical test can then be applied. As the data are count-based and not normally-distributed a t-test would not be appropriate. Most tests are based on a Poisson or negative-binomial distribution; negative binomial in the case of DESeq2. Although you might not be familiar with the negative binomial, the results should be in a familiar form with fold-changes and p-values for each gene.

dds <- nbinomWaldTest(dds)

It may seem like there is a lot to remember, but fortunately there is one convenient function that will apply the three steps (DESeq). The messages printed serve as reminders of the steps included.

Acknowledgements

Original Authors: Belinda Phipson, Anna Trigos, Matt Ritchie, Maria Doyle, Harriet Dashnow, Charity Law, Stephane Ballereau, Oscar Rueda, Ashley Sawle Based on the course RNAseq analysis in R delivered on May 11/12th 2016 and modified by Cancer Research Uk Cambridge Centre for the Functional Genomics Autumn School 2017

Resources and data files

This material has been created using the following resources:

• http://www.statsci.org/smyth/pubs/QLedgeRPreprint.pdf
• http://monashbioinformaticsplatform.github.io/RNAseq-DE-analysis-with-R/99-RNAseq_DE_analysis_with_R.html
  
• http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html
• https://bioconductor.github.io/BiocWorkshops/rna-seq-data-analysis-with-deseq2.html