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.
等待输出