针对某个基因敲除的斑马鱼，在12-25天会逐渐死亡，因此我们选择了8dpf 和16dpf 的由于进行转录组测序，根据差异基因我们绘制了如下的韦恩图，然后我们需要重点关注灰色区域的546个差异基因的部分，GO分析和KEGG分析
一、术语（Terminology）
（1）基因集（gene set）和通路（pathway）：基因集是功能相关的基因的无序集合。忽略基因间的功能关系，可以将通路解释为一个基因集。
（2）Gene Ontology（GO）：GO定义了描述基因功能的概念/类（concepts/classes）。它将功能分为三个方面：
MF: Molecular Function - molecular activities of gene products
CC: Cellular Component - where gene products are active
BP: Biological Process - pathways and larger processes made up of the activities of multiple gene products
GO terms组织在一个有向无环图（directed acyclic graph）中，terms之间的边表示父子关系（parent-child relationship）。
（3）Kyoto Encyclopedia of Genes and Genomes（KEGG）：KEGG是手工绘制的代表分子相互作用和反应网络的路径图（pathway maps）的集合。这些途径涵盖了广泛的生化过程，可分为7大类：新陈代谢、遗传和环境信息处理、细胞过程、机体系统、人类疾病和药物开发（metabolism, genetic and environmental information processing, cellular processes, organismal systems, human diseases, and drug development）。
链接：https://www.jianshu.com/p/d484003dced5
1.利用python找出这546个基因，并写入csv文件中。（因为R学的不好，不知道如何获取546个基因）
from xlrd import open_workbook
import xlrd
import requests
import re
import json
import time
import os
import xlwt
import xlwt
def together_and_apart_from(A,B,C,D,a,b,c,d): #####找到C，D集合之间的交集,并排除A,B的集合, 还需要去掉重复部分。
AandB = []
for i in C:
if i in D:
if i not in A and i not in B:
if i not in AandB:
AandB.append(i)
print(len(AandB))
return AandB
######
4个差异基因文件的路径
path_excel_8_Hv8_C = r'I:\trit1\T202111_8dpf_and_16dpf\revised 更换了个野生型16d\Summary\coverage_gene\T8_HVST8_C_Gene_differential_expression.xls'
path_excel_16_Cv8_C = r'I:\trit1\T202111_8dpf_and_16dpf\revised 更换了个野生型16d\Summary\coverage_gene\T16_CVST8_C_Gene_differential_expression.xls'
path_excel_16_Hv8_H = r'I:\trit1\T202111_8dpf_and_16dpf\revised 更换了个野生型16d\Summary\coverage_gene\T16_HVST8_H_Gene_differential_expression.xls'
path_excel_16_Hv16_C = r'I:\trit1\T202111_8dpf_and_16dpf\revised 更换了个野生型16d\Summary\coverage_gene\T16_HVST16_C_Gene_differential_expression.xls'
读取表格，选中第一列gene_name
workbook_8_Hv8_C = xlrd.open_workbook(path_excel_8_Hv8_C )
sheet1 = workbook_8_Hv8_C.sheet_by_index(0)
col_gene_name_8_Hv8_C = sheet1.col_values(1)
workbook_16_Cv8_C= xlrd.open_workbook(path_excel_16_Cv8_C)
sheet1 = workbook_16_Cv8_C.sheet_by_index(0)
col_gene_name_16_Cv8_C = sheet1.col_values(1)
workbook_16_Hv8_H = xlrd.open_workbook(path_excel_16_Hv8_H )
sheet1 = workbook_16_Hv8_H.sheet_by_index(0)
col_gene_name_16_Hv8_H = sheet1.col_values(1)
workbook_16_Hv16_C= xlrd.open_workbook(path_excel_16_Hv16_C)
sheet1 = workbook_16_Hv16_C.sheet_by_index(0)
col_gene_name_16_Hv16_C = sheet1.col_values(1)
print(len(col_gene_name_8_Hv8_C ),len(col_gene_name_16_Cv8_C),len(col_gene_name_16_Hv8_H),len(col_gene_name_16_Hv16_C))
A=col_gene_name_8_Hv8_C
B=col_gene_name_16_Cv8_C
C=col_gene_name_16_Hv8_H
D=col_gene_name_16_Hv16_C
sets=together_and_apart_from(A, B, C, D, a, b, c, d)
print(sets))
写入表格
style = xlwt.XFStyle()
font = xlwt.Font()
font.name = 'SimSun'
style.font = font
w = xlwt.Workbook(encoding='utf-8')
# 添加个sheet
ws = w.add_sheet('sheet 1', cell_overwrite_ok=True)
i=0
for x in sets:
ws.write(i, 0, x)
i=i+1
w.save('四维韦恩图挑选的交集.xls')
打开保存的xls文件，在第一行插入加入gene_name
转到R
setwd("I:/trit1/T202111_8dpf_and_16dpf/revised 更换了个野生型16d/Summary/coverage_gene")
T16<-read.csv("T16_HVST16_C_Gene_differential_expression.csv") ##1360行，，22列
four_lane<-read.csv("四维韦恩图挑选的交集.csv") ##546行
four_lane
T=merge(T16,four_lane,by = 'gene_name')
## 554 entries, 22 total columns
T1=merge(T16,four_lane,by = 'gene_name',all=TRUE)
## 1,360 entries, 22 total columns
T2=merge(T16,four_lane,by = 'gene_name',all=FALSE)
##554 entries, 22 total columns
所以T和T2得到的交集是所需要的表
T2=merge(T16,T8,by = c('EntrezID','transcript_id'))##两个矩阵以某多列进行合并，共有的序列
2. GO分析
安装需要的包
if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
BiocManager::install("clusterProfiler")
BiocManager::install("topGO")
BiocManager::install("pathview")
BiocManager::install("org.Dr.eg.db")
BiocManager::install("ggnewscale")
library(ggnewscale)
library(clusterProfiler) #用来做富集分析
library(topGO)#画GO图用的
library(pathview) #看KEGG pathway的
library(enrichplot)
library(org.Hs.eg.db)#这个包里存有人的注释文件
library(org.Dr.eg.db)#这个包里存有斑马鱼的注释文件
library(ggplot2)
library(DOSE)
library(GO.db)
在Rstiduo中安装自己的orgDb包
OrgDB注释包
image.png
载入R包并大致查看一下R包的内容
library(org.Dr.eg.db)#这个包里存有斑马鱼的注释文件
columns(org.Dr.eg.db)
image.png
image.png
columns(org.Dr.eg.db)
columns(org.Hs.eg.db)
可用keys（）命令查看一下这个包大致有哪些genes。
image.png
可用select（）命令查一下某个基因对应的GO和Pathway信息
select(org.Dr.eg.db,keys = '30068',columns=c('GO','PATH'))
image.png
做GO分析是不能直接用基因名的，必须得先转化成entre id, 公司的测序数据有ENtrezID 列的数据，但是后来发现很多基因名没有对应的entre id，所以利用返回来数据中的gene_id 重新生成entre id。
获取最新entrez ID，利用gene_id进行转换
ensembl_gene_id=T$gene_id
id <-bitr(ensembl_gene_id, fromType = "ENSEMBL",
toType = c("ENTREZID"),
OrgDb = org.Dr.eg.db,drop = FALSE )
ENTREZ_ID = id$ENTREZID
获取entrez ID
gene_name=T[,1] #基因名
log2foldchange=T[,17] #log2foldchange列
ENTREZ_ID = id$ENTREZID #EntrezID列
T的数据格式
image.png
readable=TRUE代表将基因ID转换为gene symbol
BP层面上的富集分析：
go_bp<-enrichGO(gene =ENTREZ_ID ,OrgDb = org.Dr.eg.db, keyType='ENTREZID', ont = "BP", pAdjustMethod = "BH",pvalueCutoff = 0.05, qvalueCutoff = 0.05, readable=TRUE)
write.csv(go_bp@result,"go_bp.csv")
dim(go_bp)# 148 9
CC层面上的富集分析：
go_cc<-enrichGO(gene = ENTREZ_ID,OrgDb = org.Dr.eg.db,keyType = 'ENTREZID', ont = "CC", pAdjustMethod = "BH",pvalueCutoff = 0.05, qvalueCutoff = 0.05)
把结果导出保存
write.csv(go_bp@result,"go_bp.csv")
dim(go_cc)##33 9
go_MF <- enrichGO(gene =ENTREZ_ID, OrgDb= org.Dr.eg.db, keyType = 'ENTREZID', ont = "MF",pAdjustMethod = "BH",pvalueCutoff = 0.05,qvalueCutoff = 0.05)
write.csv(go_MF@result,"go_mf.csv")
dim(go_MF)##6 9
all包括三个组分的分析
go_all<-enrichGO(gene = ENTREZ_ID,OrgDb = org.Dr.eg.db,keyType = 'ENTREZID', ont = "ALL", pAdjustMethod = "BH",pvalueCutoff = 0.05, qvalueCutoff = 0.05,readable = TRUE)
write.csv(go_all@result,"go_all.csv")
dim(go_all)# 187 10
sum(go_all$ONTOLOGY=="BP") #148
sum(go_all$ONTOLOGY=="CC") ##33
sum(go_all$ONTOLOGY=="MF")##6
go_cc 保存结果image.png
reduce redundancy of enriched GO terms
GO以parent-child结构组织，因此父术语与所有子术语有很大比例的重叠。这可能导致冗余的结果。为了解决这一问题，clusterProfiler实现了简化方法simplify，以减少enrichGO 和gseGO产生的冗余的GO术语。函数内部称为GOSemSim (Yu et al. 2010)，用于计算GO项之间的语义相似度，并通过保留一个代表性项来去除高度相似的项。
egosimp <- simplify(go_bp,cutoff=0.7,by="p.adjust",select_fun=min,measure="Wang")
write.csv(egosimp@result,"gosimp.csv")
head(egosimp);
dim(egosimp)# 57 9
方法1:基于它们的共有父条目的注释统计,计算语义相似性得分
包含Resnik、Lin、Jiang 和Schlicker四种方法
方法2:基于GO图形结构,Wang
进行GO terms集的相似性分析时一般采取基于Resnik和Lin两种方法的综合方法,简称为simRel方法
链接：https://www.jianshu.com/p/d484003dced5
可视化
做成气泡图形式
dotplot(object,x="GeneRatio",color="p.adjust",showCategory=10,
size=NULL,split=NULL,font.size=12,title="",...)
横轴为GeneRatio，代表该GO term下富集到的基因个数占列表基因总数的比例
纵轴为富集到的GO Terms的描述信息，showCategory指定展示的GO Terms的个数
dotplot(egoall,showCategory=10)
dotplot(go_bp,show,Category=20)
dotplot(go_cc)
dotplot(go_MF)
做成条形图形式
barplot(go_bp, showCategory = 12,title="The GO_BP enrichment analysis of all DEGs ", drop=TRUE)
三个图汇聚在一块
dotplot(go_all,title='Top5 GO terms of each sub-class',showCategory=5,split='ONTOLOGY')+facet_grid(ONTOLOGY~.,scale="free")
image.png
前五条通路的共同基因的circle图
对于基因和富集的GO terms之间的对应关系进行展示图中灰色的点代表基因，黄色的点代表富集到的GO terms,如果一个基因位于一个GO Terms下，则将该基因与GO连线,黄色节点的大小对应富集到的基因个数，默认画top5富集到的GO terms
geneList=log2foldchange
names(geneList)=gene_name
geneList=sort(geneList,decreasing = T)
#(3)展示top5通路的共同基因，要放大看。
#Gene-Concept Network
对于基因和富集的GO terms之间的对应关系进行展示，如果一个基因位于一个GO Terms下，则将该基因与GO连线，用法如下
cnetplot(go_bp, showCategory=5)
cnetplot(go_bp, categorySize="pvalue", foldChange=geneList,colorEdge = TRUE)
#圆形布局，给线条上色
cnetplot(go_bp, foldChange=geneList, circular = TRUE, colorEdge = TRUE)
图中灰色的点代表基因，黄色的点代表富集到的GO terms, 默认画top5富集到的GO terms, GO 节点的大小对应富集到的基因个数。更多用法和细节请参考官方文档。
通路之间的关系
goplot(go_bp)
Heatmap-like functional classification
基因展示的heatmap
heatplot(go_bp, foldChange = geneList)
太多基因就会糊。可通过调整比例或者减少基因来控制。
pdf("heatplot.pdf",width = 14,height = 5)
heatplot(ego_bp,foldChange = geneList)
dev.off()
readable=TRUE 代表将基因ID转换为gene symbol
go_bp<-enrichGO(gene =ENTREZ_ID ,OrgDb = org.Dr.eg.db, keyType='ENTREZID', ont = "BP", pAdjustMethod = "BH",pvalueCutoff = 0.05, qvalueCutoff = 0.05, readable=TRUE)
cnetplot(go_bp, categorySize="pvalue", foldChange=geneList,colorEdge = TRUE)
加上readable=TRUE 就可以显示基因的名字image.png
GO有向无环图
矩形代表富集到的top10个GO terms, 颜色从黄色过滤到红色，对应p值从大到小。
plotGOgraph(go_bp)
emapplot(go_bp,showCategory = 30)
对于富集到的GO terms之间的基因重叠关系进行展示,每个节点是一个富集到的GO term，默认画top30个富集到的GO terms。节点大小对应该GO terms下富集到的基因个数，节点的颜色对应p.adjust的值，红色小蓝色大，如果两个GO terms的差异基因存在重叠，说明这两个节点存在overlap关系，用线条连接起来
emapplot(egoMF,showCategory=10)
对于富集到的GO terms之间的基因重叠关系进行展示，如果两个GO terms系的差异基因存在重叠，说明这两个节点存在overlap关系，在图中用线条连接起来，用法如下
emapplot(go_bp,showCategory = 30)会报错，需要经过以下代码运行下
x2<-pairwise_termsim(go_bp)
emapplot(x2,showCategory = 30)
每个节点是一个富集到的GO term, 默认画top30个富集到的GO terms, 节点大小对应该GO terms下富集到的差异基因个数，节点的颜色对应p.adjust的值，从小到大，对应蓝色到红色。
emapplot(x2,showCategory = 10)
image.png
upsetplot(go_bp)
upset展现多个数据集合的交互关系的图形，类似于传统的venn图
山脊线图 ridgeline plot for expression distribution of GSEA result
KEGG分析
KEGG 分析
search_kegg_organism("human", by="common_name")
search_kegg_organism("zebrafish", by="common_name")
gene_kegg<-enrichKEGG(gene =ENTREZ_ID ,organism = 'dre',keyType='kegg', pAdjustMethod = "BH",pvalueCutoff = 0.05, qvalueCutoff = 0.05,use_internal_data = FALSE)
dim(gene_kegg)
dotplot(gene_kegg)
barplot(gene_kegg)
enrichMap(gene_kegg)
cnetplot(gene_kegg, showCategory=5)
#将ENTREZID转化为可读的gene symbol
eKEGG <- setReadable(gene_kegg, OrgDb = org.Dr.eg.db, keyType="ENTREZID")
cnetplot(eKEGG, showCategory=5)
转换位基因名字后
gene.kegg <- bitr_kegg(ENTREZ_ID,fromType="ncbi-geneid",toType="kegg",organism='dre')
head(gene.kegg)
ekegg <- enrichKEGG(ENTREZ_ID, organism='dre',keyType="ncbi-geneid",pvalueCutoff=0.05,pAdjustMethod='BH',qvalueCutoff=0.2, minGSSize=10,maxGSSize=500,use_internal_data=F)
write.csv(ekegg@result,"ekegg.csv")
head(ekegg,1);dim(ekegg)
