TCGA_RNA-Seq下游分析
Hhuang
2019/11/20
前言
先期通过TCGA下载得到HT－Seq的测序数据，主要是COUNT的数据，经过整理变换后得到Rdata数据文件（已整合）。此文为接续下来的分析步骤。
一、加载数据
load(file="expr_df.Rdata") ## 此处根据实际写入数据库路径这个数据就是测序后的counts矩阵，可以通过测序公司取得相关数据。
可以观察到如下数据格式
Rdata由一列组成，第一列是ensemble的ID，剩下的6列是3Vs3的样本，也就是通常转录组需要的6个样本。(要确定一下这个数据的类型是data.frame，如果你喜欢用data.table::fread读入文件，那么读进来的是tibble，需要使用data.frame转换一下，要不然Deseq2会报错)
构建metadata文件
因为我们需要用Deseq2这个包来分析差异，那么就要先制作一个metadata文件指示分组信息,他是一个数据框，由两列组成。
!
metadata <- data.frame(sample_id = colnames(expr_df)[-1])
sample <- rep(c("BRCA","CNT"),each=10)
metadata$sample <- relevel(factor(sample),"BRCA")构建dds对象
这一步由DESeqDataSetFromMatrix这个函数来完成，他需要输入我们的表达矩阵，制作好的metadata，还要制定分组的列，在这里是sample，最后一个tidy的意思是，我们第一列是基因ID，需要自动处理。
## results='hide' 用来隐藏运行结果。
## eval=FALSE 用来禁止代码运行
## message=FALSE warning=FALSE 关闭信息提示与警告信息显示
## BiocManager::install("DESeq2)
## 需要提前安装好 DESeq2, GenomeicRanges, GenomeInfoDb
library(ggplot2)
library(DESeq2)
dds <-DESeqDataSetFromMatrix(countData=expr_df,
colData=metadata,
design=~sample,
tidy=TRUE)
## 第一列有名称，所以tidy=TRUE看一下有多少行
nrow(dds)## [1] 60483查看一下有行名么？有的，如果没有就不要往下做了，应该是你的数据不是data.frame格式
rownames(dds)数据过滤
别看有这么多行，不是每一个就要都要表达的，如果一行基因在所有样本中的counts数小于等于1，我们就把他删掉，实际上，不做这一步，对差异分析的结果没有影响，可能会对GSEA的结果有影响。
···{r filter,echo=TRUE,error=FALSE,message=FALSE,warning=FALSE} dds <- dds[rowSums(counts(dds))>1,] nrow(dds) ##再来看他的行数,是不是变少了 ···
样本聚类
注意一定要分清实验组与对照组的分类标签。 用vst来标化数据，实际上还有rlog方法，或者就是log2的方法，官网推荐30个样本用rlog，大于30个样本用vst，速度更快，这里我们不要计较那么多了，就用vst，因为真实的TCGA数据，样本往往大于30个。
##vsd <- vst(dds, blind = FALSE)
vsd <- rlog(dds,blind=FALSE,fitType = "parametric") ## 此处也可以使用 vst或 log2
## 再用dist这个函数计算样本间的距离
sampleDists <- dist(t(assay(vsd)))
## 用hclust来进行层次聚类
hc <- hclust(sampleDists, method = "ward.D2") ## 聚类分析，也可以用 single,最短距离法
## 最后作图
plot(hc, hang = -1)
re1<-rect.hclust(hc, k=2, border="red") ## 分成2类，加边框线
***
横向图来表示
library(factoextra)
res <- hcut(sampleDists, k = 2, stand = TRUE,hc_method="ward.D2",hc_func="hclust",graph=TRUE)
## Visualize
fviz_dend(res,
# 加边框
rect = TRUE,
# 边框颜色
rect_border="cluster",
# 边框线条类型
rect_lty=2,
# 边框线条粗细
lwd=1.2,
# 边框填充
rect_fill = T,
# 字体大小
cex = 1,
# 字体颜色
color_labels_by_k=T,
# 平行放置
horiz=T)
fviz_cluster(res)
二、Deseq2 计算
主程序是Deseq这个函数，里面顺序执行了一系列函数，每一步都可以单独运行。这一步，只有6个样本基本上就是10s以内，如果是1000个样本，小电脑跑不过去，跑过去也需要5个小时以上，很耗时间。
dds <- DESeq(dds)
!
得到dds之后，我们可以通过counts这个函数得到能作图的标注化后的counts数据，他矫正了样本间测序的深度，使得样本间可以直接比较。
normalized_counts <- as.data.frame(counts(dds, normalized=TRUE))
!
选取其中一个基因来作图
Deseq2内置了一个画图函数，可以方便地制定基因作图
plotCounts(dds, gene = "ENSG00000223572", intgroup=c("sample"))
这个功能本质上没有什么用，但是，可以提前确定实验的质量。比如，你的两组是敲减某个基因以及对照组，通过制定那个基因作图，就可以看出实验有没有成功，如果这个基因没有任何改变，也可以不用往下做了。回去重新做实验送样本吧。
继续美化一下，有个内置的参数returnData可以返回作图数据。
一旦返回数据，我们就可以用ggplot2自己简单画一下。
plotdata <- plotCounts(dds, gene = "ENSG00000172183", intgroup=c("sample"),returnData = T)
## install.packages("ggthemes") #安装配色主题包
library(ggthemes)
##普通主题绘图命令
#ggplot(plotdata,aes(x=sample,y=count,col=sample))+
#geom_point()+
#theme_bw()
##改造后可叠加图形的命令
ggplot(plotdata,aes(x=sample,y=count,col=sample))+
geom_point(size=2,alpha=0.8,shape=21)+geom_point(aes(size=count))+
#geom_violin(fill="lightblue") + geom_jitter(shape = 18)+geom_rug(color="blue")+
geom_boxplot(width=2)+
theme_stata()+theme(legend.position = "right")
现在看起来平淡无奇，如果样本多，可以画出点图配boxplot，如果是配对样本，那么还可以画出配对的图。
有一点要强调一下，vst，以及log2的标准化，跟标准化的counts不是一个概念。前者是为了以后的聚类，热土，PCA分析，比如，我们计算样本间的距离就是用的vst 标化的数据，而标准化的counts是为了差异作图，你看纵坐标就会发现，他的数值一般很大。
三、LogFC的矫正
这一步，对于依赖logFC变化值的分析，很重要，比如GSEA分析。我们画一个MAplot图，看图说话。
contrast <- c("sample", "BRCA", "CNT")
dd1 <- results(dds, contrast=contrast, alpha = 0.05)
plotMA(dd1, ylim=c(-2,2))
MA-plot上的点是每一个基因。横坐标是标准化后的counts平均值，纵坐标是log后的变化值。红色的点是p.adjust的值小于0.05的基因，由counts函数中的alpha 参数指定。 我们发现在左侧，有很多counts很小的基因，发生了很大的变化，但是没有明显意义。类似于从(1,3,9)变成了（20,12,3）他们的counts很小，波动性很大，对logFC产生了很大的影响。 GSEA分析中，排序就是按照logFC来进行的。按照这个结果往下做，GSEA那里，富集不到任何条目。
那就需要矫正。用的函数是lfcShrink，这是用来收缩的。
dd2 <- lfcShrink(dds, contrast=contrast, res=dd1,type="ashr") ## 此处 type 还可以是 normal,另外可以使用 coef=2 代替 contrast，但二者不能同时存在，coef=2 可以配合 type= "apeglm" "ashr"
## 这个lfc校正只会对LFC值 影响，不会影响P值 。
plotMA(dd2, ylim=c(-5,5)) ## 此处不用2，为的是让图形收缩显示效果更好。
这样，原先那些波动性很大的基因，就被矫正了。而此时有LogFC以，红色的点为主。
四、差异表达基因分析的结果
这个结果实际上已经通过counts函数获得了，我们不在担心，处理组和对照组完全相反这种情况，因为他内置了参数设定比较组。比如，我们有5个处理组，我们不需要做5次Deseq，我们在results中指定即可。
用summary这个函数，可以看到差异分析的结果,高表达和低表达的比例。低丰度基因所占的比例。
summary(dd2, alpha = 0.05)##
## out of 48485 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up) : 2391, 4.9%
## LFC < 0 (down) : 5361, 11%
## outliers [1] : 0, 0%
## low counts [2] : 14872, 31%
## (mean count < 1)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results再把差异分析的结果转化成data.frame的格式
library(dplyr)
library(tibble)
res <- dd2 %>%
data.frame() %>%
rownames_to_column("gene_id")
!
五、基因ID转换
以前我们从gtf文件转换的, 但是我们需要gtf文件来提取mRNA以及lncRNA，就顺手做了ID转换，而且，mRNA和lncRNA是分别做的Deseq2，这从原理上来讲，是有问题的，Deseq2矫正了测序的深度，而这个深度应该是所有基因算在一起的深度，不应该分开来算。 TGCA数据的标准化以及差异分析 用两个包来转换，得到ENTREZID用于后续分析，得到SYMBOL便于识别。
library(AnnotationDbi)
library(org.Hs.eg.db)
res$symbol <- mapIds(org.Hs.eg.db,
keys=res$gene_id,
column="SYMBOL",
keytype="ENSEMBL",
multiVals="first")
res$entrez <- mapIds(org.Hs.eg.db,
keys=res$gene_id,
column="ENTREZID",
keytype="ENSEMBL",
multiVals="first")
!
有了这个文件，里面有logFC，p值，还有基因名称，我们可以完成GO，KEGG，热图，火山图所有操作。这一部分内容参考 最有诚意的GEO数据库教程
六、制作genelist
我们在那个帖子里面并没有讲GSEA分析，今天来展示一下。原理略过。我们这里还是用Y叔的神包clusterprofier，神包虽好，记得引用。
使用这个包做GSEA，要制作一个genelist，这个是一个向量，他的内容是排序后的logFC值，他的名称是ENTREZID，而这两个我们都是不缺的，在上一步得到的差异结果中。
library(dplyr)
gene_df <- res %>%
dplyr::select(gene_id,log2FoldChange,symbol,entrez) %>%
## 去掉NA
filter(entrez!="NA") %>%
## 去掉重复
distinct(entrez,.keep_all = T)
!
制作genelist的三部曲
## 1.获取基因logFC
geneList <- gene_df$log2FoldChange
## 2.命名
names(geneList) = gene_df$entrez
## 3.排序很重要
geneList = sort(geneList, decreasing = TRUE)
head(geneList) ## 看看这个genelist## 58503 100505879 126393 2277 5104 3952
## 4.264878 4.257538 3.749835 3.513019 3.352578 3.349009七、GSEA分析
以下完成KEGG的GSEA分析
library(clusterProfiler)
gseaKEGG <- gseKEGG(geneList = geneList,
organism = 'hsa',
nPerm = 1000,
minGSSize = 20,
pvalueCutoff = 0.05,
verbose = FALSE)
## 作图，气泡图展示KEGG
library(ggplot2)
dotplot(gseaKEGG,showCategory=12,split=".sign")+facet_grid(~.sign)
这时候，我们看到有一些通路是被激活的，有一些通路是被抑制的。比如Type I diabetes mellitus是被抑制的，我们可以选取单个通路来作图。
把富集的结果转换成data.frame,找到Type I diabetes mellitus的通路ID是hsa04940
gseaKEGG_results <- gseaKEGG@result ## 转换富集结果为数据框
!
使用gseaplot2把他画出来
library(enrichplot)
pathway.id = "hsa04940"
gseaplot2(gseaKEGG,
color = "blue",
geneSetID = pathway.id,
pvalue_table = T)
同样，我们也可以画一个激活的通路
library(enrichplot)
pathway.id = "hsa00510"
gseaplot2(gseaKEGG,
color = "red",
geneSetID = pathway.id,
pvalue_table = T)
八、Pathview展示
我们现在知道cell cycle是被抑制的，如果还想看一下这个通路里面的基因是如何变化的，应该怎么办呢，pathview 可以帮到我们。
library(pathview)
pathway.id = "hsa04940"
pv.out <- pathview(gene.data = geneList,
pathway.id = pathway.id,
species = "hsa")
pathway.id = "hsa00510"
pv.out <- pathview(gene.data = geneList,
pathway.id = pathway.id,
species = "hsa")改变参数，得到另一种构图
library(pathview)
pathway.id = "hsa04940"
pv.out <- pathview(gene.data = geneList,
pathway.id = pathway.id,
species = "hsa",
kegg.native = F)
pathway.id = "hsa00510"
pv.out <- pathview(gene.data = geneList,
pathway.id = pathway.id,
species = "hsa",
kegg.native = F)## [,1] [,2]
## [1,] "101" "209"
## [2,] "40" "209"
## [3,] "41" "209"
## [4,] "187" "210"
## [5,] "95" "210"
## [6,] "100" "211"
## [7,] "39" "211"
## [8,] "40" "212"
## [9,] "41" "212"
## [10,] "101" "212"
## [11,] "39" "214"
## [12,] "100" "214"
## [13,] "95" "217"
## [14,] "187" "217"一眼看过去，都是绿的，说明这个通路确实是被抑制了，还可以在图上缕一缕，哪些是核心分子，一般说来，越往上游越核心。