转录组学习八（功能富集分析）

转录组学习一（软件安装）
转录组学习二（数据下载）
转录组学习三（数据质控）
转录组学习四（参考基因组及gtf注释探究）
转录组学习五（reads的比对与samtools排序）
转录组学习六（reads计数与标准化）
转录组学习七（差异基因分析）
转录组学习八（功能富集分析）

任务

• 选择p<0.05而且abs(log2FC)大于1的基因为显著差异表达基因集，对这个基因集用R包做KEGG/GO超几何分布检验分析。
• 把表达矩阵和分组信息分别作出cls和gct文件，导入到GSEA软件分析。

<font color=orange>GO富集分析</font>

一、bioconductor注释数据库的探究

JIMMY_数据包library(org.Hs.eg.db)简介

• 官网上已有19个注释数据库，网址bioconductor.org/packages/release/ 看了一下植物只有一个拟南芥的。
• 如果不在注释数据库里，使用AnnotationHub来构建自己的Org.db。

library(AnnotationHub)
hub <- AnnotationHub()
# 可以用query（）函数来查找你要的物种注释信息
# 选择的格式是OrgDb.
query（hub, "Solanum lycopersicum"）
sl <- hub[["AH55774"]] 

• 对ORG.Mm.eg.db进行简单的探究

library(org.Mm.eg.db)
keytypes(org.Mm.eg.db) ##查看有哪些数据类型，包含着各大主流数据库的数据。
###用select函数，就可以把任意公共数据库的数据进行一一对应。
### keys是原始的ID，columns是转换之后的ID，keytype是要指定的原始ID类型
select(org.Mm.eg.db,keys = "ENSMUSG00000031762",columns = c("SYMBOL","GENENAME","UNIGENE","REFSEQ"),keytype = "ENSEMBL")

• 对ENSEMBL-ID进行转换

diff_gene_DESeq_raw <- subset(res_deseq, padj<0.05 & (log2FoldChange > 1 | log2FoldChange < -1))
diff_gene_DESeq_name <- row.names(diff_gene_DESeq_raw)
diff_gene_DESeq_transID<- select(org.Mm.eg.db, keys= diff_gene_DESeq_name, columns= c("SYMBOL", "GENENAME", "UNIGENE", "REFSEQ"), keystype="ENSEMBL")

image

• 结果可知几大数据库的基因对应关系，其中有部分的数据是未知的？还不知道是什么原因，查找了下原始的小鼠gtf文件，发现注释的基因是最近的？所以说可能是数据还不全导致的吧。

二、基因GO分析：

利用clusterProfiler的R包进行GO分析

enrichGO(gene, OrgDb, keytype = "ENTREZID", ont = "MF",
  pvalueCutoff = 0.05, pAdjustMethod = "BH", universe, qvalueCutoff = 0.2,
  minGSSize = 10, maxGSSize = 500, readable = FALSE, pool = FALSE)

• 基本参数：gene:即基因ID；OrgDb:物种注释数据；keytype:ID的类型 ont 三个层面来阐述基因功能，生物学过程（BP），细胞组分（CC），分子功能（MF）；pAdjustMethod:P值校正方法
• 进行GO分析的基本代码

ego <- enrichGO(gene = row.names(diff_gene_deseq2),  OrgDb = org.Mm.eg.db, keytype = "ENSEMBL", ont = "MF")

###气泡图
dotplot(ego, font)

### 网络图
enrichMap(ego, vertex.label.cex=1.2, layout=igraph::layout.kamada.kawai)

###GO图
plotGOgraph(ego)

image
image
image

三、基因KEGG分析：

• KEGG介绍：KEGG数据库（Kyoto Encyclopedia of Genes and Genomes）京都基因与基因组百科全书：KEGG PATHWAY数据库由一系列经人手绘制而成KEGG代谢路径图的构成，以代表关于代谢以及其他细胞与生物机能的实验成果。每一张路径图俱包括了一个分子互动与反应的网络图，为连系基因组内之基因及路径所示过程生成的基因产物（以蛋白质为主）而设。行使相同功能的基因聚在一起

• 目前KEGG等数据库，收录的是已有的研究结果。但这些pathway的信息，远没有到达完善的水准。大部分通路只是了解1个大概的调控途径，而中间有什么转录因子参与、是否还有其他代谢物的生成，都是不知道的。这些通路的完整性，也会影响pathway富集分析结果。例如，基因A发生变化了，看起来下游基因没有变化。也许是还有其他的调控在起作用，只是这些调控作用现在还不知道而已。

• pathway 和 GO富集分析结果的解读，应该从生物学意义的角度出发，P value 和 Q value只是个参考而已，那些不显著的通路也值得解读（从功能注释的角度解读，而不是从富集分析的角度解读）。只要结果可以解释，有意义，不用太迷信P value。

• 缩写查看官网：http://www.genome.jp/kegg/catalog/org_list.html

• KEGG基本过程

diff_gene_deseq2_transID_kegg <- diff_gene_deseq2_transID[,4]
ekegg <- enrichKEGG(diff_gene_deseq2_transID_kegg,keyType = "kegg",organism = "mmu",pvalueCutoff = 0.05,pAdjustMethod = "BH",qvalueCutoff = 0.1)
### 画气泡图：
dotplot(ekegg,font.size=8)

### 显示通路图
browseKEGG(ekegg,'mmu01100')

image
image

<font color=orange>基因集富集分析GSEA</font>

参考文章GSEA分析是个什么鬼？（上）, GSEA分析是个什么鬼？（下）。文章将GSEA分析做了详细的揭示，目前仅对看懂的部分做记录，知道有这个分析方法，以后有需要再做详细学习吧

一、基因富集分析概念

• 基因富集分析(Gene set enrichment analysis, GSEA)的基本思想是使用预定义的基因集（通常是来自于功能注释或先前实验的结果），将基因按照在两类样本中的差异表达程度排序，然后检验预先设定的基因集合是否在这个排序表的顶端或者底端富集。基因集合的富集分析是检测基因集合而不是单个基因的表达变化，因此可以包含这些细微的表达变化，预期得到更为理想的结果。
  • 分析的基因集合而不是单个基因
  • 将基因与预定义的基因集进行比较
  • 富集分析

二、与通常富集方法GO和KEGG的比较：

• GO和 KEGG 侧重于比较两组间的基因表达差异，集中关注少数几个显著上调或下调的基因，这容易遗漏部分差异表达不显著却有重要生物学意义的基因，忽略一些基因的生物特性、基因调控网络之间的关系及基因功能和意义等有价值的信息。

• GSEA分析不需要指定明确的差异基因阈值，算法会根据实际数据的整体趋势， 为研究者们提供了一种合理地解决目前芯片分析瓶颈问题的方法,即使在没有先验经验存在的情况下也能在表达谱整体层次上对数条基因进行分析,从而从数理统计上把表达谱芯片数据与生物学意义很好地衔接起来,使得研究者们能够更轻松、更合理地解读芯片结果。

三、通常做差异分析设定阈值与后续KEGG与GO分析的问题**：

• 差异基因列表后，都会在此之上提供给客户Pathway 以及GO 富集分析。这里这些筛选出来的差异基因本身就会存在一些问题。

• 在实际应用于生物学高通量数据时，对于差异基因检出的阈值，异常的敏感，客户需要给出差异基因的一个明确的定义(阈值)，例如abs(FC) ≧2.0 & p ≦ 0.05。这种一刀切的阈值，对于发现真正的生物学效应，许多时候是一种障碍，因为实际通过芯片观测到的RNA 表达变化，往往是层层的负反馈调控后的结果，并且不同组织对于表达差异的敏感度是不同的：在神经递质系统内，一个1.2 倍的表达差异即可能产生及其显著的效应。

四、GSEA富集过程的基本步骤

1. 计算富集分数（Enrichment Score）
2. 估计富集分数的显著性水平
3. 矫正多重假设检验
   image
   image

五、基本GSEA分析过程

GSEA_genelist <- diff_gene_deseq2_raw$log2FoldChange ### 对diff_gene的结果进行分析
names(GSEA_genelist)<- rownames(diff_gene_deseq2_raw) ### 设置名字
GSEA_genelist<- sort(GSEA_genelist,decreasing = TRUE) ### 排序


gsem_gene <- gseGO(geneList = GSEA_genelist,OrgDb = org.Mm.eg.db, keyType = "ENSEMBL", ont = "MF")

gseaplot(gsem_gene,geneSetID = "GO:0000977")

image
看不懂在说啥，以后再慢慢研究这一类的图吧。

总结：从2017年10月7日~12月4号，软件安装——数据下载——原始数据的质控——参考基因组与注释GTF文件的探究——READS比对，排序——计数，标准化——差异基因分析——功能富集分析。Linux基本操作、shell脚本编写、Perl脚本编写、软件参数的具体含义了解。两个月的时间终于跟着大神们的步伐将转录组的流程给学习了一遍。途中遇到了不少艰难的事，也花了不少深夜与周末的时间。好在终于把一个个知识点给攻克，一章章的任务分析给坚持了下来。仍然有许多待学习的地方，比如一些软件的进阶参数的选择，结果的更加准确解读，R语言的语法、各种可视化图片的绘制以及如何解读，还有越来越觉得重要的一个大坑：统计学背景知识。这些都将是后续学习的方向。这不是结束，以后还会根据学习到的各种新知识来更好的完善这个分析流程大框架。加油吧。

参考文章

1. 【基因富集分析_学习笔记】https://mp.weixin.qq.com/s?__biz=MzIwNTEwMTUyOQ==&mid=2649693906&idx=1&sn=341682dad10a9b52f3290239042c30f5&chksm=8f2dbe64b85a3772d3bc439498560ec22638783cb321be71e7ebb383a4de0186b3ea9b384475&scene=21#wechat_redirect
2. 【PANDA姐的转录组入门(8)：差异基因结果注释】https://mp.weixin.qq.com/s?__biz=MzIwNTEwMTUyOQ==&mid=2649694917&idx=1&sn=a318f8cf98f306d46963986011c73600&chksm=8f2d8273b85a0b65270a986f7bb28e8efa7fd15309505a5a026bbedaafaff7e30e4d4f13dd02&scene=21#wechat_redirect
3. 【(伪）从零开始学转录组（8）：富集分析】https://mp.weixin.qq.com/s?__biz=MzI1MjU5MjMzNA==&mid=2247484528&idx=1&sn=3af297d4163a70b049f861fc28c85bc2&chksm=e9e02dd1de97a4c79ab04a2c012ee6fe4b0aae5e1df59bdf2bf6b19e6ac1cb08ce379a6694e8&scene=21#wechat_redirect
4. 【差异基因结果注释】http://www.jianshu.com/p/4910d7cec5c8
5. 【GSEA分析是个什么鬼？（上）】https://mp.weixin.qq.com/s?__biz=MzAwMzY4MTYxNw==&mid=2655753566&idx=2&sn=5b5b2c93a7618a69da2cbc6638f03da0&chksm=80884960b7ffc076af53ae74caadb5dbb25d240c31660792e8727964d0177d6a17af7ca5fc5c&mpshare=1&scene=1&srcid=0816ADpKId3sPzgbYfubrFCf#rd
6. 【GSEA是个什么鬼?(下)】https://mp.weixin.qq.com/s?__biz=MzAwMzY4MTYxNw==&mid=2655754973&idx=1&sn=3b87d5cb8ddd2d5d77e413e9a87342da&chksm=808846e3b7ffcff5a6b41985b707f52170f20eabe15fc43264b3d14a3ccf4100263789eab856&mpshare=1&scene=1&srcid=0816gHxusewlJeILw0fWxgi3#rd
7. 【RNA-seq结果图片如何解读？（第一弹）】http://mp.weixin.qq.com/s/OFuP7nGGM3V9ghZ6lI1QuA
8. 【RNA-seq中GO、KEGG结果图如何解读】http://mp.weixin.qq.com/s/UowQnL4bD7QUFHIXQ_JQKQ