devbox
繁依Fanyi0
这个屌丝很懒，什么也没留下！
关注作者
热门标签
热门文章
当前位置:   article > 正文

手把手教你R语言CIBERSORT计算免疫浸润+Rproject的使用

作者：繁依Fanyi0 | 2024-03-12 04:10:18

赞

踩

cibersort

写在开头：

本文介绍了CIBERSORT两种使用方法，大家可以自行选择，方法二简单些，方法一原始些

本文顺便倡议大家使用Rproject来管理代码，感谢生信技能树jimmy老师让我知道了这么方便的玩意，再也不用拼命setwd()和getwd()了，不想看这部分可以直接下滑。

CIBERSORTx是原版网站，建议大家去学习，并且学习他们发的经典文章

鸣谢：生信技能树jimmy老师和 Biomamba 生信基地 BIOMAMBA老师

使用Rproject管理R项目

1.提前在你想要储存代码的地方建一个文件夹，然后打开Rstudio 中选择 NEW Project

2.选择Existing Directory，因为刚刚已经提前建好了文件夹。如果刚刚没有建，选NEW Directory就会给你建一个，属于个人习惯，都可以。

 3.建立Project，比如我在zhangming这个文件夹下建立，creat后即可建立

 4.在zhangming文件夹中出现了相关的project文件

 5.双击这个文件，就会直接定位在这个工作路径，而不需要切换setwd()

 6.再新建script,新建的code都会在同一个工作环境/工作文件夹内，非常方便，尤其是对于大量的代码学习，可以很好的节省时间。当然如果你不喜欢这样，也可以直接进入下面的CIBERSORT学习。具体不懂的地方还可以看其他人对project的解释，亲测很好上手。

下面进入CIBERSORT的学习

首先安装包

  1. # install packages 这三个安装不成功的话，就安后面的bseqsc包也行
  2. install.packages('e1071')
  3. install.pacakges('parallel')
  4. install.packages('preprocessCore')
  5. library(e1071)
  6. library(preprocessCore)
  7. library(parallel)
  8.  
  9. install.packages('devtools')
  10. library(devtools)
  11. devtools::install_github('shenorrlab/bseqsc')
  12. library(bseqsc)#这个包携带大量CIBERSORT的依赖，前三个安装不好可以安装他

方法一：自行创造函数法，较复杂。新手建议方法二

此法使用Cibersort工具需要三个文件：
1、sourcecibersort.R
2、LM22.txt
3、genes_exp.txt

1.sourcecibersort.R

直接把下列代码新建一个script，然后保存，保存名字为sourcecibersort.R

  1. #' CIBERSORT R script v1.03 (last updated 07-10-2015)
  2. #' Note: Signature matrix construction is not currently available; use java version for full functionality.
  3. #' Author: Aaron M. Newman, Stanford University (amnewman@stanford.edu)
  4. #' Requirements:
  5. #'       R v3.0 or later. (dependencies below might not work properly with earlier versions)
  6. #'       install.packages('e1071')
  7. #'       install.pacakges('parallel')
  8. #'       install.packages('preprocessCore')
  9. #'       if preprocessCore is not available in the repositories you have selected, run the following:
  10. #'           source("http://bioconductor.org/biocLite.R")
  11. #'           biocLite("preprocessCore")
  12. #' Windows users using the R GUI may need to Run as Administrator to install or update packages.
  13. #' This script uses 3 parallel processes.  Since Windows does not support forking, this script will run
  14. #' single-threaded in Windows.
  15. #'
  16. #' Usage:
  17. #'       Navigate to directory containing R script
  18. #'
  19. #'   In R:
  20. #'       source('CIBERSORT.R')
  21. #'       results <- CIBERSORT('sig_matrix_file.txt','mixture_file.txt', perm, QN)
  22. #'
  23. #'       Options:
  24. #'       i)  perm = No. permutations; set to >=100 to calculate p-values (default = 0)
  25. #'       ii) QN = Quantile normalization of input mixture (default = TRUE)
  26. #'
  27. #' Input: signature matrix and mixture file, formatted as specified at http://cibersort.stanford.edu/tutorial.php
  28. #' Output: matrix object containing all results and tabular data written to disk 'CIBERSORT-Results.txt'
  29. #' License: http://cibersort.stanford.edu/CIBERSORT_License.txt
  30. #' Core algorithm
  31. #' @param X cell-specific gene expression
  32. #' @param y mixed expression per sample
  33. #' @export
  34. CoreAlg <- function(X, y){
  35.   
  36.   #try different values of nu
  37.   svn_itor <- 3
  38.   
  39.   res <- function(i){
  40.     if(i==1){nus <- 0.25}
  41.     if(i==2){nus <- 0.5}
  42.     if(i==3){nus <- 0.75}
  43.     model<-e1071::svm(X,y,type="nu-regression",kernel="linear",nu=nus,scale=F)
  44.     model
  45.   }
  46.   
  47.   if(Sys.info()['sysname'] == 'Windows') out <- parallel::mclapply(1:svn_itor, res, mc.cores=1) else
  48.     out <- parallel::mclapply(1:svn_itor, res, mc.cores=svn_itor)
  49.   
  50.   nusvm <- rep(0,svn_itor)
  51.   corrv <- rep(0,svn_itor)
  52.   
  53.   #do cibersort
  54.   t <- 1
  55.   while(t <= svn_itor) {
  56.     weights = t(out[[t]]$coefs) %*% out[[t]]$SV
  57.     weights[which(weights<0)]<-0
  58.     w<-weights/sum(weights)
  59.     u <- sweep(X,MARGIN=2,w,'*')
  60.     k <- apply(u, 1, sum)
  61.     nusvm[t] <- sqrt((mean((k - y)^2)))
  62.     corrv[t] <- cor(k, y)
  63.     t <- t + 1
  64.   }
  65.   
  66.   #pick best model
  67.   rmses <- nusvm
  68.   mn <- which.min(rmses)
  69.   model <- out[[mn]]
  70.   
  71.   #get and normalize coefficients
  72.   q <- t(model$coefs) %*% model$SV
  73.   q[which(q<0)]<-0
  74.   w <- (q/sum(q))
  75.   
  76.   mix_rmse <- rmses[mn]
  77.   mix_r <- corrv[mn]
  78.   
  79.   newList <- list("w" = w, "mix_rmse" = mix_rmse, "mix_r" = mix_r)
  80.   
  81. }
  83. #' do permutations
  84. #' @param perm Number of permutations
  85. #' @param X cell-specific gene expression
  86. #' @param y mixed expression per sample
  87. #' @export
  88. doPerm <- function(perm, X, Y){
  89.   itor <- 1
  90.   Ylist <- as.list(data.matrix(Y))
  91.   dist <- matrix()
  92.   
  93.   while(itor <= perm){
  94.     #print(itor)
  95.     
  96.     #random mixture
  97.     yr <- as.numeric(Ylist[sample(length(Ylist),dim(X)[1])])
  98.     
  99.     #standardize mixture
  100.     yr <- (yr - mean(yr)) / sd(yr)
  101.     
  102.     #run CIBERSORT core algorithm
  103.     result <- CoreAlg(X, yr)
  104.     
  105.     mix_r <- result$mix_r
  106.     
  107.     #store correlation
  108.     if(itor == 1) {dist <- mix_r}
  109.     else {dist <- rbind(dist, mix_r)}
  110.     
  111.     itor <- itor + 1
  112.   }
  113.   newList <- list("dist" = dist)
  114. }
  115.  
  116. #' Main functions
  117. #' @param sig_matrix file path to gene expression from isolated cells
  118. #' @param mixture_file heterogenous mixed expression
  119. #' @param perm Number of permutations
  120. #' @param QN Perform quantile normalization or not (TRUE/FALSE)
  121. #' @export
  122. CIBERSORT <- function(sig_matrix, mixture_file, perm=0, QN=TRUE){
  123.   
  124.   #read in data
  125.   X <- read.table(sig_matrix,header=T,sep="\t",row.names=1,check.names=F)
  126.   Y <- read.table(mixture_file, header=T, sep="\t", row.names=1,check.names=F)
  127.   
  128.   X <- data.matrix(X)
  129.   Y <- data.matrix(Y)
  130.   
  131.   #order
  132.   X <- X[order(rownames(X)),]
  133.   Y <- Y[order(rownames(Y)),]
  134.   
  135.   P <- perm #number of permutations
  136.   
  137.   #anti-log if max < 50 in mixture file
  138.   if(max(Y) < 50) {Y <- 2^Y}
  139.   
  140.   #quantile normalization of mixture file
  141.   if(QN == TRUE){
  142.     tmpc <- colnames(Y)
  143.     tmpr <- rownames(Y)
  144.     Y <- preprocessCore::normalize.quantiles(Y)
  145.     colnames(Y) <- tmpc
  146.     rownames(Y) <- tmpr
  147.   }
  148.   
  149.   #intersect genes
  150.   Xgns <- row.names(X)
  151.   Ygns <- row.names(Y)
  152.   YintX <- Ygns %in% Xgns
  153.   Y <- Y[YintX,]
  154.   XintY <- Xgns %in% row.names(Y)
  155.   X <- X[XintY,]
  156.   
  157.   #standardize sig matrix
  158.   X <- (X - mean(X)) / sd(as.vector(X))
  159.   
  160.   #empirical null distribution of correlation coefficients
  161.   if(P > 0) {nulldist <- sort(doPerm(P, X, Y)$dist)}
  162.   
  163.   #print(nulldist)
  164.   
  165.   header <- c('Mixture',colnames(X),"P-value","Correlation","RMSE")
  166.   #print(header)
  167.   
  168.   output <- matrix()
  169.   itor <- 1
  170.   mixtures <- dim(Y)[2]
  171.   pval <- 9999
  172.   
  173.   #iterate through mixtures
  174.   while(itor <= mixtures){
  175.     
  176.     y <- Y[,itor]
  177.     
  178.     #standardize mixture
  179.     y <- (y - mean(y)) / sd(y)
  180.     
  181.     #run SVR core algorithm
  182.     result <- CoreAlg(X, y)
  183.     
  184.     #get results
  185.     w <- result$w
  186.     mix_r <- result$mix_r
  187.     mix_rmse <- result$mix_rmse
  188.     
  189.     #calculate p-value
  190.     if(P > 0) {pval <- 1 - (which.min(abs(nulldist - mix_r)) / length(nulldist))}
  191.     
  192.     #print output
  193.     out <- c(colnames(Y)[itor],w,pval,mix_r,mix_rmse)
  194.     if(itor == 1) {output <- out}
  195.     else {output <- rbind(output, out)}
  196.     
  197.     itor <- itor + 1
  198.     
  199.   }
  200.   
  201.   #save results
  202.   write.table(rbind(header,output), file="CIBERSORT-Results.txt", sep="\t", row.names=F, col.names=F, quote=F)
  203.   
  204.   #return matrix object containing all results
  205.   obj <- rbind(header,output)
  206.   obj <- obj[,-1]
  207.   obj <- obj[-1,]
  208.   obj <- matrix(as.numeric(unlist(obj)),nrow=nrow(obj))
  209.   rownames(obj) <- colnames(Y)
  210.   colnames(obj) <- c(colnames(X),"P-value","Correlation","RMSE")
  211.   obj
  212. }

2.LM22.txt（文末网盘），用来作参考的免疫浸润数据集，也可以官网下载，是一样的。

3.genes_exp.txt（文末网盘），实际上就是你要处理的自己的数据或者TCGA的数据

4.启动

  1. source("sourcecibersort.R") #启动这个函数，必须在哦那个一个文件夹内才可哟
  2. results <- CIBERSORT(sig_matrix ="LM22.txt", mixture_file ="genes_exp.txt", perm = 1000, QN = T)
  3. # perm置换次数=1000，QN分位数归一化=TRUE
  4. # 文件名可以自定义
  5. # 得到的结果可以用来绘制热图等等

5.之后的分析和方法二一样，都是对results进行绘图

方法二：使用打包的函数法，这个需要下载R包，但不需要自己创造函数了

1.还是运行以下代码，多安了个包

  1. # install packages 这三个安装不成功的话，就安后面的bseqsc包也行
  2. install.packages('e1071')
  3. install.pacakges('parallel')
  4. install.packages('preprocessCore')
  5. library(e1071)
  6. library(preprocessCore)
  7. library(parallel)
  8.  
  9. install.packages('devtools')
  10. library(devtools)
  11. devtools::install_github('shenorrlab/bseqsc')
  12. library(bseqsc)#这个包携带大量CIBERSORT的依赖，前三个安装不好可以安装他
  13.  
  14. ################安装CIBERSORT包##########################################################
  15. if(!require(CIBERSORT))devtools::install_github("Moonerss/CIBERSORT")
  16. library(CIBERSORT)
  17. # 包全部安装完成
  18.  
  19. # 画热图的包
  20. install.packages("pheatmap")
  21. install.packages("ComplexHeatmap")
  22. library(ggplot2)
  23. library(pheatmap)
  24. library(ComplexHeatmap)

安装好以后就可以使用cibersort函数了

2.它的好处在于自带了LM22文件和测试文件，你不需要额外去下载和整理了。我本人已经测试过这个LM22和自己读取的LM22完全一样。直接运行下列code即可出图

  1. # 同时准备好LM22的TXT文件，注意自己以后的文件要和这个TXT的格式一样
  2. # 加载CIBERSORT包成功后，系统内部会自带data(LM22)
  3. data(LM22) 
  4. data(mixed_expr)#TCGA的演示数据，正式情况下就用自己的数据
  5.  
  6. # 正式开始探索
  7. # 看5*5的数据
  8. LM22[1:5,1:5]
  9. mixed_expr[1:5,1:5]
  10.  
  11. # 分别定义signature矩阵LM22和我的数据（演示）矩阵mixed_expr
  12. results <- cibersort(sig_matrix = LM22, mixture_file = mixed_expr)
  13.  
  14. # 理解一下results的结果
  15. # 你可以理解为返回了一个列名为细胞类型、行名为样本名的细胞浸润程度（占比）的矩阵
  16. # 此外result中还会多出三列：
  17. # P-value: 用来展示去卷积的结果在所有细胞类群中是否具有差异
  18. # Correlation:参考矩阵与输入矩阵的特征基因相关性
  19. # RMSE: Root mean squared error，参考矩阵与输入矩阵的特征基因标准差
  20.  
  21. # heatmap
  22. # 按行（样本内部）标准化可以看出在各类样本内部，M2浸润程度（占比）最高
  23. rowscale <- results[,1:ncol(LM22)]#只是相当于备份了一下results
  24. rowscale <- rowscale[,apply(rowscale, 2, function(x){sum(x)>0})]#删除全是0的列
  25. pheatmap(rowscale,
  26.          scale = 'row',#按行标准化，不标准化就会按绝对值显示，很诡异
  27.          cluster_col=T,#是否对列聚类，不聚类，坐标轴就按照原来的顺序显示
  28.          cluster_row=F,#是否对行聚类
  29.          angle_col = "315")#调整X轴坐标的倾斜角度
  30.  
  31. # 各类样本之间也具有自己占比高的特异性免疫细胞
  32. columnscale <- results[,1:ncol(LM22)]
  33. columnscale <- columnscale[,apply(columnscale, 2, function(x){sum(x)>0})]#删除全是0的列
  34. pheatmap(columnscale,
  35.          scale = 'column',
  36.          cluster_col=F,
  37.          cluster_row=T,
  38.          angle_col = "315")
  39.  
  40. # 堆积比例图
  41. my36colors <-c('#E5D2DD', '#53A85F', '#F1BB72', '#F3B1A0', '#D6E7A3', '#57C3F3', '#476D87','#E95C59', '#E59CC4', '#AB3282', '#23452F', '#BD956A', '#8C549C', '#585658','#9FA3A8', '#E0D4CA', '#5F3D69', '#C5DEBA', '#58A4C3', '#E4C755', '#F7F398','#AA9A59', '#E63863', '#E39A35', '#C1E6F3', '#6778AE', '#91D0BE', '#B53E2B', '#712820', '#DCC1DD', '#CCE0F5',  '#CCC9E6', '#625D9E', '#68A180', '#3A6963','#968175'
  42. )
  43. cellnum <- results[,1:ncol(LM22)]
  44. cell.prop<- apply(cellnum, 1, function(x){x/sum(x)})
  45. data4plot <- data.frame()
  46. for (i in 1:ncol(cell.prop)) {
  47.   data4plot <- rbind(
  48.     data4plot,
  49.     cbind(cell.prop[,i],rownames(cell.prop),
  50.           rep(colnames(cell.prop)[i],nrow(cell.prop)
  51.           )
  52.     )
  53.   )
  54. }
  55.  
  56. colnames(data4plot)<-c('proportion','celltype','sample')
  57. data4plot$proportion <- as.numeric(data4plot$proportion)
  58. ggplot(data4plot,aes(sample,proportion,fill=celltype))+
  59.   geom_bar(stat="identity",position="fill")+
  60.   scale_fill_manual(values=my36colors)+#自定义fill的颜色
  61.   ggtitle("cell portation")+
  62.   theme_bw()+
  63.   theme(axis.ticks.length=unit(0.5,'cm'),axis.title.x=element_text(size=1))+
  64.   theme(axis.text.x = element_text(angle = 45, hjust = 0.5, vjust = 0.5))+#把x坐标轴横过来
  65.   guides(fill=guide_legend(title=NULL))

 

LM22：

#########链接：https://pan.baidu.com/s/1eQSEekekozS5osgydwzk1w
#@####提取码：fk88
 

  1. LM22read <- read.csv("LM22.csv",header  = T)
  2. gene <- LM22read[,1]
  3. rownames(LM22read) <- gene
  4. LM22read <- LM22read[,-1]
  5. data(LM22)
  6. all(LM22==LM22read)#可以看到TURE，说明两个文件完全一样了；LM22是上文提到的安装CIBERSORT包之后自带的data

鸣谢：生信技能树jimmy老师和 Biomamba 生信基地 BIOMAMBA老师

有疑问可以邮件联系我，会尽力帮忙：yunbk@mail2.sysu.edu.cn

2023.03.05更新

结合不少小伙伴的私信和邮件更新几个点

1.自己的表达矩阵格式：按照作者的文档和示例数据，应该是不取log的，也不能做其他的normalizie处理。当然如果是log处理过的，cibersort会自己去log，所以也不必担心。数据应该是标准化后的数据比如FPKM TPM CPM这样，不应该用count。

2.在数据量比较大，样本多的时候，请下载R包使用cibersort，因为source函数的读取速度会变得很慢。

3.计算前需要排除空值（NA）值，不然会报错

4.最新的版本是cibersortX,但是好像现在维护有点问题，大家可以自行搜索一下这个方法

声明：本文内容由网友自发贡献，不代表【wpsshop博客】立场，版权归原作者所有，本站不承担相应法律责任。如您发现有侵权的内容，请联系我们。转载请注明出处：https://www.wpsshop.cn/w/繁依Fanyi0/article/detail/221978
推荐阅读