前言:
微博参与话题 #给你四年时间你也学不会生信#
先前的富集分析教程[1]主要是以模式物种人为例子,展开的分析,今天在B站看了孟浩巍视频教程[2],学习新的技能,豁然开朗,欣然记之。
本文主要针对非模式物种,但是有参考基因组可用
1. R包安装及database下载
# non-model, but have the genome
> source("https://bioconductor.org/biocLite.R")
> biocLite("AnnotationHub")
> biocLite("biomaRt")
# load package
> library(AnnotationHub)
> library(biomaRt)
# make a orgDb
> hub <- AnnotationHub::AnnotationHub()
这里以桔小实蝇为例
# fruit fly = bactrocera dorsalis
> query(hub, "bactrocera")
搜索后结果如下:
> query(hub, "bactrocera")
AnnotationHub with 9 records
# snapshotDate(): 2018-04-30
# $dataprovider: ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/
# $species: Bactrocera (Bactrocera)_dorsalis, Bactrocera (Bactrocera)_latifrons, Bactrocera (Dacul...
# $rdataclass: OrgDb
# additional mcols(): taxonomyid, genome, description, coordinate_1_based, maintainer,
# rdatadateadded, preparerclass, tags, rdatapath, sourceurl, sourcetype
# retrieve records with, e.g., 'object[["AH62538"]]'
title
AH62538 | org.Bactrocera_(Bactrocera)_latifrons.eg.sqlite
AH62539 | org.Bactrocera_latifrons.eg.sqlite
AH62542 | org.Bactrocera_(Daculus)_oleae.eg.sqlite
AH62543 | org.Bactrocera_(Dacus)_oleae.eg.sqlite
AH62544 | org.Bactrocera_oleae.eg.sqlite
AH62568 | org.Bactrocera_(Zeugodacus)_cucurbitae.eg.sqlite
AH62569 | org.Bactrocera_cucurbitae.eg.sqlite
AH62581 | org.Bactrocera_(Bactrocera)_dorsalis.eg.sqlite
AH62582 | org.Bactrocera_dorsalis.eg.sqlite
我们选择AH62582 | org.Bactrocera_dorsalis.eg.sqlite
并下载它
> Bactrocera.OrgDb <- hub[["AH62582"]]
如果报错,可能是缺少依赖的安装包,可以按照提示依次下载,两种方法
- install.packages("packages")
- source("https://bioconductor.org/biocLite.R")
biocLite("pacakges")
2. 查看注释信息
> columns(Bactrocera.OrgDb)
[1] "ACCNUM" "ALIAS" "CHR" "ENTREZID" "EVIDENCE" "EVIDENCEALL" "GENENAME"
[8] "GID" "GO" "GOALL" "ONTOLOGY" "ONTOLOGYALL" "PMID" "REFSEQ"
[15] "SYMBOL"
> Bactrocera.OrgDb
OrgDb object:
| DBSCHEMAVERSION: 2.1
| DBSCHEMA: NOSCHEMA_DB
| ORGANISM: Bactrocera dorsalis
| SPECIES: Bactrocera dorsalis
| CENTRALID: GID
| Taxonomy ID: 27457
| Db type: OrgDb
| Supporting package: AnnotationDbi
Please see: help('select') for usage information
# 查看注释信息的每一列
> head(keys(Bactrocera.OrgDb,keytype = "ALIAS"))
[1] "AAA62341.1" "AAA62342.1" "AAA62343.1" "AAA62344.1" "AAF22478.1" "AAL17758.1"
实际上,ALIAS内包含了“omitted 17518 entries”
3. GO富集分析
# 对BP(Biological process)进行富集分析
# 只需将OrgDb数据库替换为我们下载好的非模式物种库即可。
> enrich.go.BP = enrichGO(gene = DEG.gene_symbol,
OrgDb = Bactrocera.OrgDb,
keyType = 'ENTREZID',ont= "BP",
pvalueCutoff = 0.01,
qvalueCutoff = 0.05,
readable = T)
> barplot(enrich.go.BP)
> dotplot(enrich.go.BP)
p_value: 富集显著性,统计显著性要去小于0.01;
q_value: 对p_value的修正,在多次统计检验时,需要有修正值;
q_value一定大于p_value
4. KEGG富集分析
# 只需将OrgDb数据库替换为我们下载好的非模式物种库即可。
> enrichKEGG(gene = DEG.gene_symbol,
OrgDb = Bactrocera.OrgDb,
keyType = 'ENTREZID',
ont = "DO",
pvalueCutoff = 0.01,
qvalueCutofF = 0.05,
readable = T)
5. GO出图解读
Rplot22.png纵轴为GO中每一term,例如Legionellosis;
横轴为GeneRatio,即输入的基因,term在整体基因中所占的百分数;
圆圈大小表示count的数目;
p.adjust:p越小,圆越大,结果越可靠;
网友评论