读CIBERSORT免疫浸润的代码实现、cibersoft使用SVM算法实现去卷积有感,知识连贯性即系统性最终决定你能否得到自己想要的结果。
cybersort
Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. 2015 May;12(5):453-7. doi: 10.1038/nmeth.3337. Epub 2015 Mar 30. PMID: 25822800; PMCID: PMC4739640.
二、输入数据类型
Importantly, all expression data should be non-negative, devoid of missing values, and represented in non-log linear space.
For Affymetrix microarrays, a custom chip definition file (CDF) is recommended (see Subheading 3.2.2) and should be normalized with MAS5 or RMA.
Illumina Beadchip and single color Agilent arrays should be processed as described in the limma package.
Standard RNA-Seq expression quantification metrics, such as frag- ments per kilobase per million (FPKM) and transcripts per kilobase million (TPM), are suitable for use with CIBERSORT. –《Profiling Tumor Infiltrating Immune Cells with CIBERSORT》
首先,所有的数据要求非负,无缺失值,未经过log转化
1、Affymetrix microarrays芯片数据要求 MAS5 or RMA标准化
2、Illumina Beadchip 需要经limma package处理
3、RNA-seq数据可以使用FPKM和TPM
TCGA中肝癌数据进行演示:
rm(list = ls())
library(tidyverse)
library(ggplot2)
library(stringr)
load("../3.DEGs/tpm_exp.Rdata")
load("../3.DEGs/DEG_limma.Rdata")
exp <- 2^tpm_exp-1
exp <- exp[,str_sub(colnames(exp),1,15)%in%group$Sample.ID]
exp <- exp[,match(group$Sample.ID,str_sub(colnames(exp),1,15))] %>%
rownames_to_column()
write.table(exp,file = "exp.txt",sep = "\t",row.names = F,quote = F)
我直接在UCSC中下载了FPKM的数据,由于我之前已经转化未log2(TPM+1),这里为了满足要求,我直接转化未TPM,将行名转换为列,保存未txt文件备用。
rm(DEG_limma,exp,tpm_exp)
if(F){
library(e1071)
source("CIBERSORT.R")
TME.results = CIBERSORT("LM22.txt",
"exp.txt" ,
perm = 1,
QN = T)
save(TME.results,file = "ciber.Rdata")
}
load("ciber.Rdata")
我本想先使用了374个样本运行代码,再挑选我需要的样本免疫细胞浸润情况,无奈运行时间有点长,让我误以为自己电脑坏了/////////所以,这里我挑选了73个样本直接运行,速度还可以,不到一分钟搞定,得到矩阵。
三、数据展示
1、免疫细胞丰度热图
#--热图---------------------------------------------------------------------------------------------
re <- TME.results[,-(23:25)]
library(pheatmap)
k <- apply(re,2,function(x) {sum(x == 0) < nrow(TME.results)/2})
table(k)
## k
## FALSE TRUE
## 9 13
re2 <- as.data.frame(t(re[,k]))
an = data.frame(group = group$Disease.Free.Status,
row.names = colnames(re2))
pheatmap(re2,scale = "row",
show_colnames = F,cluster_cols = F,
annotation_col = an,breaks = seq(-2,2,length.out = 50),
color = colorRampPalette(c("navy", "white", "firebrick3"))(50))
2、直方图
#---barplot---------------------------------------------------------------------------------------------
library(RColorBrewer)
mypalette <- colorRampPalette(brewer.pal(8,"Set1"))
dat <- re %>% as.data.frame() %>%
rownames_to_column("Sample") %>%
gather(key = Cell_type,value = Proportion,-Sample)
ggplot(dat,aes(Sample,Proportion,fill = Cell_type)) +
geom_bar(stat = "identity") +
labs(fill = "Cell Type",x = "",y = "Estiamted Proportion") +
theme_bw() +
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
legend.position = "bottom") +
scale_y_continuous(expand = c(0.01,0)) +
scale_fill_manual(values = mypalette(22))
3、 箱线图
#------------------------------------------------------------------------
a = dat %>%
group_by(Cell_type) %>%
summarise(m = median(Proportion)) %>%
arrange(desc(m)) %>%
pull(Cell_type)
dat$Cell_type = factor(dat$Cell_type,levels = a)
ggplot(dat,aes(Cell_type,Proportion,fill = Cell_type)) +
geom_boxplot(outlier.shape = 21,color = "black") +
theme_bw() +
labs(x = "Cell Type", y = "Estimated Proportion") +
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
legend.position = "bottom") +
scale_fill_manual(values = mypalette(22))
#----------------------------------------------------------------------------
dat$Group = factor(rep(group$Disease.Free.Status,times=22),levels = c(0,1))
library(ggpubr)
ggplot(dat,aes(Cell_type,Proportion,fill = Group)) +
geom_boxplot(outlier.shape = 21,color = "black") +
theme_bw() +
labs(x = "Cell Type", y = "Estimated Proportion") +
theme(legend.position = "top") +
theme(axis.text.x = element_text(angle=60,hjust =1))+
scale_fill_manual(values = mypalette(22)[c(6,1)])+
stat_compare_means(aes(group = Group,label = ..p.signif..),method = "t.test")