RNA-seq的下游分析，人人都要会！
现在转录组测个序大概是1000-2000块钱一个样，这个波动的范围取决于各个公司提供的服务，以及你们那个地区生信的普及程度。既然这么便宜，那么每个看到明确现象的实验团队都改尝试一下RNA-seq，说不定就给课题开了新的思路。
转录组测序的分析分为上游分析和下游分析，简单区分就是，你有没有服务器。如果有，那就把上游分析给包了，这在以前不可想象，但是因为生信技能树这样的团体存在，推动了我国生物信息技术的普及，让生物信息不在遥不可及，而这本是国之重器们该做的事情。
假如你没有服务器，也不要紧，Y叔的出现，使得下游分析变得十分简单。仅仅使用R语言，使用Y叔的神包clusterprofiler，坐在家里喝着咖啡，也可以搞定所有下游分析。
我们今天的任务就是展示一下这个过程，抛砖引玉。
首先，加载数据。
load(file = "expr_df.Rdata")
这个数据就是测序后的counts矩阵，跟公司要到这个数据就可以脱离服务器自己往下走。如果没有数据也不要紧，在TCGA的教程中，我们也获得过这样的数据。
从GDC下载TCGA肿瘤数据库的数据
把GDC下载的多个TCGA文件批量读入R
我们来看一看这个数据：
他由一列组成，第一列是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("con","treat"),each=3)
metadata$sample <- relevel(factor(sample),"con")
构建dds对象
这一步由DESeqDataSetFromMatrix这个函数来完成，他需要输入我们的表达矩阵，制作好的metadata，还要制定分组的列，在这里是sample，最后一个tidy的意思是，我们第一列是基因ID，需要自动处理。
library(DESeq2)
## 第一列有名称，所以tidy=TRUE
dds <-DESeqDataSetFromMatrix(countData=expr_df,
colData=metadata,
design=~sample,
tidy=TRUE)
看一下有多少行呢？58734行
nrow(dds)
查看一下有行名么？有的，如果没有就不要往下做了，应该是你的数据不是data.frame格式
rownames(dds)
数据过滤
别看有这么多行，不是每一个就要都要表达的，如果一行基因在所有样本中的counts数小于等于1，我们就把他删掉，实际上，不做这一步，对差异分析的结果没有影响，可能会对GSEA的结果有影响。
dds <- dds[rowSums(counts(dds))>1,]
这时候我们再来看他的行数已经变成了27134
nrow(dds)
样本聚类
有时候，我们会把实验样本的的标签搞错，比如，明明是加药组，但是有一组没忘记加药，如果对他聚类，他会被划归为正常组，这个是需要我们知道的。
用vst来标化数据，实际上还有rlog方法，或者就是log2的方法，官网推荐< 30个样本用rlog，大于30个样本用vst，速度更快，这里我们不要计较那么多了，就用vst，因为真实的TCGA数据，样本往往大于30个。
vsd <- vst(dds, blind = FALSE)
再用dist这个函数计算样本间的距离
sampleDists <- dist(t(assay(vsd)))
用hclust来进行层次聚类
hc <- hclust(sampleDists, method = "ward.D2")
然后画图
plot(hc, hang = -1)
从图上看，两类样本分开了，不需要特殊处理, 这个图可能没什么用，但是如果要给老板呈现，我觉得态度要端正，可以尝试一下别的方法。
library(factoextra)
res <- hcut(sampleDists, k = 2, stand = 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)
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 = "ENSG00000172183", intgroup=c("sample"))
这个功能本质上没有什么用，但是，可以提前确定实验的质量。比如，你的两组是敲减某个基因以及对照组，通过制定那个基因作图，就可以看出实验有没有成功，如果这个基因没有任何改变，也可以不用往下做了。回去重新做实验送样本吧。
我们说过，这种图自己看就可以，给老板看就算了，你得美化一下，那么他里面有个内置的参数returnData可以返回作图数据。
一旦返回数据，我们就可以用ggplot2自己简单画一下。
plotdata <- plotCounts(dds, gene = "ENSG00000172183", intgroup=c("sample"),returnData = T)
library(ggplot2)
ggplot(plotdata,aes(x=sample,y=count,col=sample))+
geom_point()+
theme_bw()
现在看起来平淡无奇，如果样本多，可以画出点图配boxplot，如果是配对样本，那么还可以画出配对的图。
有一点要强调一下，vst，以及log2的标准化，跟标准化的counts不是一个概念。前者是为了以后的聚类，热土，PCA分析，比如，我们计算样本间的距离就是用的vst 标化的数据，而标准化的counts是为了差异作图，你看纵坐标就会发现，他的数值一般很大。
LogFC的矫正
这一步，对于依赖logFC变化值的分析，很重要，比如GSEA分析。我们画一个MAplot图，看图说话
contrast <- c("sample", "treat", "con")
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)
plotMA(dd2, ylim=c(-2,2))
这样，原先那些波动性很大的基因，就被矫正了。而此时有LogFC以，红色的点为主。
差异分析的结果
这个结果实际上已经通过counts函数获得了，我们不在担心，处理组和对照组完全相反这种情况，因为他内置了参数设定比较组。比如，我们有5个处理组，我们不需要做5次Deseq，我们在results中指定即可。
用summary这个函数，可以看到差异分析的结果,高表达和低表达的比例。低丰度基因所占的比例。
summary(dd2, alpha = 0.05)
再把差异分析的结果转化成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数据库教程
这个帖子目前已经有接近4000次访问，这在我们这样一个小号是不容易的，靠的是真诚。
制作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)
看一下这个genelist，增加感性的理解
head(geneList)
GSEA分析
完成KEGG的GSEA分析
library(clusterProfiler)
## 没有富集到任何数据
gseaKEGG <- gseKEGG(geneList = geneList,
organism = 'hsa',
nPerm = 1000,
minGSSize = 20,
pvalueCutoff = 0.05,
verbose = FALSE)
作图展示富集分布图
library(ggplot2)
dotplot(gseaKEGG,showCategory=12,split=".sign")+facet_grid(~.sign)
这时候，我们看到有一些通路是被激活的，有一些通路是被抑制的。比如Cell cycle是被抑制的，我们可以选取单个通路来作图。
把富集的结果转换成data.frame,找到Cell cycle的通路ID是hsa04110
gseaKEGG_results <- gseaKEGG@result
使用gseaplot2把他画出来
library(enrichplot)
pathway.id = "hsa04110"
gseaplot2(gseaKEGG,
color = "red",
geneSetID = pathway.id,
pvalue_table = T)
也可以画出一个激活的
pathway.id = "hsa04060"
gseaplot2(gseaKEGG,
color = "red",
geneSetID = pathway.id,
pvalue_table = T)
pathview 展示
我们现在知道cell cycle是被抑制的，如果还想看一下这个通路里面的基因是如何变化的，应该怎么办呢，pathview 可以帮到我们。
library(pathview)
pathway.id = "hsa04110"
pv.out <- pathview(gene.data = geneList,
pathway.id = pathway.id,
species = "hsa")
改变一下参数，可以得到另外一种构图
library(pathview)
pathway.id = "hsa04110"
pv.out <- pathview(gene.data = geneList,
pathway.id = pathway.id,
species = "hsa",
kegg.native = F)
一眼看过去，都是绿的，说明这个通路确实是被抑制了，还可以在图上缕一缕，哪些是核心分子，一般说来，越往上游越核心。
总结
写到这里，GEO的分析，TCGA的基本分析，RNA-seq的基本分析都写完了。这几个帖子可以把大部分的培训班给搞定。里面的内容，只要会一点R语言，就可以重复。说到底，R语言的培训，最应该培训的是基本技能。
但是你有没有发现，虽然图做的这么好看，我们总觉得还欠缺了些什么。这也是目前生物信息分析和实验结合的痛点。
我如果开题，你这一通分析，还是没有让我确定要研究的核心基因。
而这个，是一篇帖子，或者普通培训班不能给予的。
如果未来的生信培训要出点什么彩，这个一定是个方向。
生信技术是通用的，优点就是可被重复，可以被写成教程，但是挖掘能做实验可发文章的核心基因，目前并没有系统教程，类似于玄学。这需要我们长年累月地积累，才能有点感觉，而我觉得这是科研人员做生信分析的核心技能。在这方面，我也是个初学者。
否则，我们就是个能做实验会做分析的技术人员，谈不上一个独立的科研人。
