RNAseq in Bioconductor exercises (part 3)

In todays session we will work with some of the RNAseq data of adult mouse tissues from Bing Ren’s lab, Liver and Heart.

• More information on liver data can be found here
• More information on heart data can be found here
• More information on Kidney data can be found here

Preprocessed counts for Tissue data for this practical session can be found in the data directory. data/tissueCounts.RData

Exercises

1. PCA

Load in the DESeq2 object for Tissues experiment, produce an rlog transformation and plot the PCA for these genes.

library(DESeq2)
library(ggplot2)
load("data/gC_TissueFull.RData")
dds <- DESeqDataSet(geneCounts, design = ~Tissue)
dds <- DESeq(dds)

rlogTissue <- rlog(dds)
myPlot <- DESeq2::plotPCA(rlogTissue, intgroup = "Tissue")
myPlot+theme_bw()

2. LRT test

Filter to genes significant (padj < 0..01) from our LRT test between all groups and produce a heatmap.

dds2 <- DESeq(dds,test="LRT",reduced = ~1)
AllChanges <- results(dds2)
rlogMatrix <- assay(rlogTissue)
sigChanges <- rownames(AllChanges)[AllChanges$padj < 0.01 & !is.na(AllChanges$padj)]
sigMat <- rlogMatrix[rownames(rlogMatrix) %in% sigChanges,]
annoDF <- as.data.frame(colData(rlogTissue)[,1,drop=FALSE])
library(pheatmap)
pheatmap(sigMat,
         scale="row",
         show_rownames = FALSE,
         annotation_col = annoDF)

3. PCA loadings

Extract the influence of genes to PC2 and produce a heatmap as above but now ranked by PC2 rotation/loadings.

pcRes <- prcomp(t(rlogMatrix))
PC2_rnk <- sort(pcRes$rotation[,2],decreasing = TRUE)
PC2_mat <- sigMat[match(names(PC2_rnk),rownames(sigMat),nomatch = 0),]
pheatmap(PC2_mat,
         scale="row",
         cluster_rows=FALSE,
         show_rownames = FALSE,annotation_col = annoDF
         )

4. PCA of a GO term

Produce a PCA plot of the GO “heart development” genes (“GO:0007507”). You will need to extract all genes within this group from the org.Mm.eg.db package.

library(org.Mm.eg.db)
HeartDevelopment <- AnnotationDbi::select(org.Mm.eg.db, keytype = "GOALL",
                              keys = "GO:0007507", columns = "ENTREZID")
HeartDevelopment_Entrez <- unique(HeartDevelopment$ENTREZID)
DESeq2::plotPCA(rlogTissue[rownames(rlogTissue) %in% HeartDevelopment_Entrez],intgroup = "Tissue") + theme_bw()

5. Heatmap of a GO term

Create a Heatmap of all genes and include row annotation showing members of heart development, liver development and kidney development GO set.

library(org.Mm.eg.db)
HeartDevelopment <- AnnotationDbi::select(org.Mm.eg.db,keytype = "GOALL",
                              keys = "GO:0007507",columns = "ENTREZID")
HeartDevelopment_Entrez <- unique(HeartDevelopment$ENTREZID)
library(org.Mm.eg.db)
LiverDevelopment <- AnnotationDbi::select(org.Mm.eg.db,keytype = "GOALL",
                              keys = "GO:0001889",columns = "ENTREZID")
LiverDevelopment_Entrez <- unique(LiverDevelopment$ENTREZID)
library(org.Mm.eg.db)
KidneyDevelopment <- AnnotationDbi::select(org.Mm.eg.db,keytype = "GOALL",
                              keys = "GO:0001822",columns = "ENTREZID")
KidneyDevelopment_Entrez <- unique(KidneyDevelopment$ENTREZID)
annoRow <- data.frame(HeartDev=factor(rownames(sigMat) %in% HeartDevelopment_Entrez),
           LiverDev=factor(rownames(sigMat) %in% LiverDevelopment_Entrez),
           KidneyDev=factor(rownames(sigMat) %in% KidneyDevelopment_Entrez))
rownames(annoRow) <- rownames(sigMat)

ann_colors = list(
    HeartDev = c("FALSE"="white","TRUE"="green"),
    LiverDev = c("FALSE"="white","TRUE"="red"),
    KidneyDev = c("FALSE"="white","TRUE"="blue")
)

pheatmap(sigMat,
         scale="row",
         show_rownames = FALSE,
         annotation_col = annoDF,
         annotation_row = annoRow,
         annotation_colors = ann_colors)

6. Clustering

Take the LRT filtered rlog counts and cluster using pheatmap into 3 clusters,

set.seed(153)
k <-   pheatmap(sigMat,
           scale="row",kmeans_k = 3)

7. Heatmap ordered by cluster

names(k$kmeans)

## [1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss"
## [6] "betweenss"    "size"         "iter"         "ifault"

clusterDF <- as.data.frame(factor(k$kmeans$cluster))
colnames(clusterDF) <- "Cluster"
clusterDF[1:10,,drop=FALSE]

##           Cluster
## 20671           1
## 27395           1
## 18777           2
## 21399           3
## 108664          2
## 319263          3
## 76187           2
## 70675           3
## 73824           3
## 100039596       2

OrderByCluster <- sigMat[order(clusterDF$Cluster),]

pheatmap(OrderByCluster,
           scale="row",annotation_row = clusterDF,
           show_rownames = FALSE,cluster_rows = FALSE)