最近學習了生信菜鳥團的純R代碼實現(xiàn)ssGSEA算法評估腫瘤免疫浸潤程度,想復(fù)制作者的流程猿挚，但是發(fā)現(xiàn)了幾個不一樣的地方，所以重新整理流程驶鹉，代碼主要來自原作者Juan_NF绩蜻。
文章來源：Local mutational diversity drives intratumoral immune heterogeneity in non-small cell lung cancer
方法來源：
Pan-cancer Immunogenomic Analyses Reveal Genotype-Immunophenotype Relationships and Predictors of Response to Checkpoint Blockade
研究背景：方法基于使用metagenes，即代表特定免疫細胞亞群的非重疊基因組室埋，并且既不在CRC細胞系中也不在正常組織中表達办绝。然后使用這些組的metagenes的表達來使用基因集富集分析（GSEA）分析統(tǒng)計富集。優(yōu)點是該方法的穩(wěn)健性姚淆，這是由于兩個特征：
（1）使用一組基因而不是代表一個免疫亞群的單個基因孕蝉，因為使用單個基因作為免疫亞群的標記可以是誤導(dǎo)因為許多基因在不同的細胞類型中表達;
（2）評估一組基因相對于樣品中所有其他基因表達的相對表達變化。
為28個免疫細胞亞群定義了一組泛癌癥腌逢，并將分析擴展到實體癌癥降淮。對具有> 8,000個腫瘤樣品的20個實體癌癥的TCGA數(shù)據(jù)進行了免疫原性表征，并提供了腫瘤內(nèi)免疫浸潤的細胞組成的全面視圖搏讶。此外骤肛，還得出了癌癥抗原以及個體樣本的遺傳特征（腫瘤異質(zhì)性和克隆性）纳本，以便能夠?qū)γ庖咛卣饕约澳[瘤的遺傳特征進行綜合分析。開發(fā)了數(shù)據(jù)庫TCIA（癌癥免疫組圖譜）腋颠。基于反卷積方法識別免疫亞群的分數(shù)（CIBERSORT;
紐曼等人吓笙，2015年）在TCIA網(wǎng)站上提供了GSEA和反卷積數(shù)據(jù)淑玫。
1.提取矩陣和表型信息，需要手動從GEO下載
下載GSE112996_merged_fpkm_table.txt.gz
GSE112996_series_matrix.txt.gz面睛，這兩個文件絮蒿，對
GSE112996_series_matrix.txt.gz進行解壓，把這兩個文件放到Rproject創(chuàng)建的文件夾叁鉴。
rm(list=ls())
a <- read.table('GSE112996_merged_fpkm_table.txt.gz',
header = T,
row.names=1)
raw_data<- a[,-1]
###表型信息提取
pheno <- read.csv(file = 'GSE112996_series_matrix.txt')
pheno <- data.frame(num1 = strsplit(as.character(pheno[42,]),split='\t')[[1]][-1],
num2 = gsub('patient: No.','P',strsplit(as.character(pheno[51,]),split='\t')[[1]][-1]))
{
####數(shù)據(jù)過濾
data<- a[!apply(raw_data,1,sum)==0,]
####去除重復(fù)基因名的行土涝，歸一化
data$median=apply(data[,-1],1,median)
data=data[order(data$GeneName,data$median,decreasing = T),]
data=data[!duplicated(data$GeneName),]
rownames(data)=data$GeneName
uni_matrix <- data[,grep('\\d+',colnames(data))]
uni_matrix <- log2(uni_matrix+1)
colnames(uni_matrix)<- gsub('X','',gsub('\\.','\\-',colnames(uni_matrix)))
uni_matrix<- uni_matrix[,order(colnames(uni_matrix))]
}
save(uni_matrix,pheno,file = 'uni_matrix.Rdata')
2.進行ssGSEA分析
只是用到了處理后的矩陣和基因集兩個內(nèi)容；對score結(jié)果歸一化后進行熱圖繪制幌墓。
獲取免疫細胞的metagenes基因集,得到一個名為mmc3.xlsx的文件但壮，刪除前兩行，保存為mmc3.csv常侣。
##加載包
{
library(genefilter)
library(GSVA)
library(Biobase)
library(stringr)
}
##載入數(shù)據(jù)
load('uni_matrix.Rdata')
gene_set<-read.csv("mmc3.csv")[, 1:2]
head(gene_set)
list<- split(as.matrix(gene_set)[,1], gene_set[,2])
gsva_matrix<- gsva(as.matrix(uni_matrix), list,method='ssgsea',kcdf='Gaussian',abs.ranking=TRUE)
library(pheatmap)
gsva_matrix1<- t(scale(t(gsva_matrix)))
gsva_matrix1[gsva_matrix1< -2] <- -2
gsva_matrix1[gsva_matrix1>2] <- 2
anti_tumor <- c('Activated CD4 T cell', 'Activated CD8 T cell', 'Central memory CD4 T cell', 'Central memory CD8 T cell', 'Effector memeory CD4 T cell', 'Effector memeory CD8 T cell', 'Type 1 T helper cell', 'Type 17 T helper cell', 'Activated dendritic cell', 'CD56bright natural killer cell', 'Natural killer cell', 'Natural killer T cell')
pro_tumor <- c('Regulatory T cell', 'Type 2 T helper cell', 'CD56dim natural killer cell', 'Immature dendritic cell', 'Macrophage', 'MDSC', 'Neutrophil', 'Plasmacytoid dendritic cell')
anti<- gsub('^ ','',rownames(gsva_matrix1))%in%anti_tumor
pro<- gsub('^ ','',rownames(gsva_matrix1))%in%pro_tumor
non <- !(anti|pro)
gsva_matrix1<- rbind(gsva_matrix1[anti,],gsva_matrix1[pro,],gsva_matrix1[non,])
normalization<-function(x){
return((x-min(x))/(max(x)-min(x)))}
nor_gsva_matrix1 <- normalization(gsva_matrix1)
annotation_col = data.frame(patient=pheno$num2)
rownames(annotation_col)<-colnames(uni_matrix)
bk = unique(c(seq(0,1, length=100)))
pheatmap(nor_gsva_matrix1,
show_colnames = F,
cluster_rows = F,cluster_cols = F,
annotation_col = annotation_col,
breaks=bk,cellwidth=5,cellheight=5,
fontsize=5,gaps_row = c(12,20),
filename = 'ssgsea.pdf',width = 8)
save(gsva_matrix,gsva_matrix1,pheno,file = 'score.Rdata')
3.計算score加和后蜡饵，ggplot2進行繪圖
rm(list=ls())
anti_tumor <- c('Activated CD4 T cell', 'Activated CD8 T cell', 'Central memory CD4 T cell', 'Central memory CD8 T cell', 'Effector memeory CD4 T cell', 'Effector memeory CD8 T cell', 'Type 1 T helper cell', 'Type 17 T helper cell', 'Activated dendritic cell', 'CD56bright natural killer cell', 'Natural killer cell', 'Natural killer T cell')
pro_tumor <- c('Regulatory T cell', 'Type 2 T helper cell', 'CD56dim natural killer cell', 'Immature dendritic cell', 'Macrophage', 'MDSC', 'Neutrophil', 'Plasmacytoid dendritic cell')
load('score.Rdata')
anti<- as.data.frame(gsva_matrix1[gsub('^ ','',rownames(gsva_matrix1))%in%anti_tumor,])
pro<- as.data.frame(gsva_matrix1[gsub('^ ','',rownames(gsva_matrix1))%in%pro_tumor,])
anti_n<- apply(anti,2,sum)
pro_n<- apply(pro,2,sum)
patient <- pheno$num2[match(colnames(gsva_matrix1),pheno$num1)]
library(ggplot2)
data <- data.frame(anti=anti_n,pro=pro_n,patient=patient)
anti_pro<- cor.test(anti_n,pro_n,method='pearson')
gg<- ggplot(data,aes(x = anti, y = pro),color=patient) +
xlim(-20,15)+ylim(-15,10)+
labs(x="Anti-tumor immunity", y="Pro-tumor suppression") +
geom_point(aes(color=patient),size=3)+geom_smooth(method='lm')+
annotate("text", x = -5, y =7.5,label=paste0('R=',round(anti_pro$estimate,4),'\n','p<0.001'))
ggsave(gg,filename = 'cor.pdf', height = 6, width = 8)
基本一致，但是細節(jié)之處需要調(diào)節(jié)胳施，不當之處請指正溯祸！
參考文獻：
- 純R代碼實現(xiàn)ssGSEA算法評估腫瘤免疫浸潤程度
- Pan-cancer Immunogenomic Analyses Reveal Genotype-Immunophenotype Relationships and Predictors of Response to Checkpoint Blockade
- Local mutational diversity drives intratumoral immune heterogeneity in non-small cell lung cancer
生信技能樹公益視頻合輯：學習順序是linux，r舞肆，軟件安裝焦辅，geo，小技巧椿胯，ngs組學筷登！
B站鏈接
YouTube鏈接
生信工程師入門最佳指南
學徒培養(yǎng)
