RNA-seq实战(下)_数据可视化、GO/KEGG分析

2019年7月6日 1,184次阅读 来源: 谢俊飞

前言：
写这篇文章的目的是为了梳理一下学习思路，按部就班地仿生信菜鸟团和简书：Y大宽
教程大纲，做归纳整理，即便再次运行，仍然遇到代码出错或者软件使用等诸多问题，是以为之记。

MA plot

详细关于Bland-Altman图具体含义，请移步至

  1. R语言：Bland-Altman分析
  2. 连续变量的一致性评价
> plotMA(res,ylim=c(-2,2))
> topGene <- rownames(res)[which.min(res$padj)]
> with(res[topGene, ], {
+   points(baseMean, log2FoldChange, col="dodgerblue", cex=6, lwd=2)
+   text(baseMean, log2FoldChange, topGene, pos=2, col="dodgerblue")
+ })

结果如下：

《RNA-seq实战(下)_数据可视化、GO/KEGG分析》 Rplot01.png

经过lfcShrink 收缩log2 fold change

> res_order<-res[order(row.names(res)),]
> res = res_order
> res.shrink <- lfcShrink(dds, contrast = c("condition","treat","control"), res=res)
> plotMA(res.shrink, ylim = c(-5,5))
> topGene <- rownames(res)[which.min(res$padj)]
> with(res[topGene, ], {
+   points(baseMean, log2FoldChange, col="dodgerblue", cex=2, lwd=2)
+   text(baseMean, log2FoldChange, topGene, pos=2, col="dodgerblue")
+ })
> plotMA(res.shrink, ylim = c(-5,5))
> topGene <- rownames(res)[which.min(res$padj)]
> with(res[topGene, ], {
+   points(baseMean, log2FoldChange, col="dodgerblue", cex=2, lwd=2)
+   text(baseMean, log2FoldChange, topGene, pos=2, col="dodgerblue")
+ })

结果如下：

《RNA-seq实战(下)_数据可视化、GO/KEGG分析》 Rplot02.png

2. Plot counts

# 不画图，只显示数据
plotCounts(dds, gene=which.min(res$padj), intgroup="condition", returnData=TRUE)
#只画图，不显示数据
plotCounts(dds, gene="ENSG00000002586", intgroup="condition", returnData=FAULSE)

下面用ggplot2来画CD99的box图和point图

  • boxplot
# Plot it
> plotCounts(dds, gene="ENSG00000002586", intgroup="condition", returnData=TRUE) %>% 
+   ggplot(aes(condition, count)) + geom_boxplot(aes(fill=condition)) + scale_y_log10() + ggtitle("CD99")

《RNA-seq实战(下)_数据可视化、GO/KEGG分析》 Rplot03.png

  • point plot
> d <- plotCounts(dds, gene="ENSG00000002586", intgroup="condition", returnData=TRUE)
> ggplot(d, aes(x=condition, y=count)) + 
+   geom_point(aes(color= condition),size= 4, position=position_jitter(w=0.5,h=0)) + 
+   scale_y_log10(breaks=c(25,100,400))+ ggtitle("CD99")

《RNA-seq实战(下)_数据可视化、GO/KEGG分析》 Rplot04.png

3. PCA(principal components analysis)

> vsdata <- vst(dds, blind=FALSE)
> plotPCA(vsdata, intgroup="condition")

《RNA-seq实战(下)_数据可视化、GO/KEGG分析》 Rplot05.png

4. 热图：两部分

4.1 count matrix 热图

> library("pheatmap")
> select<-order(rowMeans(counts(dds, normalized = TRUE)),
+               decreasing = TRUE)[1:20]
> df <- as.data.frame(colData(dds)[,c("condition","sizeFactor")])
 # this gives log2(n + 1)
> ntd <- normTransform(dds)
> pheatmap(assay(ntd)[select,], cluster_rows=FALSE, show_rownames=FALSE,
+          cluster_cols=FALSE, annotation_col=df)

《RNA-seq实战(下)_数据可视化、GO/KEGG分析》 Rplot06.png

4.2 sample-to-sample distances热图

> sampleDists <- dist(t(assay(vsdata)))
> library("RColorBrewer")
> sampleDistMatrix <- as.matrix(sampleDists)
> rownames(sampleDistMatrix) <- paste(vsdata$condition, vsdata$type, sep="-")
> colnames(sampleDistMatrix) <- NULL
> colors <- colorRampPalette( rev(brewer.pal(9, "Blues")) )(255)
> pheatmap(sampleDistMatrix,
+          clustering_distance_rows=sampleDists,
+          clustering_distance_cols=sampleDists,
+          col=colors)

《RNA-seq实战(下)_数据可视化、GO/KEGG分析》 Rplot07.png

5. 基因集(gene set)富集分析

# enrichment analysis using clusterprofiler package created by yuguangchuang
> library(clusterProfiler)
> library(DOSE)
> library(org.Hs.eg.db)
> sig.gene <- read.csv(file="readcount_all.csv")# 我用DEG_treat_vs_control.csv后续分析没有富集到，所以换用这个
> head(sig.gene)
> gene <- sig.gene[,1] #注意是ID所在的列
> head(gene)
> gene.df <- bitr(gene, fromType = "ENSEMBL", 
              toType = c("SYMBOL","ENTREZID"),
              OrgDb = org.Hs.eg.db)
> head(gene.df)

5.1 GO富集分析：

#Go classification
#Go enrichment
> library("stringr")  #后面出图的参数str_wrap需要这个包
ego_cc <- enrichGO(gene       = gene.df$ENSEMBL,
                 OrgDb      = org.Hs.eg.db,
                 keyType    = 'ENSEMBL',
                 ont        = "CC",
                 pAdjustMethod = "BH",
                 pvalueCutoff = 0.01,
                 qvalueCutoff = 0.05)
ego_bp <- enrichGO(gene       = gene.df$ENSEMBL,
                 OrgDb      = org.Hs.eg.db,
                 keyType    = 'ENSEMBL',
                 ont        = "BP",
                 pAdjustMethod = "BH",
                 pvalueCutoff = 0.01,
                 qvalueCutoff = 0.05)
> barplot(ego_cc,showCategory = 18,title="The GO_CC enrichment analysis of all DEGs ")+ 
+   scale_size(range=c(2, 12))+
+   scale_x_discrete(labels=function(ego_cc)  str_wrap(ego_cc,width = 25))
> barplot(ego_bp,showCategory = 18,title="The GO_BP enrichment analysis of all DEGs ")+ 
+   scale_size(range=c(2, 12))+
+   scale_x_discrete(labels=function(ego_bp) str_wrap(ego_bp,width = 25))

《RNA-seq实战(下)_数据可视化、GO/KEGG分析》 Rplot10.png

《RNA-seq实战(下)_数据可视化、GO/KEGG分析》 Rplot11.png

5.2 KEGG enrichment

> install.packages("stringr")
> library(stringr)
> kk<-enrichKEGG(gene      =gene.df$ENTREZID,
               organism = 'hsa',
               pvalueCutoff = 0.05)
> kk[1:30]
> barplot(kk,showCategory = 25, title="The KEGG enrichment analysis of all DEGs")+
  scale_size(range=c(2, 12))+
  scale_x_discrete(labels=function(kk) str_wrap(kk,width = 25))
> dotplot(kk,showCategory = 25, title="The KEGG enrichment analysis of all DEGs")+
  scale_size(range=c(2, 12))+
  scale_x_discrete(labels=function(kk) str_wrap(kk,width = 25))

《RNA-seq实战(下)_数据可视化、GO/KEGG分析》 Rplot12.png
《RNA-seq实战(下)_数据可视化、GO/KEGG分析》 Rplot13.png

6. Gene Set Enrichment Analysis（GSEA）

我用的数据库这部分没有找到差异的，所以只是罗列操作步骤，待以后有数据再做处理

# 获取按照log2FC大小来排序的基因列表
> sig.gene <- read.csv(file="DEG_treat_vs_control.csv")
> genelist <- sig.gene$log2FoldChange
> names(genelist) <- sig.gene[,1]
> genelist <- sort(genelist, decreasing = TRUE)
# GSEA分析
> gsemf <- gseGO(genelist,
               OrgDb = org.Hs.eg.db,
               keyType = "ENSEMBL",
               ont="BP"
)
# 查看信息
head(gsemf)
# 画出GSEA图
gseaplot(gsemf, geneSetID="GO:0001819")
    原文作者：谢俊飞
    原文地址: https://www.jianshu.com/p/f73de8022cf9
    本文转自网络文章，转载此文章仅为分享知识，如有侵权，请联系博主进行删除。
点赞