很有诚意！人人可做的转录组数据下游分析

现在转录组测个序大概是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语言的培训，最应该培训的是基本技能。

但是你有没有发现，虽然图做的这么好看，我们总觉得还欠缺了些什么。这也是目前生物信息分析和实验结合的痛点。

我如果开题，你这一通分析，还是没有让我确定要研究的核心基因。

而这个，是一篇帖子，或者普通培训班不能给予的。

如果未来的生信培训要出点什么彩，这个一定是个方向。

生信技术是通用的，优点就是可被重复，可以被写成教程，但是挖掘能做实验可发文章的核心基因，目前并没有系统教程，类似于玄学。这需要我们长年累月地积累，才能有点感觉，而我觉得这是科研人员做生信分析的核心技能。在这方面，我也是个初学者。

否则，我们就是个能做实验会做分析的技术人员，谈不上一个独立的科研人。