RNA-seq Data Analysis

Load the packages

suppressPackageStartupMessages({
    library(here)
    library(dplyr)
    library(ggplot2)
    library(pheatmap)
    library(RColorBrewer)
    library(EnsDb.Hsapiens.v86)
    library(DESeq2)
})

Load the dataset

The transcript quantification results have been produced by salmon. The tximport package can import and summarize the information on gene level.

library(tximport)
tx2gene <- transcripts(EnsDb.Hsapiens.v86,
        return.type = "data.frame", columns = c("tx_id", "gene_id"))
salmon_dirs <- list.files(here("result/salmon_quant"), full.names = TRUE)
salmon_files <- sapply(salmon_dirs, function(x) file.path(x, "quant.sf"))
names(salmon_files) <- basename(salmon_dirs)
txi <- tximport(salmon_files, type = "salmon", tx2gene = tx2gene, ignoreTxVersion = TRUE)

## reading in files with read_tsv

## 1 2 3 4 5 6 7 
## transcripts missing from tx2gene: 7499
## summarizing abundance
## summarizing counts
## summarizing length

names(txi)

## [1] "abundance"           "counts"              "length"             
## [4] "countsFromAbundance"

We will use DESeq2 for further analysis. Convert the data to DESeqDataSet class.

samples <- data.frame(condition = c(rep.int("control", 3), rep.int("treat", 4)),
                      sample_name = names(salmon_files))
rownames(samples) <- names(salmon_files)
samples

##              condition  sample_name
## control_rep1   control control_rep1
## control_rep2   control control_rep2
## control_rep3   control control_rep3
## treat_rep1       treat   treat_rep1
## treat_rep2       treat   treat_rep2
## treat_rep3       treat   treat_rep3
## treat_rep4       treat   treat_rep4

dds <- DESeqDataSetFromTximport(txi, colData = samples, design = ~ condition)

## using counts and average transcript lengths from tximport

dds

## class: DESeqDataSet 
## dim: 37361 7 
## metadata(1): version
## assays(2): counts avgTxLength
## rownames(37361): ENSG00000000003 ENSG00000000005 ...
##   ENSG00000283698 ENSG00000283700
## rowData names(0):
## colnames(7): control_rep1 control_rep2 ... treat_rep3 treat_rep4
## colData names(2): condition sample_name

Remove the empty rows.

nrow(dds)

## [1] 37361

dds <- dds[ rowSums(counts(dds)) > 0, ]
nrow(dds)

## [1] 25001

Regularized-logarithm transformation (rlog)

We will use the regularized log transformation for PCA plot and SD-mean plot. It will reduce the variation on lowly expressed genes.

rld <- rlog(dds, blind = FALSE)

## using 'avgTxLength' from assays(dds), correcting for library size

head(assay(rld), 3)

##                 control_rep1 control_rep2 control_rep3 treat_rep1
## ENSG00000000003    11.300466    11.231395    11.104771  12.783082
## ENSG00000000419    11.171577    11.132015    11.321105  11.167904
## ENSG00000000457     8.207076     8.139017     7.899893   8.840509
##                 treat_rep2 treat_rep3 treat_rep4
## ENSG00000000003  12.758421  12.810201  12.789334
## ENSG00000000419  10.826120  10.920410  11.130422
## ENSG00000000457   8.943653   8.884305   8.846911

Heatmap showing sample distance (rlog transformed)

The samples can be clearly clustered by the control/treat condition.

sampleDists <- dist(t(assay(rld)))
sampleDistMatrix <- as.matrix( sampleDists )
colnames(sampleDistMatrix) <- NULL
colors <- colorRampPalette( rev(brewer.pal(9, "Blues")) )(255)
pheatmap(sampleDistMatrix,
         clustering_distance_rows = sampleDists,
         clustering_distance_cols = sampleDists,
         col = colors)

PCA plot (rlog transformed)

Previously we have conducted PCA analysis on samples. Here do PCA plot on genes.

plotPCA(rld)

Differential expression analysis

Run the differential expression pipeline.

dds <- DESeq(dds)

## estimating size factors

## using 'avgTxLength' from assays(dds), correcting for library size

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

Then build the results table, we use an adjusted p-value cutoff by 0.05.

res <- results(dds, alpha = 0.05)
head(res[order(res$pvalue), ], n = 3)

## log2 fold change (MLE): condition treat vs control 
## Wald test p-value: condition treat vs control 
## DataFrame with 3 rows and 6 columns
##                         baseMean    log2FoldChange             lfcSE
##                        <numeric>         <numeric>         <numeric>
## ENSG00000241343 10887.8871376054 -4.06463847625495 0.108594662387474
## ENSG00000185950 7290.71556646402  3.71238174924945 0.103648480871635
## ENSG00000183054 3308.40662655646  5.76746803913769 0.166493932393504
##                              stat                pvalue
##                         <numeric>             <numeric>
## ENSG00000241343 -37.4294499093518 1.29721089573246e-306
## ENSG00000185950  35.8170396520052 5.99696386772001e-281
## ENSG00000183054  34.6407100620727 6.16464222513363e-263
##                                  padj
##                             <numeric>
## ENSG00000241343  2.4815644435362e-302
## ENSG00000185950 5.73609593947419e-277
## ENSG00000183054 3.93098685889354e-259

mcols(res, use.names = TRUE)

## DataFrame with 6 rows and 2 columns
##                        type
##                 <character>
## baseMean       intermediate
## log2FoldChange      results
## lfcSE               results
## stat                results
## pvalue              results
## padj                results
##                                                       description
##                                                       <character>
## baseMean                mean of normalized counts for all samples
## log2FoldChange log2 fold change (MLE): condition treat vs control
## lfcSE                  standard error: condition treat vs control
## stat                   Wald statistic: condition treat vs control
## pvalue              Wald test p-value: condition treat vs control
## padj                                         BH adjusted p-values

summary(res)

## 
## out of 25001 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up)       : 5536, 22%
## LFC < 0 (down)     : 5371, 21%
## outliers [1]       : 54, 0.22%
## low counts [2]     : 5817, 23%
## (mean count < 1)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results

MA plot

Points are be colored red if the adjusted p value is less than 0.05.

plotMA(res)

Distribution of p values

hist(res$pvalue[res$baseMean > 1], breaks = 0:20/20,
     col = "grey50", border = "white")

Annotating the gene names

library(AnnotationDbi)
library(org.Hs.eg.db)
res$symbol <- mapIds(org.Hs.eg.db, keys=row.names(res),
    column="SYMBOL", keytype="ENSEMBL", multiVals="first")

## 'select()' returned 1:many mapping between keys and columns

res <- res[, c(ncol(res), 1:(ncol(res)-1))]
head(res[order(res$pvalue),], n = 3)

## log2 fold change (MLE): condition treat vs control 
## Wald test p-value: condition treat vs control 
## DataFrame with 3 rows and 7 columns
##                      symbol         baseMean    log2FoldChange
##                 <character>        <numeric>         <numeric>
## ENSG00000241343      RPL36A 10887.8871376054 -4.06463847625495
## ENSG00000185950        IRS2 7290.71556646402  3.71238174924945
## ENSG00000183054       RGPD6 3308.40662655646  5.76746803913769
##                             lfcSE              stat                pvalue
##                         <numeric>         <numeric>             <numeric>
## ENSG00000241343 0.108594662387474 -37.4294499093518 1.29721089573246e-306
## ENSG00000185950 0.103648480871635  35.8170396520052 5.99696386772001e-281
## ENSG00000183054 0.166493932393504  34.6407100620727 6.16464222513363e-263
##                                  padj
##                             <numeric>
## ENSG00000241343  2.4815644435362e-302
## ENSG00000185950 5.73609593947419e-277
## ENSG00000183054 3.93098685889354e-259

Generate a table of differentially expressed genes

We will select the genes that have adjusted p value < 0.05 and order them by their fold change.

Also, since there are a lot of genes differentially expressed, we may be interested in ones that have larger fold change. We will select the ones that have abs(log2foldchange) > 2.

selected <- res[!is.na(res$padj), ]
selected <- selected[selected$padj < 0.05, ]

We selected 10907 out of 25001 genes in the step above.

The following is the distribution of log2foldchange in genes with adjusted p value < 0.05. The yellow bars are ones with fold change < 2, which will be filtered.

hist(sort(selected$log2FoldChange), breaks = 30,
     col = c(rep("white", 14), rep("yellow", 2), rep("white", 14)), xlim = c(-15, 15))

selected <- selected[abs(selected$log2FoldChange) > 2, ]
selected <- selected[order(abs(selected$log2FoldChange), decreasing = TRUE), ]

After filtering by fold change, we finally selected 1885 genes.

Here we generate a table showing the top 100 rows.

selected <- as.data.frame(selected)
selected$log10pvalue <- log10(selected$pvalue)
selected$log10padj <- log10(selected$padj)
selected <- selected[ , c("symbol", "baseMean", "log2FoldChange", "pvalue",
                         "log10pvalue", "padj", "log10padj")]

knitr::kable(selected[1:100,])

Gene clustering

library(genefilter)

## 
## Attaching package: 'genefilter'

## The following objects are masked from 'package:matrixStats':
## 
##     rowSds, rowVars

mat <- assay(rld)[rownames(selected), ]
anno <- as.data.frame(colData(rld)[, c("condition", "sample_name")])
pheatmap(mat, annotation_col = anno, show_rownames = FALSE,
         main = "Heatmap of rlog-transformed gene expression levels")

mat <- mat - rowMeans(mat)
pheatmap(mat, annotation_col = anno, show_rownames = FALSE,
         main = "Heatmap of relative rlog-transformed values across samples")

library(genefilter)
topVarGenes <- head(order(rowVars(assay(rld)), decreasing = TRUE), 30)
mat <- assay(rld)[topVarGenes, ]
mat <- mat[rownames(mat) %in% rownames(selected), ]
mat <- mat - rowMeans(mat)
anno <- as.data.frame(colData(rld)[, c("condition", "sample_name")])
pheatmap(mat, annotation_col = anno,
         main = "Heatmap of relative rlog-transformed values across samples")

PCA plot of genes

We will use rlog-transformed values. The following is a PCA plot of all selected genes.

library(ggfortify)
autoplot(prcomp(assay(rld)[rownames(assay(rld)) %in% rownames(selected), ]), size = 1, alpha = 0.8)

It is clear that the genes can be separated by two parts by the second dimension. If we plot all the genes like the following, it will be clear for us that the second dimension aligns with the foldchange – since we have filtered the genes with low foldchange. Thus the two clusters are up-regulated genes and down-regulated genes.

autoplot(prcomp(assay(rld)[rownames(assay(rld)), ]), size = 1, alpha = 0.8)

topVarGenes <- head(order(rowVars(assay(rld)), decreasing = TRUE), 50)
tmp <- assay(rld)[topVarGenes, ]
tmp <- tmp[rownames(tmp) %in% rownames(selected), ]
autoplot(prcomp(tmp))

Pathway analysis

The package ReactomePA offers the possibility to test enrichment of specific pathways using the free, open-source, curated and peer reviewed pathway Reactome pathway database.

library(ReactomePA)
entrez_ids <- mapIds(org.Hs.eg.db, keys = rownames(res),
                     column = "ENTREZID", keytype = "ENSEMBL", multiVals = "first")

## 'select()' returned 1:many mapping between keys and columns

entrez_ids <- entrez_ids[rownames(selected)]
entrez_ids <- entrez_ids[!is.na(entrez_ids)]

reactome_enrich <- enrichPathway(gene = entrez_ids, organism = "human", readable = TRUE)

Then we visualize the enriched terms.

barplot(reactome_enrich, showCategory = 25)

dotplot(reactome_enrich, showCategory = 25)

ReactomePA::emapplot(reactome_enrich, showCategory = Inf)

up_regulated <- selected[selected$log2FoldChange > 0, ]
up_regulated <- up_regulated[rownames(up_regulated) %in% names(entrez_ids), ]
tryCatch(
    ReactomePA::emapplot(
        enrichPathway(gene = entrez_ids[rownames(up_regulated)], organism = "human"),
        showCategory = Inf),
    error = function(e) print(e$message)
)

## [1] "no enriched term found..."

down_regulated <- selected[selected$log2FoldChange < 0, ]
down_regulated <- down_regulated[rownames(down_regulated) %in% names(entrez_ids), ]
ReactomePA::emapplot(
    enrichPathway(gene = entrez_ids[rownames(down_regulated)], organism = "human"),
    showCategory = Inf)