测序了，然后呢（四）OrgDb与富集分析

生物_医药_科研 2019-01-22

展开全文

今天是生信星球陪你的第256天

大神一句话，菜鸟跑半年。我不是大神，但我可以缩短你走弯路的半年~

就像歌儿唱的那样，如果你不知道该往哪儿走，就留在这学点生信好不好~

这里有豆豆和花花的学习历程，从新手到进阶，生信路上有你有我！

豆豆写于19.1.20-21
为什么要搞全基因组测序（一）
测序了，然后呢（二）基因功能注释
测序了，然后呢（三）功能注释整合流程
要进行GO或者KEGG富集分析，就需要知道每个基因对应什么样的GO/KEGG分类，OrgDb就是存储不同数据库基因ID之间对应关系，以及基因与GO等注释的对应关系的 R 软件包
如果自己研究的物种不在http:///packages/release/BiocViews.html#___OrgDb 之列，很大可能就需要自己构建OrgDb，然后用clusterProfiler分析
所以本次的内容都是基于非模式生物
--------------------------------------------------------------------------
发推送的时候花花和豆豆就在看我们十二期的学员们在简书第七天的感想，看的好感动，谢谢你们的支持，看到你们真的没有白付出，我们很欣慰。
今天其实很累，但是写推送让我重回活力，我爱写推送，我爱和你们交流
感谢所有关注我们的小伙伴，我们会更加努力，招财猫送你们好运！！！

第一类：利用AnnotationHub得到org.db

下面我以棉铃虫为例

首先下载并加载AnnotationHub

options('repos'= c(CRAN='https://mirrors.tuna./CRAN/'))
options(BioC_mirror='http://mirrors.ustc.edu.cn/bioc/')

if (!requireNamespace('BiocManager', quietly = TRUE))
    install.packages('BiocManager')
BiocManager::install('AnnotationHub', version = '3.8')
library(AnnotationHub)

然后加载现有的物种database

hub <>#这一步受限于网速，不成功的话多试几次
# 调用图形界面查看物种
display(hub) 
# 或者根据物种拉丁文名称查找
query(hub,'helicoverpa')
#   'object[['AH66950']]' 

            title                                      
  AH66950 | org.Helicoverpa_armigera.eg.sqlite         
  AH66951 | org.Heliothis_(Helicoverpa)_armigera.eg....
# 这里AH66950是我们需要的
# 然后下载这个sqlite数据库
ha.db <>'AH66950']]
#查看前几个基因（Entrez命名）
head(keys(ha.db))
#查看包含的基因数
length(keys(ha.db)) 
#查看包含多少种ID
columns(ha.db)
#查看前几个基因的ID
select(ha.db, keys(ha.db)[1:3], 
       c('REFSEQ', 'SYMBOL'), #想获取的ID
       'ENTREZID')
#得到结果
  ENTREZID         REFSEQ SYMBOL
1  9977712 YP_004021052.1   COX1
2  9977713 YP_004021053.1   COX2
3  9977714 YP_004021054.1   ATP8
#保存到文件
saveDb(ha.db, 'Harms-AH66950.sqlite')
#之后再使用直接加载进来
maize.db <>'Harms-AH66950.sqlite')

参考：https://mp.weixin.qq.com/s/lHKZtzpN2k9uPN7e6HjH3w

第二类：利用AnnotationForge得到org.db

通过上次的推送（功能注释整合流程），我们知道了怎样利用eggnog-mapper去人工构建注释、富集桥梁

上次选取的是一个病毒的小例子，这次可以用芝麻（Sesame）做演示

先下载蛋白序列

wget http://www./SesameFG/BLAST_search/G608_contig_2014-08-29.FgeneSH.pep.rar
# 解压后上传

因为统计了下有38406条序列，因此使用diamond来进行功能注释

# 前提是自己安装好eggnog-mapper并且下载好相应的数据库
emapper.py -m diamond \
           -i sesame.fa \
           -o diamond \
           --cpu 19
# 得到如下信息，然后进行处理，只保留表头query_name这一行的注释信息，去掉头尾的# 等信息
sed -i '/^# /d' diamond.emapper.annotations 
sed -i 's/#//' diamond.emapper.annotations

关于结果解释：https://github.com/jhcepas/eggnog-mapper/wiki/Results-Interpretation
其中关于COG的介绍：倒数第二列
COG functional categories: COG functional category inferred from best matching OG 【会给出一个大写字母，每一个大写字母都有自己的解释：COG_explanation】

STEP1：自己构建的话，首先需要读入文件

egg_f <>'diamond.emapper.annotations'
egg <>'\t')
egg[egg=='']<>NA #这个代码来自花花的指导(将空行变成NA，方便下面的去除)

STEP2: 从文件中挑出基因query_name与eggnog注释信息

gene_info <- egg %>%
        dplyr::select(GID = query_name, GENENAME = `eggNOG annot`) %>% na.omit()

STEP3-1：挑出query_name与GO注释信息

gterms <- egg %>%
        dplyr::select(query_name, GO_terms) %>% na.omit()

STEP3-2：我们想得到query_name与GO号的对应信息

# 先构建一个空的数据框(弄好大体的架构，表示其中要有GID =》query_name，GO =》GO号， EVIDENCE =》默认IDA)
# 关于IEA：就是一个标准，除了这个标准以外还有许多。IEA就是表示我们的注释是自动注释，无需人工检查http://wiki./index.php/Inferred_from_Electronic_Annotation_(IEA)
# 两种情况下需要用IEA：1. manually constructed mappings between external classification systems and GO terms； 2.automatic transfer of annotation to orthologous gene products.
gene2go <>
                         GO = character(),
                         EVIDENCE = character())

# 然后向其中填充：注意到有的query_name对应多个GO，因此我们以GO号为标准，每一行只能有一个GO号，但query_name和Evidence可以重复
for (row in 1:nrow(gterms)) {
        gene_terms <>'GO_terms'], ',', simplify = FALSE)[[1]]  
        gene_id <>'query_name'][[1]]
        tmp <>
                              GO = gene_terms,
                              EVIDENCE = rep('IEA', length(gene_terms)))
        gene2go <>
    }

STEP4-1: 挑出query_name与KEGG注释信息

gene2ko <- egg %>%
        dplyr::select(GID = query_name, KO = KEGG_KOs) %>%
        na.omit()

STEP4-2: 得到pathway2name, ko2pathway

这一步不需要管代码什么意思，只需要知道它可以帮我们得到以上两个文件就好

if(F){
    # 需要下载 json文件(这是是经常更新的)
    # https://www./kegg-bin/get_htext?ko00001
    # 代码来自：http://www./course/225/task/4861/show
    library(jsonlite)
    library(purrr)
    library(RCurl)

    update_kegg <>function(json = 'ko00001.json') {
        pathway2name <>
        ko2pathway <>

        kegg <>

        for (a in seq_along(kegg[['children']][['children']])) {
            A <>'children']][['name']][[a]]

            for (b in seq_along(kegg[['children']][['children']][[a]][['children']])) {
                B <>'children']][['children']][[a]][['name']][[b]] 

                for (c in seq_along(kegg[['children']][['children']][[a]][['children']][[b]][['children']])) {
                    pathway_info <>'children']][['children']][[a]][['children']][[b]][['name']][[c]]

                    pathway_id <>'ko[0-9]{5}')[1]
                    pathway_name <>' \\[PATH:ko[0-9]{5}\\]', '') %>% str_replace('[0-9]{5} ', '')
                    pathway2name <>

                    kos_info <>'children']][['children']][[a]][['children']][[b]][['children']][[c]][['name']]

                    kos <>'K[0-9]*')[,1]

                    ko2pathway <>
                }
            }
        }

        save(pathway2name, ko2pathway, file = 'kegg_info.RData')
    }

    update_kegg(json = 'ko00001.json')

}

STEP5: 利用GO将gene与pathway联系起来，然后挑出query_name与pathway注释信息

load(file = 'kegg_info.RData')
gene2pathway <- gene2ko %>% left_join(ko2pathway, by = 'KO') %>% 
        dplyr::select(GID, Pathway) %>%
        na.omit()

STEP6： 制作自己的Orgdb

# 查询物种的Taxonomy，例如要查sesame
# https://www.ncbi.nlm./taxonomy/?term=sesame
tax_id = '4182'
genus = 'Sesamum' 
species = 'indicum'

makeOrgPackage(gene_info=gene_info,
                   go=gene2go,
                   ko=gene2ko,
                   pathway=gene2pathway,
                   version='0.0.1',
                   outputDir = '.',
                   tax_id=tax_id,
                   genus=genus,
                   species=species,
                   goTable='go')
sesame.orgdb <>'org.', str_to_upper(str_sub(genus, 1, 1)) , species, '.eg.db', sep = '')

有了Orgdb就可以做富集分析了

# enrichGO最主要的目的就是将基因编号转换成GO号
ego <>
                #模式物种
                #OrgDb = org.Mm.eg.db,
                #非模式物种，例如芝麻
                OrgDb = sesame.orgdb,
                ont = 'BP', #或MF或CC
                pAdjustMethod = 'BH',
                #pvalueCutoff  = 0.01,
                qvalueCutoff  = 0.01) 
# 同理也能做enrichKEGG