美文网首页NGS
基因集的转录因子富集分析

基因集的转录因子富集分析

作者: 医学小白学生信 | 来源:发表于2021-04-05 12:27 被阅读0次

    基因集的转录因子富集
    为什么是AUC值而不是GSEA来挑选转录因子呢

    library("RcisTarget")
    # 官网的genelist
    geneList1 <- read.table(file.path(system.file('examples', package='RcisTarget'), "hypoxiaGeneSet.txt"), 
                            stringsAsFactors=FALSE)[,1]
    geneLists <- list(hypoxia=geneList1 ) 
    geneLists
    
    
    #motifAnnotations_hgnc里面包含了16万行,可以看到每个基因有多少个motif
    data(motifAnnotations_hgnc)
    motifAnnotations_hgnc
    
    
    #前面我们下载好的 hg19-tss-centered-10kb-7species.mc9nr.feather 文件,需要存储在当前工作目录文件夹里面:
    #这个文件下载了,在生信技能树的SCENIC文件夹里面
    motifRankings <- importRankings("C:/Users/18700/Desktop/123/hg19-tss-centered-10kb-7species.mc9nr.feather")
    # Motif enrichment analysis:geneSets必须是list
    #最后结果可以看到100个motif被富集到,同时还计算了AUC和nes值。
    motifEnrichmentTable_wGenes <- cisTarget(geneLists, motifRankings,
                                             motifAnnot=motifAnnotations_hgnc)
    
    
    
    #因为通过RcisTarget包的 cisTarget()函数,一句代码完成的转录因子富集分析,等价于下面的3个步骤。
    # 1. Calculate AUC
    motifs_AUC <- calcAUC(geneLists, motifRankings)
    
    # 2. Select significant motifs, add TF annotation & format as table
    motifEnrichmentTable <- addMotifAnnotation(motifs_AUC, 
                                               motifAnnot=motifAnnotations_hgnc)
    # 3. Identify significant genes for each motif
    motifEnrichmentTable_wGenes <- addSignificantGenes(motifEnrichmentTable, 
                                                       geneSets=geneLists,
                                                       rankings=motifRankings, 
                                                       nCores=1,
                                                       method="aprox")
    

    可以看到富集到了100个motif,可以看到每一个motif富集到的 highly ranked 基因

    image.png
    image.png
    #首先批量计算AUC
    motifs_AUC <- calcAUC(geneLists, motifRankings, nCores=1)
    motifs_AUC#查看这个基因集有多少个motif参与了计算
    
    
    #挑选统计学显著的motif
    auc <- getAUC(motifs_AUC)[1,]
    #画图判断auc是否符合正态分布
    hist(auc, main="hypoxia", xlab="AUC histogram",
         breaks=100, col="#ff000050", border="darkred")
    #一般来说,对正态分布,我们会挑选 mean+2sd范围外的认为是统计学显著,
    #但这里选择的是 mean+3sd,更加严格 
    nes3 <- (3*sd(auc)) + mean(auc)
    abline(v=nes3, col="red")
    
    data(motifAnnotations_hgnc) #读入所有的motif注释文件
    motifAnnotations_hgnc
    cg=auc[auc>nes3]#显示符合要求的有多少个motif
    names(cg)
    cgmotif=motifAnnotations_hgnc[match(names(cg),motifAnnotations_hgnc$motif),]
    cgmotif=na.omit(cgmotif)#去除NA,剩余的就是挑选得到的motif
    
    image.png 最后得到82个motif
    #高级分析之可视化motif
    motifEnrichmentTable_wGenes
    #addLogo函数添加可视化图表
    motifEnrichmentTable_wGenes_wLogo <- addLogo(motifEnrichmentTable_wGenes)
    library(DT)
    datatable(motifEnrichmentTable_wGenes_wLogo[,-c("enrichedGenes", "TF_lowConf"), with=FALSE], 
              escape = FALSE, # To show the logo
              filter="top", options=list(pageLength=5))
    
    
    
    ##高级分析之网络图
    #数据库直接注释和同源基因推断的TF为高可信靶基因TF_highConf,
    #使用motif序列相似性注释的TF为低可信度结果。
    #将高可信靶基因拉出来,共51个
    anotatedTfs <- lapply(split(motifEnrichmentTable_wGenes$TF_highConf,
                                motifEnrichmentTable$geneSet),
                          function(x) {
                            genes <- gsub(" \\(.*\\). ", "; ", x, fixed=FALSE)
                            genesSplit <- unique(unlist(strsplit(genes, "; ")))
                            return(genesSplit)
                          })
    
    anotatedTfs$hypoxia
    signifMotifNames <- motifEnrichmentTable$motif[1:4]#取前4个motif绘制网络图
    
    #Identify which genes (of the gene-set) are highly ranked for each motif.
    #incidMatrix里会有enrStatsh和incidMatrix两个,选择incidMatrix
    #incidMatrix是靶基因和哪个TF相关的信息
    incidenceMatrix <- getSignificantGenes(geneLists$hypoxia, 
                                           motifRankings,
                                           signifRankingNames=signifMotifNames,
                                           plotCurve=TRUE, maxRank=5000, 
                                           genesFormat="incidMatrix",
                                           method="aprox")$incidMatrix
    
    #将绘图需要的edge和node文件进行创建
    library(reshape2)
    edges <- melt(incidenceMatrix)
    edges <- edges[which(edges[,3]==1),1:2]
    colnames(edges) <- c("from","to")
    library(visNetwork)
    motifs <- unique(as.character(edges[,1]))
    genes <- unique(as.character(edges[,2]))
    nodes <- data.frame(id=c(motifs, genes),   
                        label=c(motifs, genes),    
                        title=c(motifs, genes), # tooltip 
                        shape=c(rep("diamond", length(motifs)), rep("elypse", length(genes))),
                        color=c(rep("purple", length(motifs)), rep("skyblue", length(genes))))
    visNetwork(nodes, edges) %>% visOptions(highlightNearest = TRUE, 
                                            nodesIdSelection = TRUE)
    
    
    image.png image.png

    相关文章

      网友评论

        本文标题:基因集的转录因子富集分析

        本文链接:https://www.haomeiwen.com/subject/caniiktx.html