对富集结果进行贝叶斯网络推断

安装补充包

1
2
3
4
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’)

准备数据

折腾exp矩阵的基因ID

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
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')

折腾通路的名称大小写

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

推断基因

准备画图

1
2
3
4
5
6
7
8
9
10
11
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")

另一种风格

1
2
3
4
5
6
7
8
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)

推断通路

  • 改回折腾通路的名称大小写中通路的名字
1
2
3
bnpathplot(results = pway, exp = vsted, expSample = incSample, R = 200,
           nCategory = 5,
           expRow='ENSEMBL', orgDb=org.Hs.eg.db)

Reactome

1
2
3
4
5
pway@ontology = 'Reactome'
bnpathplot(results = pway, exp = vsted, expSample = incSample, R = 200,
           nCategory = 5,
           compareRef=T, 
           expRow='ENSEMBL', orgDb=org.Hs.eg.db)
uncategorized
对富集结果进行贝叶斯网络推断
https://b.limour.top/1994.html
Author
Limour
Posted on
August 22, 2022
Licensed under