occdn.limour.top/source/_posts/对富集结果进行贝叶斯网络推断.md

---

title: 对富集结果进行贝叶斯网络推断 tags: [] id: '2275' categories:

• - 通路富集 date: 2022-08-22 18:37:09 ---

安装补充包

• conda activate clusterprofiler
• ~/dev/xray/xray -c ~/etc/xui2.json &
• wget -e "https_proxy=http://127.0.0.1:20809" https://github.com/noriakis/CBNplot/archive/refs/heads/main.zip -O CBNplot-master.zip
• conda install -c conda-forge r-rmpfr -y
• # conda install -c bioconda bioconductor-depmap -y

• BiocManager::install("depmap", force = TRUE)

  dep = depmap::depmap_crispr()
  rnai = depmap::depmap_rnai()
  depMeta = depmap::depmap_metadata()
  save(dep, rnai, depMeta, file = '~/upload/zl_liu/gsea/CBNplot_depmap.rdata')
  

• devtools::install_local('CBNplot-master.zip')

准备数据

• 差异基因、有批次的差异基因
• 通路富集Demo
• 表达矩阵标准化，emm，这个包有bug，exp矩阵的基因不能是SYMBOL，还得折腾。

折腾exp矩阵的基因ID

library(SummarizedExperiment)
tmp <- load('../../../DEG/TCGA/PRAD_tp.rda')
counts <- data@assays@data$unstranded
colnames(counts) <- data@colData$patient
f_rm_duplicated <- function(NameL, reverse=F){
    tmp <- data.frame(table(NameL))
    if(reverse){
        tmp <- tmp$NameL[tmp$Freq > 1]
    }else{
        tmp <- tmp$NameL[tmp$Freq == 1]
    }
    which(NameL %in% as.character(tmp))
}
counts <- counts[,f_rm_duplicated(colnames(counts))]
geneInfo <- as.data.frame(data@rowRanges)[c('gene_id','gene_type','gene_name')]
f_dedup_IQR <- function(df, rowNn, select_func='IQR'){
    if(typeof(select_func) == 'character'){
        select_func = get(select_func)
    }
    # 拆出无重复的数据，后续不进行处理
    noDup <- f_rm_duplicated(rowNn)
    tmp <- rowNn[noDup]
    noDup <- df[noDup,]
    rownames(noDup) <- tmp
    # 拆除有重复的数据
    Dup <- f_rm_duplicated(rowNn, T)
    rowNn <- rowNn[Dup]
    Dup <- df[Dup,]
    rownames(Dup) <- NULL
    # 处理重复的数据
    lc_tmp = by(Dup,
         rowNn,
         function(x){rownames(x)[which.max(apply(X = x, FUN = select_func, MARGIN = 1))]})
    lc_probes = as.integer(lc_tmp)
    Dup = Dup[lc_probes,]
    rownames(Dup) <- rowNn[lc_probes]
    # 合并数据并返回
    return(rbind(noDup,Dup))
}
require(stringr)
rininiang <- t(as.data.frame(str_split(geneInfo$gene_id, '\\.')))
rownames(rininiang) <- NULL
rininiang <- f_dedup_IQR(as.data.frame(counts),unlist(rininiang[,1]))
rininiang
f_counts2VST <- function(countsMatrix){
    require(DESeq2)
    conditions <- factor(rep("Control",ncol(countsMatrix)))
    colData_b <- data.frame(row.names = colnames(countsMatrix), conditions)
    dds <- DESeqDataSetFromMatrix(countData = countsMatrix,
                              colData = colData_b,
                              design = ~ 1)
    vsd <- vst(object=dds, blind=T) 
    assay(vsd)
}
VST <- f_counts2VST(rininiang)
VST
saveRDS(VST, 'rininiang.rds')

折腾通路的名称大小写

• ~/dev/xray/xray -c ~/etc/xui2.json
• Sys.setenv(http_proxy = "http://127.0.0.1:20809")

• Sys.setenv(https_proxy = "http://127.0.0.1:20809")

  rininiang <- graphite::pathways('hsapiens','reactome')
  pway@result[1, "Description"]
  rininiang[['Cell Cycle']]
  pway@result[1, "Description"] <- 'Cell Cycle'
  

推断基因

准备画图

require(org.Hs.eg.db)
set.seed(123)
bngeneplot(results = pway, exp = vsted, expSample = incSample, R = 200, cl = cl,
           pathNum = 1,
           strThresh = 0.8,
           showDir = T,
           hub=5, 
#             sizeDep = T, dep = dep, cellLineName = "UMUC3_URINARY_TRACT", showLineage = T, depMeta = depMeta,
           pathDb = "reactome",  compareRef = T, compareRefType = "difference",
           labelSize=7, shadowText=TRUE, 
           expRow='ENSEMBL', convertSymbol=T, orgDb=org.Hs.eg.db, sp = "hsapiens")

另一种风格

set.seed(123)
bngeneplotCustom(results = pway, exp = vsted, expSample = incSample, R = 200, cl = cl,
           pathNum = 1,
           strThresh = 0.8,
           showDir = T,
           hub=5, 
           expRow='ENSEMBL', convertSymbol=T, orgDb=org.Hs.eg.db,
           layout="nicely", fontFamily="sans", strType="normal", glowEdgeNum=5,  labelSize=7)

推断通路

• 改回折腾通路的名称大小写中通路的名字

  bnpathplot(results = pway, exp = vsted, expSample = incSample, R = 200,
         nCategory = 5,
         expRow='ENSEMBL', orgDb=org.Hs.eg.db)
  

Reactome

pway@ontology = 'Reactome'
bnpathplot(results = pway, exp = vsted, expSample = incSample, R = 200,
           nCategory = 5,
           compareRef=T, 
           expRow='ENSEMBL', orgDb=org.Hs.eg.db)