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.

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(42)
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)