---
title: 对富集结果进行贝叶斯网络推断
tags: []
id: '2275'
categories:
  - - 通路富集
date: 2022-08-22 18:37:09
---

## 安装补充包

*   [conda activate clusterprofiler](https://occdn.limour.top/2126.html)
*   ~/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)

```R
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')

## 准备数据

*   [差异基因](https://occdn.limour.top/2132.html)、[有批次的差异基因](https://occdn.limour.top/2197.html)
*   [通路富集Demo](https://occdn.limour.top/2190.html)
*   [表达矩阵标准化](https://occdn.limour.top/2195.html)，emm，这个包有bug，exp矩阵的基因不能是SYMBOL，还得折腾。

### 折腾exp矩阵的基因ID

```R
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")

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

## 推断基因

### 准备画图

```R
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")
```

### 另一种风格

```R
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)
```

## 推断通路

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

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

### Reactome

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