CIBERSORT计算免疫细胞丰度

最新推荐文章于 2023-08-14 00:26:30 发布

arrhythmia08

最新推荐文章于 2023-08-14 00:26:30 发布

阅读量2.2k

点赞数 2

文章标签： r语言

本文链接：https://blog.csdn.net/arrhythmia08/article/details/128376730

版权

我也是第一次做cibersort分析，刚好看见前辈们写的教程特别好，所以自己总结一下我所需要的，写成教程主要是方便自己。

cibersort计算免疫细胞丰度总共只需要一个函数+一行代码，主要工作在于文件的准备。总共需要准备三个文件。

第一，函数代码。cibersort没有开发成一个包，所以需要复制代码到R，保存为"CIBERSORT.R"的脚本到工作目录下，需要使用时引用该代码即可。该代码脚本网上资源很多，百度就能搜到。

"CIBERSORT.R"脚本前面有一段注释，使用该需要安装几个包。自行检查，该安装该加载按要求做。

第二，准备22个免疫基因的表达文件“LM22.txt”。

Robust enumeration of cell subsets from tissue expression profiles | Nature Methods在原文章下下载supplementary table1，如下如所示:

找到第一个工作簿，从“gene symbol”列复制，将22个免疫细胞和基因表达文件复制到新的excel表，并保存为“LM22.txt”的的文本文件到工作目录下，如图片所示：

第三，接下来准备我们自己的表达矩阵文件

cibersort对输入的表达矩阵有一定要求，原文如下：

1，不可有负值，缺失值；没有经过log处理。

2，affymetrix的芯片，推荐经过RAM或者MAS5.0标准化。

3，illumina的beadchip和单色安捷伦芯片需要limma包处理。（这部分处理方式我暂不了解，有需要的推荐孙小洁老师的教程，贴个链接）illumina beadchip 芯片原始数据处理 - 简书

4，测序数据可以使用FPKM和TPM。

我接下来要使用的示例数据GSE10846,affymetrix芯片数据，经过MAS5.0标准化。作者都会在data_process一列中标注数据处理方式。mitrx文件我已经在GEO下载好，下面我们用代码来提取表达矩阵：

#清空环境变量
rm(list = ls())

#获取eSet
library(GEOquery)
eSet <- getGEO(GEO = "GSE10846", 
               filename = "GSE10846_series_matrix.txt.gz",
               destdir = '.', 
               getGPL = F)

#提取表达矩阵和测序平台GPL570
exp <- exprs(eSet);dim(exp)
gpl_number <- eSet@annotation;gpl_number

#ID转换
library(tinyarray)
library(AnnoProbe)
ids <- AnnoProbe::idmap('GPL570')
exp <- trans_array(exp,ids)
exp[1:4,1:4]

#提取临床信息
pd <- pData(eSet)
View(pd)

#小习惯，随时检查表达矩阵和临床信息是否一致
identical(colnames(exp),rownames(pd))

save(exp,pd,file="exp_pd.Rdata")

从pd的data_precossing这一列我们可以看到该数据经过MAS5.0标准化，可以用于cibersort分析。

#加载刚保存的变量
rm(list = ls())
load("exp_pd.Rdata")

#检查是否有缺失值和负值
table(is.na(exp))
# FALSE 
#8478960 

range(exp_origin)
#[1]  0.000 17.825

#如果有缺失值
exp <- na.omit(exp)

#此处将行名基因symbol变为表达矩阵的一列，并存为txt文件
library(dplyr)
exp <- as.data.frame(exp)
exp <- rownames_to_column(exp)
write.table(exp,file = "exp.txt",row.names = F,quote = F,sep = "\t")

用excel打开刚刚写入的表达矩阵，将第一列的列名由"rowname",改为“Gene symbol”，和LM22.txt第一列的列名完全一致，否则等会儿可能报错。

我们三个文件都准备完成，可以计算免疫细胞丰度了。

#第一个文件这里需要用到
source("CIBERSORT.R")

#这里计算时间比较久，耐心等待一下
TME.results = CIBERSORT("LM22.txt", 
                        "exp.txt" , 
                        perm = 1000, 
                        QN = T)

colnames(TME.results)  #可以看到前22列为免疫细胞，后面三列不需要
immu_cell <- TME.results[,-(23:25)]
immu_cell[1:4,1:4]

#          B cells naive B cells memory Plasma cells T cells CD8
#GSM274895    0.01736808     0.12164276            0  0.07648905
#GSM274896    0.03402937     0.04344233            0  0.09497394
#GSM274897    0.00000000     0.16178958            0  0.07226354
#GSM274898    0.22320799     0.32301388            0  0.03555337

这样就将22种免疫细胞丰度计算出来，接下来想画箱线图，热图，直方图全凭需要。

参考：孙小洁老师CIBERSORT 免疫浸润 · 语雀

【块】Cibersort.R计算22种免疫细胞浸润分数 - 简书