Cibersort实战

rm(list=ls())

load(file='../Rdata/@step3_DEG_LN.Rdata')
# exprSet=new_exprSet
# phe<- data.frame(ID=pheno$geo_accession,gender=pheno$`gender:ch1`,age=pheno$`age:ch1`,
#                  differentiation=pheno$`his type:ch1`,T_stage= pheno$`depth:ch1`,
#                  N_stage=pheno$`lymph node metastasis:ch1`)

# phe<- phe[is.na(phe$gender)==F,]
# phe<- phe[phe$gender!='NA',]
# phe<- phe[phe$N_stage!='NA',]
phe$N_stage<- factor(as.character(phe$N_stage))
exprSet<- exprSet[,colnames(exprSet)%in%phe$GSM_ID]
dim(exprSet)
colnames(phe)
# exprSet<- exprSet[,1:5]
exprSet<- exprSet[order(row.names(exprSet)),]
# exprSet[1:5,1:5]
match(colnames(exprSet),phe$GSM_ID)
##### 基因如果存在一些数字开头的名字 可能会出错
# exprSet<- exprSet[-(1:2),]

#### 如果实在不行就只提取500个基因进行计算
# ##### cibersort 计算只需要500多个基因就可以了。
# LM22 <- read.delim("~/R Project/ESCC-undermine/GSE47404/LN_related_genes/LM22.txt")
# exprSet<- exprSet[rownames(exprSet)%in%LM22$Gene.symbol,]
# write.table(exprSet,file = '../Rdata/for_cibersort_mitrix.txt',sep = '\t',quote=FALSE)

###### 输入结果后作图
library(reshape2)
library(ggthemes)
library(ggpubr)
library("scales")
cibersort_result<- read.csv(file = '../Rdata/CIBERSORT.Output_Job13.csv')
colnames(cibersort_result)
match(cibersort_result$Input.Sample,phe$GSM_ID)
cibersort_result<- cibersort_result[,-c(1,24:26)]
### 只挑选一部分细胞
cibersort_result<- cibersort_result[,c(1,4,7,9,14,15,16,18,19,20,22)]

gene_name<- 'GPD1L'
exprSet<- as.matrix(exprSet)
cibersort_result$group<- ifelse(exprSet[gene_name,]>quantile(exprSet[gene_name,],0.5),'High','Low')
cibersort_result_L<-melt(cibersort_result,id.vars = c('group'),variable.name="TILs",value.name="value")
library(ggplot2)
library(stringr)
p=ggplot(cibersort_result_L,aes(x=TILs,y=value,fill=group))+geom_boxplot()+
  stat_compare_means(method = "wilcox.test",label="p.signif")
  # scale_x_discrete(labels=colnames(cibersort_result))+
  # theme(axis.title  = element_text(,angle = 45,hjust = .5, vjust = .5))
print(p)
ggsave('../figure/TILs_ggplot.pdf', p)

ggplot(cibersort_result_L,aes(group,value,fill=TILs,colour = TILs))+
  geom_bar(stat="identity",position="fill")+
  ggtitle("LM22")+
  theme_classic()+
  xlab(gene_name)+
  scale_fill_manual(values=c("#999999", "#E69F00", "#56B4E9","black","red","chocolate","deeppink","darkmagenta",
                             "#000000", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7",'green','#FFB6C1',
                             '#DC143C','#DDA0DD','#800080','#6A5ACD','#87CEEB','#BDB76B','#FF4500'))+
  scale_colour_manual(values=c("#999999", "#E69F00", "#56B4E9","black","red","chocolate","deeppink","darkmagenta",
                               "#000000", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7",'green','#FFB6C1',
                               '#DC143C','#DDA0DD','#800080','#6A5ACD','#87CEEB','#BDB76B','#FF4500'))