一、生信分析的目的
这里使用R语言IOBR包对患者进行免疫浸润分析，其中包括CIBERSORT，mcpcounter，xcell，timer,quantiseq多种算法。
1）estimate
2）TIMER
3)ABIS
4)ConsensusTME
5)xCell
6)EPIC
7)quanTIseq
8)CIBERSORT
1、CIBERSORT是一种用于计算免疫细胞相对丰度的计算工具，可以通过分析基因表达数据估计复杂组织或血液样本中不同免疫细胞类型的相对丰度。根据样本基因表达数据和参考基因签名，使用机器学习算法（如线性混合模型）来估计每种免疫细胞类型的相对丰度。算法通过比较样本中基因表达模式与参考基因签名之间的相似度来推断免疫细胞的相对丰度。根据计算得到的结果，可以获得不同免疫细胞类型的相对丰度信息，例如各类淋巴细胞、单核细胞、重要的调节性T细胞亚群等。这些信息有助于理解免疫细胞在正常生理状态或疾病进程中的作用和变化。
2、MCPcounter（Microenvironment Cell Populations-counter）是一种用于推测肿瘤组织中免疫细胞丰度的计算工具。它利用肿瘤组织的基因表达数据，通过分析特定的基因集合来估计不同免疫细胞亚群的相对丰度。MCPcounter利用特定基因集合的表达模式来推测肿瘤组织中免疫细胞亚群的相对丰度。它通过比较样本的基因表达数据与事先建立的基因表达模式之间的匹配程度来推测免疫细胞的相对丰度。通常，线性模型或其他机器学习算法用于这个推断过程。根据计算结果，可以获得肿瘤组织中不同免疫细胞亚群的相对丰度信息。这些信息有助于研究人员理解免疫细胞在肿瘤微环境中的组成和相互作用，并推测免疫细胞对肿瘤免疫反应的影响。
3、xCell是一种用于计算复杂组织中免疫细胞相对丰度的计算工具。它通过分析基因表达数据来估计不同免疫细胞类型在组织中的相对含量，包括肿瘤组织、正常组织等。根据以往的研究，建立不同免疫细胞类型的参考字典。参考字典是由每个免疫细胞类型特有的基因组成，这些基因在该细胞类型中显示高表达。通过使用参考字典和基因表达数据，xCell计算每个细胞类型的相对丰度。它种通过比较样本中基因表达数据与参考字典之间的相似度来推断免疫细胞的相对丰度。根据计算结果，可以获得不同免疫细胞类型的相对丰度信息，例如淋巴细胞、单核细胞、造血细胞等。这些信息可用于理解免疫细胞在组织中的分布和相互作用，并推断它们对特定生理或病理状态的影响。
4、Timer（Tumor IMmune Estimation Resource）是一种被广泛应用的计算工具，用于计算肿瘤组织中免疫细胞的相对丰度。它通过分析基因表达数据来估计肿瘤微环境中不同免疫细胞类型的存在程度。使用基因表达数据和参考数据集，Timer通过比较样本中标记基因表达的模式，推断不同免疫细胞类型的相对丰度。根据Timer的计算结果，可以得到不同免疫细胞类型在肿瘤组织中的相对丰度信息。这些信息可用于研究人员的免疫细胞分析和相关生物学研究，以理解免疫系统在肿瘤发展和治疗中的角色。
5、Quantiseq（Quantitative Sequencing）是一种用于计算免疫细胞丰度的方法，它通过基因表达数据的定量测量来推断不同免疫细胞类型在样本中的存在程度。
使用基因表达数据和基因签名，Quantiseq通过量化样本中基因签名的表达水平来推断不同免疫细胞类型的相对丰度。这一计算过程通常使用统计方法和数学模型来完成。
根据Quantiseq的计算结果，可以得到不同免疫细胞类型在样本中的相对存在程度。这些结果有助于研究人员了解免疫细胞组成和免疫系统在疾病进程中的作用。
二、具体的代码
remove(list = ls()) ##清空当前环境
setwd("D:/生信代码复现/105.多算法免疫浸润") ##💚💚💚💚💚设置路径
library(immunedeconv)
library(job)
library(tcltk)
library(tidyverse)
library(ggsci)
library(scales)
library(ggtext)
library(reshape2)
library(limma)
library(parallel)
library(ggplot2)
library(ggpubr)
library(future.apply)
library(car)
library(ridge)
library(e1071)
library(preprocessCore)
library(foreach)
library(doParallel)
library(data.table)
#生成随机颜色
randomColor <- function() {
paste0("#",paste0(sample(c(0:9, letters[1:6]), 6, replace = TRUE),collapse = ""))
}
# 生成200个随机颜色
randomColors <- replicate(200,randomColor())
# 读入表达数据
exprMatrix <- read.table("TCGA.txt",header=TRUE,sep = "\t",as.is = T,row.names = 1) #💛💛💛本地需要有TCGA.txt这个文件
exprMatrix= exprMatrix%>% dplyr::select(str_which(colnames(.), ".01"))%>%as.data.frame()
cancer="BRCA" #TCGA癌症缩写
brca_vector <- rep(cancer, ncol(exprMatrix))
{
job::job({exp=exprMatrix
res_estimate <- deconvolute(exp, method="estimate",tumor = T,indications = brca_vector)
write.table(res_estimate,"estimate.txt",sep="\t",row.names = F) #💜💜💜输出文件estimate.txt
res_estimate$cell_type=paste0(res_estimate$cell_type,"_ESTIMATE")
})
job::job({exp=exprMatrix
res_timer <- deconvolute(exp, method="timer",tumor = T,indications = brca_vector)
write.table(res_timer,"timer.txt",sep="\t",row.names = F) #💜💜💜输出文件timer.txt
res_timer$cell_type=paste0(res_timer$cell_type,"_TIMER")
})
job::job({exp=exprMatrix
res_abis <- deconvolute(exp, method="abis",tumor = T,indications = brca_vector)
write.table(res_abis,"abis.txt",sep="\t",row.names = F) #💜💜💜输出文件abis.txt
res_abis$cell_type=paste0(res_abis$cell_type,"_ABIS")
})
job::job({exp=exprMatrix
res_consensus_tme <- deconvolute(exp, method="consensus_tme",tumor = T,indications = brca_vector)
write.table(res_consensus_tme,"consensus_tme.txt",sep="\t",row.names = F) #💜💜💜输出文件consensus_tme.txt
res_consensus_tme$cell_type=paste0(res_consensus_tme$cell_type,"_ConsensusTME")
})
job::job({exp=exprMatrix
res_xcell <- deconvolute(exp, method="xcell",tumor = T,indications = brca_vector)
write.table(res_xcell,"xcell.txt",sep="\t",row.names = F) #💜💜💜输出文件xcell.txt
res_xcell$cell_type=paste0(res_xcell$cell_type,"_xCell")
})
job::job({exp=exprMatrix
res_epic <- deconvolute(exp, method="epic",tumor = T,indications = brca_vector)
write.table(res_epic,"epic.txt",sep="\t",row.names = F) #💜💜💜输出文件epic.txt
res_epic$cell_type=paste0(res_epic$cell_type,"_EPIC")
})
job::job({exp=exprMatrix
res_quantiseq <- deconvolute(exp, method="quantiseq",tumor = T,indications = brca_vector)
write.table(res_quantiseq,"quantiseq.txt",sep="\t",row.names = F) #💜💜💜输出文件quantiseq.txt
res_quantiseq$cell_type=paste0(res_quantiseq$cell_type,"_quanTIseq")
})
job::job({data=exprMatrix
v=voom(data, plot=F, save.plot=F)
out=v$E
out=rbind(ID=colnames(out), out)
write.table(out,file="uniq.symbol.txt",sep="\t",quote=F,col.names=F)
source("CIBERSORT多线程主代码.R",encoding = "utf-8") #💛💛💛确保本地存在这个代码CIBERSORT多线程主代码.R
results=CIBERSORT("LM22.txt", "uniq.symbol.txt", perm=1000, QN=TRUE) #💜💜💜输出文件LM22.txt和uniq.symbol.txt
res_CIBERSORT=read.table("CIBERSORT-Results.txt",sep="\t",header=T,row.names=1,check.names=F) #💛💛💛CIBERSORT-Results.txt这个文件也要有
res_CIBERSORT=t(res_CIBERSORT)
rownames(res_CIBERSORT)=paste0(rownames(res_CIBERSORT),"_CIBERSORT")
res_CIBERSORT=data.frame(cell_type=rownames(res_CIBERSORT),res_CIBERSORT)})
}
res_consensus_tme=res_consensus_tme[1:(nrow(res_consensus_tme)-1),]
res_xcell=res_xcell[1:(nrow(res_xcell)-3),]
res_CIBERSORT=res_CIBERSORT[1:(nrow(res_CIBERSORT)-3),]
res=rbind(res_abis,res_consensus_tme,res_epic,res_quantiseq,res_timer,res_xcell,res_CIBERSORT)
res=as.matrix(res)
rownames(res)=res[,1]
exp=res[,2:ncol(res)]
dimnames=list(rownames(exp),colnames(exp))
data=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
#输入要研究的基因
gene="TP53" #💚💚💚这里改成我们想要研究的目的基因的名称
gene_exp=exprMatrix[gene,]
gene_exp=t(gene_exp)
samesample=intersect(rownames(gene_exp),colnames(data))
gene_exp=gene_exp[samesample,]
gene_exp=as.data.frame(gene_exp)
colnames(gene_exp)=gene
data=data[,samesample]
x=as.numeric(gene_exp[,1])
outTab=data.frame()
for(i in rownames(data)){
y=as.numeric(data[i,])
if(sd(y)<0.001){next}
corT=cor.test(x, y, method="spearman")
cor=corT$estimate
pvalue=corT$p.value
if(pvalue<0.05){
outTab=rbind(outTab,cbind(immune=i, cor, pvalue))
}
}
#输出相关性结果
write.table(file="corResult.txt", outTab, sep="\t", quote=F, row.names=F) #💜💜💜输出文件corResult.txt
#绘制气泡图
corResult=read.table("corResult.txt", head=T, sep="\t")
corResult$Software=sapply(strsplit(corResult[,1],"_"), '[', 2)
corResult$Software=factor(corResult$Software,level=as.character(unique(corResult$Software[rev(order(as.character(corResult$Software)))])))
b=corResult[order(corResult$Software),]
b$immune=factor(b$immune,levels=rev(as.character(b$immune)))
colslabels=rep(hue_pal()(length(levels(b$Software))),table(b$Software)) #定义图形颜色
#输出图形
pdf(file="correlation.pdf", width=9, height=10) #💜💜💜输出文件correlation.pdf
ggplot(data=b, aes(x=cor, y=immune, color=Software))+
labs(x="Correlation coefficient",y="Immune cell")+
geom_point(size=4.1)+
theme(panel.background=element_rect(fill="white",size=1,color="black"),
panel.grid=element_line(color="grey75",size=0.5),
axis.ticks = element_line(size=0.5),
axis.text.y = ggtext::element_markdown(colour=rev(colslabels)))
dev.off()
R语言多线程代码：
CoreAlg <- function(X, y){
svn_itor <- 3
res <- function(i){
if(i==1){nus <- 0.25}
if(i==2){nus <- 0.5}
if(i==3){nus <- 0.75}
model<-svm(X,y,type="nu-regression",kernel="linear",nu=nus,scale=F)
model
}
if(Sys.info()['sysname'] == 'Windows') out <- mclapply(1:svn_itor, res, mc.cores=1) else
out <- mclapply(1:svn_itor, res, mc.cores=svn_itor)
nusvm <- rep(0,svn_itor)
corrv <- rep(0,svn_itor)
t <- 1
while(t <= svn_itor) {
weights = t(out[[t]]$coefs) %*% out[[t]]$SV
weights[which(weights<0)]<-0
w<-weights/sum(weights)
u <- sweep(X,MARGIN=2,w,'*')
k <- apply(u, 1, sum)
nusvm[t] <- sqrt((mean((k - y)^2)))
corrv[t] <- cor(k, y)
t <- t + 1
}
rmses <- nusvm
mn <- which.min(rmses)
model <- out[[mn]]
q <- t(model$coefs) %*% model$SV
q[which(q<0)]<-0
w <- (q/sum(q))
mix_rmse <- rmses[mn]
mix_r <- corrv[mn]
newList <- list("w" = w, "mix_rmse" = mix_rmse, "mix_r" = mix_r)
}
doPerm <- function(perm, X, Y){
message(paste0("进行doPerm循环,调用线程数： ",detectCores()))
itor <- 1
Ylist <- as.list(data.matrix(Y))
dist<-list()
registerDoParallel(round(detectCores()))
dist<-foreach(itor = 1:perm) %dopar% {
library(limma)
library(parallel)
library(ggplot2)
library(ggpubr)
library(future.apply)
library(car)
library(ridge)
library(e1071)
library(preprocessCore)
library(tcltk)
library(limma)
library(parallel)
library(ggplot2)
library(ggpubr)
library(future.apply)
library(car)
library(ridge)
library(e1071)
library(preprocessCore)
library(foreach)
library(doParallel)
library(limma)
library(parallel)
library(ggplot2)
library(ggpubr)
library(future.apply)
library(car)
library(ridge)
library(e1071)
library(preprocessCore)
library(tcltk)
library(limma)
library(parallel)
library(ggplot2)
library(ggpubr)
library(future.apply)
library(car)
library(ridge)
library(e1071)
library(preprocessCore)
library(foreach)
library(doParallel)
library(data.table)
CoreAlg <- function(X, y){
svn_itor <- 3
res <- function(i){
if(i==1){nus <- 0.25}
if(i==2){nus <- 0.5}
if(i==3){nus <- 0.75}
model<-svm(X,y,type="nu-regression",kernel="linear",nu=nus,scale=F)
model
}
if(Sys.info()['sysname'] == 'Windows') out <- mclapply(1:svn_itor, res, mc.cores=1) else
out <- mclapply(1:svn_itor, res, mc.cores=svn_itor)
nusvm <- rep(0,svn_itor)
corrv <- rep(0,svn_itor)
t <- 1
while(t <= svn_itor) {
weights = t(out[[t]]$coefs) %*% out[[t]]$SV
weights[which(weights<0)]<-0
w<-weights/sum(weights)
u <- sweep(X,MARGIN=2,w,'*')
k <- apply(u, 1, sum)
nusvm[t] <- sqrt((mean((k - y)^2)))
corrv[t] <- cor(k, y)
t <- t + 1
}
rmses <- nusvm
mn <- which.min(rmses)
model <- out[[mn]]
q <- t(model$coefs) %*% model$SV
q[which(q<0)]<-0
w <- (q/sum(q))
mix_rmse <- rmses[mn]
mix_r <- corrv[mn]
newList <- list("w" = w, "mix_rmse" = mix_rmse, "mix_r" = mix_r)
}
yr <- as.numeric(Ylist[sample(length(Ylist),dim(X)[1])])
yr <- (yr - mean(yr)) / sd(yr)
result <- CoreAlg(X, yr)
mix_r <- result$mix_r
dist[[itor]]<-mix_r
return(dist[[itor]])
}
stopImplicitCluster()
dist<-do.call(rbind,dist)
newList <- list("dist" = dist)
}
CIBERSORT <- function(sig_matrix, mixture_file, perm=0, QN=TRUE){
library(e1071)
library(parallel)
library(preprocessCore)
message("读取表达矩阵")
X <- fread(sig_matrix,header=T,sep="\t",check.names=F,data.table=FALSE)
rownames(X)=X[,1]
X=X[,2:ncol(X)]
Y <- fread(mixture_file, header=T, sep="\t",check.names=F,data.table=FALSE)
rownames(Y)=Y[,1]
Y=Y[,2:ncol(Y)]
X <- na.omit(X)
Y <- na.omit(Y)
X <- data.matrix(X)
Y <- data.matrix(Y)
X <- X[order(rownames(X)),]
Y <- Y[order(rownames(Y)),]
P <- perm
if(max(Y) < 50) {Y <- 2^Y}
if(QN == TRUE){
tmpc <- colnames(Y)
tmpr <- rownames(Y)
Y <- normalize.quantiles(Y)
colnames(Y) <- tmpc
rownames(Y) <- tmpr
}
Xgns <- row.names(X)
Ygns <- row.names(Y)
YintX <- Ygns %in% Xgns
Y <- Y[YintX,]
XintY <- Xgns %in% row.names(Y)
X <- X[XintY,]
X <- (X - mean(X)) / sd(as.vector(X))
if(P > 0) {nulldist <- sort(doPerm(P, X, Y)$dist)}
header <- c('Mixture',colnames(X),"P-value","Correlation","RMSE")
itor <- 1
mixtures <- dim(Y)[2]
pval <- 9999
message(paste0("进行CIBERSORT循环,调用线程数： ",detectCores()))
output=list()
registerDoParallel(round(detectCores()))
output<- foreach(itor = 1:mixtures) %dopar% {
library(limma)
library(parallel)
library(ggplot2)
library(ggpubr)
library(future.apply)
library(car)
library(ridge)
library(e1071)
library(preprocessCore)
library(tcltk)
library(limma)
library(parallel)
library(ggplot2)
library(ggpubr)
library(future.apply)
library(car)
library(ridge)
library(e1071)
library(preprocessCore)
library(foreach)
library(doParallel)
CoreAlg <- function(X, y){
svn_itor <- 3
res <- function(i){
if(i==1){nus <- 0.25}
if(i==2){nus <- 0.5}
if(i==3){nus <- 0.75}
model<-svm(X,y,type="nu-regression",kernel="linear",nu=nus,scale=F)
model
}
if(Sys.info()['sysname'] == 'Windows') out <- mclapply(1:svn_itor, res, mc.cores=1) else
out <- mclapply(1:svn_itor, res, mc.cores=svn_itor)
nusvm <- rep(0,svn_itor)
corrv <- rep(0,svn_itor)
t <- 1
while(t <= svn_itor) {
weights = t(out[[t]]$coefs) %*% out[[t]]$SV
weights[which(weights<0)]<-0
w<-weights/sum(weights)
u <- sweep(X,MARGIN=2,w,'*')
k <- apply(u, 1, sum)
nusvm[t] <- sqrt((mean((k - y)^2)))
corrv[t] <- cor(k, y)
t <- t + 1
}
rmses <- nusvm
mn <- which.min(rmses)
model <- out[[mn]]
q <- t(model$coefs) %*% model$SV
q[which(q<0)]<-0
w <- (q/sum(q))
mix_rmse <- rmses[mn]
mix_r <- corrv[mn]
newList <- list("w" = w, "mix_rmse" = mix_rmse, "mix_r" = mix_r)
}
y <- Y[,itor]
y <- (y - mean(y)) / sd(y)
result <- CoreAlg(X, y)
w <- result$w
mix_r <- result$mix_r
mix_rmse <- result$mix_rmse
if(P > 0) {pval <- 1 - (which.min(abs(nulldist - mix_r)) / length(nulldist))}
out <- c(colnames(Y)[itor],w,pval,mix_r,mix_rmse)
output[[itor]]<-out
return(output[[itor]])
}
stopImplicitCluster()
output<-do.call(cbind,output)
output<-t(output)
output<-rbind(header,output)
write.table(output, file="CIBERSORT-Results.txt", sep="\t", row.names=F, col.names=F, quote=F)
message("运行结束，结果导出")
}
三、最终的效果
①R语言中的效果：
②文件效果：
