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

安装补充包

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)

对富集结果进行贝叶斯网络推断
https://b.limour.top/1994.html
Author
Limour
Posted on
August 22, 2022
Licensed under