中文字幕理论片,69视频免费在线观看,亚洲成人app,国产1级毛片,刘涛最大尺度戏视频,欧美亚洲美女视频,2021韩国美女仙女屋vip视频

前言

關(guān)于clusterProfiler這個(gè)R包就不介紹了，網(wǎng)紅教授宣傳得很成功，功能也比較強(qiáng)大，主要是做GO和KEGG的功能富集及其可視化。簡(jiǎn)單總結(jié)下用法，以后用時(shí)可直接找來(lái)用。

首先考慮一個(gè)問(wèn)題：clusterProfiler做GO和KEGG富集分析的注釋信息來(lái)自哪里？

GO的注釋信息來(lái)自Bioconductor，提供了19個(gè)物種的org類型的GO注釋信息，如下表所示。Bioconductor中更多的注釋包可參考http://www.bioconductor.org/packages/release/data/annotation/，很亂，大多數(shù)我都不知道干啥用的。

KEGG的注釋信息clusterProfiler通過(guò)KEGG 數(shù)據(jù)庫(kù)的API來(lái)獲取，https://www.kegg.jp/kegg/rest/keggapi.html。

首先是一個(gè)物種所有基因?qū)?yīng)的pathway注釋文件，比如人的：http://rest.kegg.jp/link/hsa/pathway。
其次還需要pathway對(duì)應(yīng)的描述信息，比如人的：
http://rest.kegg.jp/list/pathway/hsa。

關(guān)于KEGG數(shù)據(jù)庫(kù)全部的物種及其簡(jiǎn)寫（三個(gè)字母）如下列表：
https://www.genome.jp/kegg/catalog/org_list.html。

因此對(duì)于以上已有pathway注釋的物種，只需要將物種簡(jiǎn)寫輸入給clusterProfiler， 它會(huì)通過(guò)聯(lián)網(wǎng)自動(dòng)獲取該物種的pathway注釋信息。

以上都是有物種信息的情況，那么對(duì)于無(wú)物種信息的項(xiàng)目怎么辦？

GO可以通過(guò)讀取外部的GO注釋文件進(jìn)行分析。關(guān)于基因的GO注釋，interproscan、eggnog-mapper和blas2go等軟件都可以做，不過(guò)輸出格式有些不同。clusterProfiler需要導(dǎo)入的GO注釋文件的格式如下：

clusterProfiler包只針對(duì)含有OrgDb對(duì)象，如果是公共數(shù)據(jù)庫(kù)中有該物種注釋信息，只是未制作成org.db數(shù)據(jù)庫(kù)（標(biāo)準(zhǔn)注釋庫(kù)），則可以不需要從頭注釋，只需手動(dòng)制作org.db數(shù)據(jù)庫(kù)類型，完成后直接使用即可，代碼如下：

source("https://bioconductor.org/biocLite.R")
if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")

BiocManager::install("AnnotationHub") # 一個(gè)包含大量注釋信息的數(shù)據(jù)庫(kù)，里面有很多物種及來(lái)源于很多數(shù)據(jù)庫(kù)的注釋信息。
BiocManager::install("biomaRt")

library(AnnotationHub) 
library(biomaRt)

hub <- AnnotationHub() #建立AnnotationHub對(duì)象（視人品，網(wǎng)不行加載不了）
# unique(hub$species) #查看AnonotationHub里面物種
hub$species[which(hub$species=="Solanum")] #看AnonotationHub里是否包含想要的物種
# Solanum是番茄的拉丁名
query(hub, "Solanum")  #查看該物種信息
hub[hub$species=="Solanum" & hub$rdataclass == "OrgDb"] #OrgDb屬于rdataclass中，因此查看下該物種有沒(méi)有OrgDb
Solanum.OrgDb <- hub[["AH59087"]]#AH59087是番茄對(duì)應(yīng)的編號(hào)
#制作為標(biāo)準(zhǔn)注釋庫(kù)，就可和模式生物一樣使用了

同樣地，對(duì)于pathway數(shù)據(jù)庫(kù)中沒(méi)有的物種，也支持讀取基因的pathway注釋文件，然后進(jìn)行分析，注釋文件的格式如下：

富集分析

通常用的富集分析有ORA、FCS和拓?fù)淙N方法。ORA簡(jiǎn)單來(lái)說(shuō)就是超幾何檢驗(yàn)或Fisher精確檢驗(yàn)，大同小異，都符合超幾何檢驗(yàn)，這也是目前用的最多的方法，優(yōu)劣不談。FCS的代表就是GSEA，即基因集富集分析，優(yōu)劣亦不談。clusterProfiler提供了這兩種富集分析方法。
1. ORA（Over-Representation Analysis）
GO富集參考代碼：

#標(biāo)準(zhǔn)富集分析
ego <- enrichGO(
          gene  = gene$entrzID,
          keyType = "ENTREZID", 
          universe = names(geneList), #背景基因集，可省
          OrgDb   = org.Hs.eg.db,
          ont     = "CC",
          pAdjustMethod = "BH",
          pvalueCutoff  = 0.01,
          qvalueCutoff  = 0.05,
          readable      = TRUE)

#通過(guò)導(dǎo)入外部注釋文件富集分析
data <- read.table("go_annotation.txt",header = T,sep = "\t")
go2gene <- data[, c(2, 1)]
go2name <- data[, c(2, 3)]
x <- enricher(gene,TERM2GENE = go2gene,TERM2NAME = go2name)

gene差異基因?qū)?yīng)的向量;
keyType指定基因ID的類型，默認(rèn)為ENTREZID, 可參考keytypes(org.Hs.eg.db)類型 ；
OrgDb指定該物種對(duì)應(yīng)的org包的名字；
ont代表GO的3大類別，BP, CC, MF，也可是全部ALL；
pAdjustMethod指定多重假設(shè)檢驗(yàn)矯正的方法，有“ holm”, “hochberg”, “hommel”, “bonferroni”, “BH”, “BY”, “fdr”, “none”中的一種；
cufoff指定對(duì)應(yīng)的閾值；
readable=TRUE代表將基因ID轉(zhuǎn)換為gene symbol。

KEGG Pathway富集參考代碼：

#標(biāo)準(zhǔn)富集分析
ego <- enrichKEGG(
          gene = gene,
          keyType = "kegg",
          organism  = 'hsa',
          pvalueCutoff  = 0.05,
          pAdjustMethod  = "BH",
          qvalueCutoff  = 0.05
)


#通過(guò)外部導(dǎo)入注釋文件富集
data <- read.table("pathway_annotation.txt",header = T,sep = "\t")
go2gene <- data[, c(2, 1)]
go2name <- data[, c(2, 3)]
x <- enricher(gene,TERM2GENE = go2gene,TERM2NAME = go2name)

默認(rèn)基因ID為kegg gene id，也可以是ncbi-geneid, ncbi-proteinid, uniprot等。
organism物種對(duì)應(yīng)的三字母縮寫，其他參數(shù)同GO富集。ID轉(zhuǎn)換函數(shù)：

library(clusterProfiler)
bitr_kegg("1",fromType = "kegg",toType = 'ncbi-proteinid',organism='hsa')

library(org.Hs.eg.db)
keytypes(org.Hs.eg.db) #支持的ID類型
bitr(gene, fromType = "ENTREZID", toType = c("ENSEMBL", "SYMBOL"), OrgDb = org.Hs.eg.db)

#以上看出ID轉(zhuǎn)換輸入時(shí)，可以向量的形式，也可以單列基因名list導(dǎo)入，也可以是內(nèi)置數(shù)據(jù)
gene <- c("AASDH","ABCB11","ADAM12","ADAMTS16","ADAMTS18")
gene  <-  data$V1 #字符串

data(geneList, package="DOSE") #富集分析的背景基因集
gene <- names(geneList)[abs(geneList) > 2]

2. GSEA（Gene Set Enrichment Analysis）
GO富集參考代碼：

#標(biāo)準(zhǔn)富集分析
ego <- gseGO(
      geneList  = geneList,
      OrgDb  = org.Hs.eg.db,
      ont  = "CC",
      nPerm  = 1000,  #置換檢驗(yàn)的置換次數(shù)
      minGSSize  = 100,
      maxGSSize  = 500,
      pvalueCutoff = 0.05,
      verbose  = FALSE)

#通過(guò)導(dǎo)入外部注釋文件富集分析參考代碼：
data <- read.table("go_annotation.txt",header = T,sep = "\t")
go2gene <- data[, c(2, 1)]
go2name <- data[, c(2, 3)]
x <- GSEA(gene,TERM2GENE = go2gene,TERM2NAME = go2name)

KEGG Pathway富集參考代碼：

#標(biāo)準(zhǔn)富集分析
kk <- gseKEGG(
  geneList  = gene,
  keyType  = 'kegg',
  organism = 'hsa',
  nPerm  = 1000,
  minGSSize = 10,
  maxGSSize = 500,
  pvalueCutoff = 0.05,
  pAdjustMethod     = "BH"
)

#通過(guò)外部導(dǎo)入注釋文件富集
data <- read.table("pathway_annotation.txt",header = T,sep = "\t")
go2gene <- data[, c(2, 1)]
go2name <- data[, c(2, 3)]
x <- GSEA(gene,TERM2GENE = go2gene,TERM2NAME = go2name)

可視化

1.GO富集分析結(jié)果可視化

#barplot
barplot(ego, showCategory = 10) #默認(rèn)展示顯著富集的top10個(gè)，即p.adjust最小的10個(gè)

#dotplot
dotplot(ego, showCategory = 10)

#DAG有向無(wú)環(huán)圖
plotGOgraph(ego)  #矩形代表富集到的top10個(gè)GO terms, 顏色從黃色過(guò)濾到紅色，對(duì)應(yīng)p值從大到小。

#igraph布局的DAG
goplot(ego)

#GO terms關(guān)系網(wǎng)絡(luò)圖（通過(guò)差異基因關(guān)聯(lián)）
emapplot(ego, showCategory = 30)

#GO term與差異基因關(guān)系網(wǎng)絡(luò)圖
cnetplot(ego, showCategory = 5)

2.Pathway富集分析結(jié)果可視化

#barplot
barplot(kk, showCategory = 10)

#dotplot
dotplot(kk, showCategory = 10)

#pathway關(guān)系網(wǎng)絡(luò)圖（通過(guò)差異基因關(guān)聯(lián)）
emapplot(kk,  showCategory = 30)

#pathway與差異基因關(guān)系網(wǎng)絡(luò)圖
cnetplot(kk, showCategory = 5)

#pathway映射
browseKEGG(kk, "hsa04934") #在pathway通路圖上標(biāo)記富集到的基因，會(huì)鏈接到KEGG官網(wǎng)