Xcell实战

xCell is a recently published method based on ssGSEA that estimates the abundance scores of 64 immune cell types, including adaptive and innate immune cells, hematopoietic progenitors, epithelial cells, and extracellular matrix cells

关于Xcell找对网址很重要，我一开始找错了地方

install.packages('Rcpp')#########安装各类程序包
devtools::install_github('dviraran/xCell')

安装的时候还会有错误。

安装好的这一刻，还是很开心的。

使用方法

第一步 计算xCell

library(xCell)
exprMatrix = read.table(file = '/Users/chenyuqiao/Desktop/TCGA-LUAD.htseq_counts.tsv',header=TRUE,row.names=1, as.is=TRUE)
xCellAnalysis(exprMatrix)

library(xCell)
exprMatrix = read.table(file = '/Users/chenyuqiao/Desktop/TCGA-LUAD.htseq_counts.tsv',header=TRUE,row.names=1, as.is=TRUE)

###exprMatrix<- exprMatrix[1:10,1:10]
Ensemble_ID<- rownames(exprMatrix)
ID<- strsplit(Ensemble_ID, "[.]")
str(ID)
IDlast<- sapply(ID, "[", 1)
exprMatrix$Ensemble_ID<- IDlast
row.names(exprMatrix)<- exprMatrix$Ensemble_ID
save(exprMatrix, file = 'TCGA.Rdata')
load(file = 'TCGA.Rdata')


####library(clusterProfiler)
library(org.Hs.eg.db)
ls("package:org.Hs.eg.db")
g2s=toTable(org.Hs.egSYMBOL);head(g2s)
g2e=toTable(org.Hs.egENSEMBL);head(g2e)
tmp=merge(g2e,g2s,by='gene_id')
head(tmp)
colnames(exprMatrix)[ncol(exprMatrix)] <- c("ensembl_id")###################重命名Ensemble_ID 便于后面merge
exprMatrix[1:4,1:4]
exprMatrix<- merge(tmp,exprMatrix,by='ensembl_id')
exprMatrix[1:4,1:4]
exprMatrix<- exprMatrix[,- c(1,2)]
exprMatrix=exprMatrix[!duplicated(exprMatrix$symbol),]
row.names(exprMatrix)<- exprMatrix[,1]
exprMatrix<- exprMatrix[,-1]
exprMatrix[1:4,1:4]
xCellAnalysis(exprMatrix)####################一句话就分析完成了
##save(results,file = 'Xcell_result.Rdata')#############需要重新修改

第二步：批量生存分析

load(file = 'Xcell_result.Rdata')
result<- as.data.frame(result)
library(dplyr)
library(tidyverse)

TCGA.LUAD.GDC_phenotype <- read.delim("TCGA-LUAD.GDC_phenotype.tsv")

#colnames(TCGA.LUAD.GDC_phenotype)
#head(TCGA.LUAD.GDC_phenotype)

LUAD_Pheno<- select(TCGA.LUAD.GDC_phenotype, "submitter_id.samples", "vital_status.diagnoses", "days_to_death.diagnoses", "days_to_last_follow_up.diagnoses", "pathologic_N", "pathologic_M", "days_to_new_tumor_event_after_initial_treatment")
LUAD_Pheno<- LUAD_Pheno[grep('01A',LUAD_Pheno$submitter_id.samples),]  #####只筛选01A的   01A代表肿瘤
LUAD_Pheno[is.na(LUAD_Pheno)]<- 0
LUAD_Pheno$PFS_status<- ifelse((LUAD_Pheno$days_to_new_tumor_event_after_initial_treatment == 0 & LUAD_Pheno$days_to_death.diagnoses == 0), 0,1)
##################################
LUAD_Pheno$OS<- ifelse(LUAD_Pheno$days_to_last_follow_up.diagnoses > LUAD_Pheno$days_to_death.diagnoses, LUAD_Pheno$days_to_last_follow_up.diagnoses,LUAD_Pheno$days_to_death.diagnoses)
LUAD_Pheno$PFS<- ifelse(LUAD_Pheno$days_to_new_tumor_event_after_initial_treatment == 0, LUAD_Pheno$OS ,LUAD_Pheno$days_to_new_tumor_event_after_initial_treatment)
LUAD_Pheno$OS_status<- as.factor(LUAD_Pheno$vital_status.diagnoses)
#############################设计好分组




#############################生存曲线

firstdata<- result  ###############expre
firstdata$ID<- rownames(firstdata)
gene<- row.names(firstdata)
#######select only gene to analysis
library(survminer)
library(survival)
library(ggplot2)
library(dplyr)
for (x in gene) {
  RNA_seq_data<-filter(firstdata, firstdata$ID == x)
  RNA_seq_data<- t(RNA_seq_data)
  RNA_seq_data<- as.data.frame(RNA_seq_data)
  # str(RNA_seq_data)
  # colnames(LUAD_Pheno)
  RNA_seq_data$submitter_id.samples<- row.names(RNA_seq_data)
  colnames(RNA_seq_data)<- c("Expressionvalue","submitter_id.samples")
  LUAD_Pheno$submitter_id.samples<- as.character(LUAD_Pheno$submitter_id.samples)
  LUAD_Pheno$submitter_id.samples<- sub('-', '.', LUAD_Pheno$submitter_id.samples)#############- replaced by .
  LUAD_Pheno$submitter_id.samples<- sub('-', '.', LUAD_Pheno$submitter_id.samples)#############- replaced by .
  LUAD_Pheno$submitter_id.samples<- sub('-', '.', LUAD_Pheno$submitter_id.samples)#############- replaced by .
  LUAD_Pheno$submitter_id.samples<- sub('-', '.', LUAD_Pheno$submitter_id.samples)#############- replaced by .
  finaldata<- inner_join(LUAD_Pheno,RNA_seq_data, by = "submitter_id.samples")
  finaldata$PFS_status<- as.character(finaldata$PFS_status)
  finaldata$PFS_status<- as.numeric(as.factor(finaldata$PFS_status))
  finaldata$Expressionvalue<- as.numeric(as.character(finaldata$Expressionvalue))
  finaldata$group<- ifelse(finaldata$Expressionvalue>median(finaldata$Expressionvalue),'high','low')
  library(survminer)
  library(survival)
  fit <- survfit(Surv(finaldata$PFS,finaldata$PFS_status)~finaldata$group, data=finaldata) 
  summary(fit)
  pp<- ggsurvplot(fit, data = finaldata, conf.int = F, pval = TRUE,
                  xlim = c(0,2000), # present narrower X axis, but not affect
                  # survival estimates. 
                  xlab = "Time in days", # customize X axis label. 
                  break.time.by = 200) # break X axis in time intervals by 500. 
  ggsave(filename = paste("plot_",x,".pdf",sep = ""))
  print(x)
}