这一部分是笔者毕业论文中的对Nguyen的基因集分析方法评估的R语言复现。
前面的准备过程可见其他的帖子。
此处以GSVA方法进行示例：
代码语言：javascript复制library(parallel)
library(doParallel)
library(foreach) #并行计算所使用的R包
library(tibble)
library(dplyr) #包含对dataframe数据的函数
load("/Users/narcissus/Desktop/bioinfomatics/Data/HealthPool_20824.RData")
#文件目录需要根据实际情况进行更改
# gsva的方法评估 - - - - - - - - - - >
NormalPool<-HealthData0
# 这一步根据具体RData文件中的名字进行更改，笔者的文件中的变量为“HealthData0”
source("/Users/narcissus/Desktop/bioinfomatics/gsaMethod/gsva_method.R")
# 引入具体的评估方法
cl <- makeCluster(4)
registerDoParallel(cl) # 开启并行，申请4核，根据电脑实际情况进行更改
GSVA_output<-foreach (1:2000,.combine =cbind,.packages ="GSVA") %dopar%{
sampleIndex=sample(dim(NormalPool)[2],size=30) # 首先采样30个样本
exprSet<-NormalPool[c(sampleIndex)]
# 对于每一个KEGG基因集，使用具体方法得到该基因集的p值列表（总共进行2000次）
GSVA_adjust(exprSet)
} # 每一次循环都会产生一个列向量(基因集数量*1)，通过foreach中的cbind参数将2000个列向量拼接在一起。
save(GSVA_output,file="/Users/narcissus/Desktop/bioinfomatics/Data/method_assess/GSVA_output5367.RData")
#保存结果
#rm(GSVA_output) #可将该变量删除节省空间
#gc() #整理一下内存空间
stopCluster(cl) #关闭并行核得到所有基因集的p值分布数据之后，进行偏性的评估。
使用样本偏度（skewness）进行某一对象（方法-基因集对）偏向0或者偏向1的判断
代码语言：javascript复制library(ggplot2)#进行图像绘制
# 数据加载 - - - - - - - - - - >
load("/Users/narcissus/Desktop/bioinfomatics/Data/method_assess/sigPathway_output20824_.Rdata") #根据实际情况进行更改
P_distribution<-GSVA_output
# 进行偏向0/1判断 - - - - - - - - - - >
SetNum<-dim(P_distribution)[1] #获得基因集数量
SkewValue<-array(0,SetNum) # 以0初始化容器
for (k in 1:SetNum){
mean<-mean(P_distribution[k,]) # 计算均值
sd<-sd(P_distribution[k,]) # 计算标准差
SkewValue[k]<-mean(( ((P_distribution[k,]-mean))/sd )^3)
} # SkewValue可能会有NA的情况，对于这种情况，人工设置为0.5（当该基因集的所有p值相同时）
Prefer0<-length(which(SkewValue>0.1)) #计算偏向0的数量
Prefer1<-length(which(SkewValue< -0.1)) #计算偏向1的数量
PreferNeither<- length(which( abs(SkewValue) <0.1)) # 计算无偏通路的数量
load("/Users/narcissus/Desktop/bioinfomatics/Data/KeggGeneSet.RData") #获得之前得到的KEGG基因集数据
gsSize<-array(0,SetNum) # 初始化容器
for (k in 1:SetNum){
gsSize[k]<-length(pathways[[k]])
}
skew_with_size<-cbind(gsSize,SkewValue) #将基因集大小与基因集数量进行拼接 
