MAIN BIOINFORMATICS ANALYZES for RNA-seq DATA
Yusuf Çağlar ODABAŞI
Data
Rett syndrome is an X-linked neurodevelopmental disorder caused by mutation in the methyl-CpG-binding protein 2 gene in the majority of cases.RNA sequencing dataset of postmortem brain tissue samples from four females clinically diagnosed with Rett syndrome and four age-matched female donors. The dataset contains transcriptomes from two brain regions, temporal and cingulate cortex (GSE128380). We are going to use count matrix, tpm_matrix and phenotype list. (RTT:mutated, WT:healthy control)
Citation: Aldinger KA, Timms AE, MacDonald JW, McNamara HK et al. Transcriptome data of temporal and cingulate cortex in the Rett syndrome brain. Sci Data 2020 Jun 19;7(1):192. PMID: 32561870
library(DESeq2)
library(pheatmap)
library(EnhancedVolcano)
library(wrMisc)
library(ggfortify)
library(ggplot2)
library(enrichplot)
library(clusterProfiler)
organism = "org.Hs.eg.db"
library(organism, character.only = TRUE)
library(stringr)
require(DOSE)
library(pathview)load("Count.Rdata")
head(Count)## GSE1283856.PM.WT. GSE1283857.PM.WT. GSE1283858.PM.WT. GSE1283859.PM.WT.
## A1BG 3 1 0 1
## A1CF 0 1 2 0
## A2M 3331 2768 2357 2145
## A2ML1 73 56 90 19
## A3GALT2 1 0 0 0
## A4GALT 31 20 22 16
## GSE1283860.PM.RTT. GSE1283861.PM.RTT. GSE1283862.PM.RTT.
## A1BG 0 0 6
## A1CF 1 2 1
## A2M 1910 3500 1341
## A2ML1 12 41 124
## A3GALT2 1 0 0
## A4GALT 65 115 75
## GSE1283863.PM.RTT. GSE1283864.PM.WT. GSE1283865.PM.WT.
## A1BG 2 0 5
## A1CF 0 4 1
## A2M 3856 3722 2144
## A2ML1 88 56 74
## A3GALT2 0 0 0
## A4GALT 57 23 20
## GSE1283866.PM.WT. GSE1283867.PM.WT. GSE1283868.PM.RTT.
## A1BG 14 2 0
## A1CF 0 2 4
## A2M 3811 3079 2812
## A2ML1 62 16 19
## A3GALT2 0 0 0
## A4GALT 26 26 92
## GSE1283869.PM.RTT. GSE1283870.PM.RTT. GSE1283871.PM.RTT.
## A1BG 4 0 5
## A1CF 2 1 0
## A2M 3426 685 5404
## A2ML1 40 107 75
## A3GALT2 2 0 0
## A4GALT 142 16 38load("TPM.Rdata")
head(TPM)## GSE1283856.PM.WT. GSE1283857.PM.WT. GSE1283858.PM.WT. GSE1283859.PM.WT.
## A1BG 0.12130266 0.07847415 0.00000000 0.06400858
## A1CF 0.00000000 0.01291718 0.02356503 0.00000000
## A2M 39.05161859 69.42452255 53.93546082 33.55941065
## A2ML1 1.02755223 1.59882204 2.46413979 0.37489325
## A3GALT2 0.06609256 0.00000000 0.00000000 0.00000000
## A4GALT 0.96055831 1.23477672 1.28496757 0.62050687
## GSE1283860.PM.RTT. GSE1283861.PM.RTT. GSE1283862.PM.RTT.
## A1BG 0.00000000 0.00000000 0.258178510
## A1CF 0.01464856 0.01236089 0.007478194
## A2M 54.29478765 42.00916693 19.471147233
## A2ML1 0.40143644 0.54438988 3.080066253
## A3GALT2 0.16022468 0.00000000 0.000000000
## A4GALT 4.29914513 3.38732599 2.879330660
## GSE1283863.PM.RTT. GSE1283864.PM.WT. GSE1283865.PM.WT.
## A1BG 0.06289642 0.00000000 0.31305443
## A1CF 0.00000000 0.02416994 0.01081742
## A2M 48.08659724 44.12291877 45.04060532
## A2ML1 1.83138746 0.81271410 2.47555093
## A3GALT2 0.00000000 0.00000000 0.00000000
## A4GALT 1.75078190 0.71363235 0.99983207
## GSE1283866.PM.WT. GSE1283867.PM.WT. GSE1283868.PM.RTT.
## A1BG 0.5046040 0.05688404 0.00000000
## A1CF 0.0000000 0.01178670 0.03401098
## A2M 46.2144643 35.24174167 46.44238194
## A2ML1 1.0646610 0.41895132 0.36875048
## A3GALT2 0.0000000 0.00000000 0.00000000
## A4GALT 0.8047624 0.78680405 3.71115648
## GSE1283869.PM.RTT. GSE1283870.PM.RTT. GSE1283871.PM.RTT.
## A1BG 0.14019499 0.00000000 0.1744825
## A1CF 0.01190009 0.00938817 0.0000000
## A2M 39.58997723 12.48651061 58.8325755
## A2ML1 0.57991689 3.04151350 0.9485134
## A3GALT2 0.13029447 0.00000000 0.0000000
## A4GALT 4.12463185 0.72310813 1.0586425phenotype <- read.csv("phenotype.csv", sep = " ")
head(phenotype)## Run Assay.Type phenotype tissue
## 1 GSE128380 RNA-Seq WT cingulate
## 2 GSE128380 RNA-Seq WT cingulate
## 3 GSE128380 RNA-Seq WT cingulate
## 4 GSE128380 RNA-Seq WT cingulate
## 5 GSE128380 RNA-Seq RTT cingulate
## 6 GSE128380 RNA-Seq RTT cingulatePreprocessing and QC
In order to reduce the noise in the data matrix and obtain unbiased expression differences between WT and RTT, average values of each row of the TPM matrix were calculated, and the rows with average value ≤1.0 were filtered. Then, log2 transformation was applied to the TPM matrix, and the distribution across samples was examined by using box plots. Coefficient of variation (CV) values of each gene were computed and 2000 genes with the highest CV values were elicited for the PCA. In addition to the PCA, we created a sample-to-sample Pearson correlation heatmap of our data to inspect similarities of each sample to others by using the top 3500 genes from the CV-sorted TPM matrix before DEG analysis.
# -------------taking mean of rows--------------------
row_mean <- rowMeans(TPM)
TPM <- as.data.frame(cbind(TPM, tpmMean = row_mean))
# ---eliminate genes that have meantpm values lower than 1----
filtered_TPM <- TPM[TPM$tpmMean>1,]
# -------------taking log2 transformation--------------
filtered_TPM <- log2(filtered_TPM[,-17] + 1)
# ---------------------plot boxplot-------------------------
par(mar=c(11,4,4,4))
boxplot(filtered_TPM, ylim =c(0,16), ylab="log2(tpm_values +1)", las=2)
# -------------------coefficient to variation, sort, get top 2000 genes--------------------
CV_values <- rowCVs(filtered_TPM)
CV <- cbind(filtered_TPM, cv_values = CV_values)
order <- order(CV[,17],decreasing=TRUE)
CV_SORT <- CV[order(CV[,17],decreasing=TRUE),]
CV_SORT_TOP2000 <- CV_SORT[c(1:2000),c(1:16)]
# ---------------------plot PCA for top 2000 genes-------------------------
T_CV_SORT_TOP2000 <- t(CV_SORT_TOP2000)
PCA <- prcomp(T_CV_SORT_TOP2000,scale. = TRUE)
pca1 <- autoplot(PCA, data = phenotype, colour ='phenotype')
pca2 <- autoplot(PCA, data = phenotype, colour ='tissue')
par(mfrow=c(1,2))
pca1
pca2
#-----------------Sample-to-Sample-Heatmap--------------------
CV_values <- rowCVs(Count)
Count <- cbind(Count, cv_values = CV_values)
Count_CV_SORT <- Count[order(Count[,17],decreasing=TRUE),]
Count_CV_TOP3500 <- Count_CV_SORT[c(1:3500),-17]
top3500 <-rownames(Count_CV_TOP3500)
cv_filtered_Count_35 <- Count[rownames(Count) %in% top3500, -17]
pheatmap::pheatmap(cor(cv_filtered_Count_35))
Differential Gene Expression (DGE) Analysis
We initially sought to determine individual gene expression differences between samples originated from patients RTT and WT Therefore, we first determined differentially expressed genes (DEG) between RTT samples and their normal counterparts within each sample type. The count matrix was analyzed with the R package “DESeq2”
#-------------------------------Deseq2 Object and Volcano Plot-----------------------------
phenotype$phenotype <- factor(phenotype$phenotype)
factor_col <- sapply(phenotype,is.factor)
factor_col## Run Assay.Type phenotype tissue
## FALSE FALSE TRUE FALSECount <- Count[,-17]
rownames(phenotype) <- colnames(Count)
dds <- DESeqDataSetFromMatrix(Count, phenotype,~phenotype)
dds$phenotype <- relevel(dds$phenotype, ref = "WT")
dds <- DESeq(dds)
resultsNames(dds)## [1] "Intercept" "phenotype_RTT_vs_WT"res <- results(dds, alpha=0.05)
res <- na.omit(res)
saveRDS(res, "practice_deseq2.rds")
up <- subset(res, res$log2FoldChange > 1 & res$padj < 0.05)
up <- up[order(-up$log2FoldChange),]
up_list <- head(rownames(up), n=10L)
down <- subset(res, res$log2FoldChange < -1 & res$padj < 0.05)
down <- down[order(down$log2FoldChange),]
down_list <- head(rownames(down), n=10L)
plot_top_20 <- c(up_list, down_list)
volc <- EnhancedVolcano(res,
lab=rownames(res),
x="log2FoldChange",
y="padj",
selectLab=plot_top_20,
drawConnectors=TRUE,
legendPosition = "none",
FCcutoff=1,
pCutoff=0.05
)
volc
Gene Ontology (GO) Analysis
To further investigate all differentially expressed genes, we conducted a Gene Ontology analysis. DEG lists were given to R package “Clusterprofiler” and gseGO() function.
up_genes <- up$log2FoldChange
names(up_genes) <- rownames(up)
down_genes <- down$log2FoldChange
names(down_genes) <- rownames(down)
down_genes <- sort(down_genes, decreasing = TRUE)
DE_genes <- c(up_genes,down_genes)
gse <- gseGO(geneList=DE_genes,
ont ="ALL",
keyType = "SYMBOL",
#nPerm = 10000,
minGSSize = 3,
maxGSSize = 800,
pvalueCutoff = 0.05,
verbose = TRUE,
OrgDb = org.Hs.eg.db,
pAdjustMethod = "none")
p1 <- dotplot(gse, showCategory=5, split=".sign", font.size = 7.5, label_format = 24) + facet_grid(.~.sign)
p1
cex.params = list(category_label = 0.1, gene_label = 0.4)
p2 <- cnetplot(gse, categorySize="pvalue", foldChange=DE_genes, showCategory = 10, circular=T, colorEdge=T, cex.params=cex.params)
p2
p3 <- ridgeplot(gse, showCategory = 10, label_format =28) + labs(x = "enrichment distribution")
p3
gse<-pairwise_termsim(gse)
p4 <- emapplot(gse, showCategory = 10)
p4
T_genes <-bitr(names(DE_genes), fromType = "SYMBOL", toType = "ENTREZID", OrgDb=organism)
EID_genes <- DE_genes
names(EID_genes) <- T_genes$ENTREZID
kegg = gseKEGG(geneList = EID_genes,
organism = "hsa",
minGSSize = 3,
maxGSSize = 800,
pvalueCutoff = 1,
pAdjustMethod = "none",
keyType = "ncbi-geneid")
sig_res <- kegg@result[kegg@result$p.adjust < 0.05,]
kegg_pathways_id <- head(sig_res$ID, n=3L)
kegg_pathways_id## [1] "hsa05150" "hsa04613" "hsa04657"fig <- pathview(gene.data=EID_genes, pathway.id=kegg_pathways_id, species = "hsa")
knitr::include_graphics(c("hsa05150.pathview.png", "hsa04270.pathview.png", "hsa04657.pathview.png"))
