Nguyen的评估代码复现

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)  #关闭并行核

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)

 #将基因集大小与基因集数量进行拼接