免费视频淫片aa毛片_日韩高清在线亚洲专区vr_日韩大片免费观看视频播放_亚洲欧美国产精品完整版

豆豆寫于19.1.20-21

為什么要搞全基因組測(cè)序（一）

測(cè)序了，然后呢（二）基因功能注釋

測(cè)序了，然后呢（三）功能注釋整合流程

要進(jìn)行GO或者KEGG富集分析，就需要知道每個(gè)基因?qū)?yīng)什么樣的GO/KEGG分類，OrgDb就是存儲(chǔ)不同數(shù)據(jù)庫(kù)基因ID之間對(duì)應(yīng)關(guān)系，以及基因與GO等注釋的對(duì)應(yīng)關(guān)系的 R 軟件包

如果自己研究的物種不在http://bioconductor.org/packages/release/BiocViews.html#___OrgDb 之列，很大可能就需要自己構(gòu)建OrgDb，然后用clusterProfiler分析

所以本次的內(nèi)容都是基于非模式生物

--------------------------------------------------------------------------

發(fā)推送的時(shí)候花花和豆豆就在看我們十二期的學(xué)員們?cè)诤?jiǎn)書第七天的感想，看的好感動(dòng)，謝謝你們的支持，看到你們真的沒有白付出，我們很欣慰。
今天其實(shí)很累，但是寫推送讓我重回活力，我愛寫推送，我愛和你們交流
感謝所有關(guān)注我們的小伙伴，我們會(huì)更加努力，招財(cái)貓送你們好運(yùn)！?。?/p>

分類

非模式生物要想找到自己的注釋包，又分成兩類：

• 一類是在AnnotationHub（https://bioconductor.org/packages/release/bioc/html/AnnotationHub.html）中存在的，例如棉鈴蟲

• 另一類是在AnnotationHub也不存在相應(yīng)物種，就需要用AnnotationForge（ https://bioconductor.org/packages/release/bioc/html/AnnotationForge.html ）來自己構(gòu)建

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

下面我以棉鈴蟲為例

首先下載并加載AnnotationHub

然后加載現(xiàn)有的物種database

hub <>#這一步受限于網(wǎng)速，不成功的話多試幾次
# 調(diào)用圖形界面查看物種
display(hub) 
# 或者根據(jù)物種拉丁文名稱查找
query(hub,'helicoverpa')
#   'object[['AH66950']]' 

            title                                      
  AH66950 | org.Helicoverpa_armigera.eg.sqlite         
  AH66951 | org.Heliothis_(Helicoverpa)_armigera.eg....
# 這里AH66950是我們需要的
# 然后下載這個(gè)sqlite數(shù)據(jù)庫(kù)
ha.db <>'AH66950']]
#查看前幾個(gè)基因（Entrez命名）
head(keys(ha.db))
#查看包含的基因數(shù)
length(keys(ha.db)) 
#查看包含多少種ID
columns(ha.db)
#查看前幾個(gè)基因的ID
select(ha.db, keys(ha.db)[1:3], 
       c('REFSEQ', 'SYMBOL'), #想獲取的ID
       'ENTREZID')
#得到結(jié)果
  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')
#之后再使用直接加載進(jìn)來
maize.db <>'Harms-AH66950.sqlite')

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

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

通過上次的推送（功能注釋整合流程），我們知道了怎樣利用eggnog-mapper去人工構(gòu)建注釋、富集橋梁

上次選取的是一個(gè)病毒的小例子，這次可以用芝麻（Sesame）做演示

先下載蛋白序列

因?yàn)榻y(tǒng)計(jì)了下有38406條序列，因此使用diamond來進(jìn)行功能注釋

# 前提是自己安裝好eggnog-mapper并且下載好相應(yīng)的數(shù)據(jù)庫(kù)
emapper.py -m diamond \
           -i sesame.fa \
           -o diamond \
           --cpu 19
# 得到如下信息，然后進(jìn)行處理，只保留表頭query_name這一行的注釋信息，去掉頭尾的# 等信息
sed -i '/^# /d' diamond.emapper.annotations 
sed -i 's/#//' diamond.emapper.annotations 

關(guān)于結(jié)果解釋：https://github.com/jhcepas/eggnog-mapper/wiki/Results-Interpretation

其中關(guān)于COG的介紹：倒數(shù)第二列

COG functional categories: COG functional category inferred from best matching OG 【會(huì)給出一個(gè)大寫字母，每一個(gè)大寫字母都有自己的解釋：COG_explanation】

STEP1：自己構(gòu)建的話，首先需要讀入文件

STEP2: 從文件中挑出基因query_name與eggnog注釋信息

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

STEP3-1：挑出query_name與GO注釋信息

STEP3-2：我們想得到query_name與GO號(hào)的對(duì)應(yīng)信息

# 先構(gòu)建一個(gè)空的數(shù)據(jù)框(弄好大體的架構(gòu)，表示其中要有GID =》query_name，GO =》GO號(hào)， EVIDENCE =》默認(rèn)IDA)
# 關(guān)于IEA：就是一個(gè)標(biāo)準(zhǔn)，除了這個(gè)標(biāo)準(zhǔn)以外還有許多。IEA就是表示我們的注釋是自動(dòng)注釋，無(wú)需人工檢查http://wiki.geneontology.org/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對(duì)應(yīng)多個(gè)GO，因此我們以GO號(hào)為標(biāo)準(zhǔn)，每一行只能有一個(gè)GO號(hào)，但query_name和Evidence可以重復(fù)
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注釋信息

STEP4-2: 得到pathway2name, ko2pathway

這一步不需要管代碼什么意思，只需要知道它可以幫我們得到以上兩個(gè)文件就好

if(F){
    # 需要下載 json文件(這是是經(jīng)常更新的)
    # https://www.genome.jp/kegg-bin/get_htext?ko00001
    # 代碼來自：http://www.genek.tv/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聯(lián)系起來，然后挑出query_name與pathway注釋信息

STEP6： 制作自己的Orgdb

# 查詢物種的Taxonomy，例如要查sesame
# https://www.ncbi.nlm.nih.gov/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就可以做富集分析了

剩下的可以參考：http://bioconductor.org/packages/release/bioc/vignettes/clusterProfiler/inst/doc/clusterProfiler.html#supported-organisms

最后，如果不想構(gòu)建orgdb還想做富集分析

滿足需求永遠(yuǎn)是第一生產(chǎn)力，如果你只想用一次，那么clusterProfiler的enricher可以學(xué)習(xí)一下

主要的兩個(gè)參數(shù)需要注意：gene 是基因ID，TERM2GENE 是GO/KEGG terms與基因ID的對(duì)應(yīng)，例如上面圖片中的GO_terms、KEGG_KOs等eggnog-mapper結(jié)果，提取出來就好。對(duì)于沒有對(duì)應(yīng)terms的基因ID，那我們就把它們?nèi)サ簟?/p>

參考：http://guangchuangyu.github.io/2015/05/use-clusterprofiler-as-an-universal-enrichment-analysis-tool/