0055-【生物数据库】-如何进行RNA差异基因KEGG注释分析-kobas在线分析

1. 有参物种使用gene ID的方法

1. 差异基因文件准备

只需要用到两列

• ENTREZ_GENE_ID
• logFC

geneNames   ENTREZ_GENE_ID  normalAve   tumorAve    logFC   pValue  qValue
CCL23   6368    95.05964624 5.566645819 -4.066608903    2.07E-31    5.99E-29
COLEC10 10584   1459.366228 83.66298626 -4.122671832    2.11E-31    6.00E-29
FAM189B 10712   383.9435808 1289.852064 1.747953745 2.17E-31    6.08E-29
CDC45   8318    12.20616678 258.9248256 4.38682126  3.59E-31    9.94E-29
RCAN1   1827    11046.97758 2309.590455 -2.257915165    3.90E-31    1.07E-28
N4BP2L1 90634   2644.65753  734.73331   -1.847750259    4.57E-31    1.23E-28
FCN3    8547    6777.184345 389.412555  -4.120767162    5.41E-31    1.44E-28
UHRF1   29128   15.89471347 327.8659692 4.353192433 5.73E-31    1.50E-28
HMMR    3161    25.23294528 407.9486624 4.008655285 8.18E-31    2.12E-28
NEK2    4751    18.88655007 390.7591103 4.36024922  9.48E-31    2.43E-28

选择基因的ID作为输入文件

6368
10584
10712
8318
1827
90634
8547
29128
3161
4751

2. 登陆kobas数据库

网站：http://kobas.cbi.pku.edu.cn/

进入 Gene-list-Enrichment
http://kobas.cbi.pku.edu.cn/anno_iden.php

输入数据类型：

• Fasta Protein Sequence ——蛋白序列
• Fasta Nucleotide Sequence——核酸序列
• Tabular BLAST Output——blast输出的表格
• Entrez Gene ID——基因ID
• UniProtKB AC
• Refseq Protein ID
• Ensembl Gene ID

3. 选择

1. 输入类型选择：Gene ID

2. 物种选择：Homo sapiens (human)

3. 粘贴Gene ID列表

4. 数据库 Clear All取消Pathway、Disease、GO全部选项，只选择KEGG Patway

点击RUN

4. 在线分析完成，输出结果

5. 输出文件说明

统计学检验方法：超几何检验、FIsher精确检验
FDR校正方法：Benjamini and Hochberg，需要补充此方法

##Statistical test method: hypergeometric test / Fisher's exact test
##FDR correction method: Benjamini and Hochberg

输出表格：

• Term KEGG的注释类
• Database 数据库类型
• ID Term对应的ID
• input number 富集到这个Term的输入基因个数
• Background number 数据库中富集到这个通路的总有基因数量
• P-value P值
• Corrected P-Value 校正后P值
• Input 输入的Gene ID，如果多个，以|号分开
• Hyperlink 网页链接

如链接：
http://www.genome.jp/kegg-bin/show_pathway?hsa04512/hsa:3161%09red

图片会将对应的Gene name标志为红色

6. 软件安装准备

由于bioconductor外网链接慢，使用conda的方法安装，同时安装依赖的包

conda install bioconductor-clusterprofiler

7. 画图

# 初始化环境
rm(list=ls())

# 安装软件
#source("https://bioconductor.org/biocLite.R")
#biocLite()
#biocLite("clusterProfiler")
#biocLite("pathview")

# 设置通路
setwd("/home/toucan/Project/001.kegg_map")

# 加载库
library("clusterProfiler")
# 读入文件，不检测name
rt=read.table("input.txt",sep="\t",header=T,check.names=F)
rt

# 构建gene id为行名称的，logFC 
geneFC=rt$logFC
geneFC
gene <- rt$ENTREZ_GENE_ID
gene
names(geneFC)=gene
geneFC

#kegg
# 保存输出文件
# 设定物种，qvalue小于0.05才输出，readable是否输出转换为gene name
kk <- enrichKEGG(gene = gene, organism = "human", pvalueCutoff = 0.05,qvalueCutoff = 0.05)
class(kk)
kk

as.data.frame(kk)

write.table(as.data.frame(kk),file="KEGG.xls",sep="\t",quote=F,row.names = F)


# 生成barplot
pdf(file="KEGG.barplot.pdf")

barplot(kk, drop = TRUE, showCategory = 12)

pdf(file="KEGG.cnetplot.pdf")
# 生成网络图，需要通路描述列、输入gene ID列组成
#cnetplot(kk,categorySize = "geneNum", foldChange = geneFC)

library("pathview")

keggxls=read.table("KEGG.xls",sep="\t",header=T)
# 联网，将map图片下载
for(i in keggxls$ID){
  pv.out <- pathview(gene.data = geneFC, pathway.id = i, species = "hsa", out.suffix = "pathview")}

7. 输出结果

输出富集的表格：

ID  Description GeneRatio   BgRatio pvalue  p.adjust    qvalue  geneID  Count
hsa04110    Cell cycle  19/199  124/7431    5.54E-10    1.37E-07    1.31E-07    8318/7272/890/1870/701/4085/4998/4171/4175/898/23594/1031/4172/4616/8317/4176/4174/9134/993 19
hsa03030    DNA replication 11/199  36/7431 1.29E-09    1.60E-07    1.53E-07    2237/4171/4175/10535/5984/4172/5558/5424/23649/4176/4174    11
hsa03440    Homologous recombination    8/199   41/7431 1.03E-05    0.000849457 0.000811238 146956/8438/5888/7517/5424/641/7516/25788   8
hsa05222    Small cell lung cancer  11/199  93/7431 3.44E-05    0.002135572 0.002039489 1870/898/3910/4616/1282/3655/1284/9134/5743/3915/1163   11

每个通路生成三个文件

• hsa03030.pathview.png
• hsa03030.png
• hsa03030.xml

输出富集的图片

输出伏击通路下载的map：
有差异基因显示，红色为正相关，绿色为负

同时，输出网站原始下载的，无颜色标注

非模式生物

以序列作为输入文件

>seq1
CTAATTTTGATGTAACAATAAGCAAATCCATCTCATTGACATGTCAACTTACCTTAATCTTTAATAAGTG
ATAAAGTCATATGTATGCCAAAAATTGCCTTAGCATTGCGTTATGACCTACCGTTAGTAGATGTCTGATT


>seq2
AGTCTCGAATACAACTTGTTGCTGCGCGGACGCGAATCGCTCAGTACGGACGTCTTGAGCTCGAATCCTC
GGCCATATCTGTGCTCTCGATCGCAGCGTTTGCTAATTCGAAGATCGTGCTAATCGAAGTACCGAGAAAT

注意，物种应选择KO,但会笔记慢

显示：
不应该超过200行的输入文件

If you choose KO, Please input no more than 200 lines at one time.

运行中：

http://kobas.cbi.pku.edu.cn/wait_kobas.php?taskid=180629456069220
Your task is still running, your task id is 180629456069220, you can get the results automatically when the task is finished.
Also you can use the task id to fetch results at the result retrive page in the future.

等待输出