一 ssGSEA
一) 以之前得到的表达矩阵exp_1为例
安装R包GSVA
if(!require("GSVA")) BiocManager::install("GSVA",update = F,ask = F)
1. 准备细胞marker
来自文献《Pan-cancer Immunogenomic Analyses Reveal Genotype-Immunophenotype Relationships and Predictors of Response to Checkpoint Blockade》 中的 cellMarker
cellmarker <- read.csv("./immune/cellMarker.csv",header = T)
colnames(cellmarker)[2] <- "celltype" ##更改第二列有空格的行名,便于后续操作
此时的cellMarker需要进行转换
cellmarker <- split(as.matrix(cellmarker$Metagene) ,cellmarker$celltype)
###用split函数按不同的细胞类型重新组合
得到一个列表的集合 : 每个免疫细胞名作为列表名,每个子列表中是相应的metagene
2.准备表达量矩阵
使用之前做过的exp_1为例,此时exp_1是行为基因，列为样本的表达矩阵
3. 使用ssGSEA量化免疫浸润
exp_1 <- as.matrix(exp_1)
library(GSVA)
gsva_data <- gsva(exp_1 , cellmarker , method = "ssgsea")
4. 作图
library(pheatmap)
x <- pheatmap(gsvadata,
cluster_rows = T,
annotation_legend = T,
show_rownames = T,
color = colorRampPalette(c("blue","white","red"))(100),
cellwidth = 20 ,cellheight =10,
fontsize=10)
二 用TCGA数据来实战ssGSEA
1 数据下载
1.1TCGA数据
cohort: TCGA Pan-Cancer (PANCAN)
dataset: gene expression RNAseq - TOIL RSEM tpm 中数据
[下载地址]
[注释文件]
1.2 临床信息
来自文章An Integrated TCGA Pan-Cancer Clinical Data Resource to Drive High-Quality Survival Outcome Analytics
[下载地址]
1.3 下载gtf文件
2 准备细胞marker
3 处理数据
3.1 加载数据
dd <- data.table::fread("tcga_RSEM_gene_tpm",data.table = F) ###加载TCGA数据
clin <- data.table::fread("Survival_SupplementalTable_S1_20171025_xena_sp",data.table = F) ####加载临床数据
colnames(clin)[2:3] <- c("TCGA_id","type") ###名字不好看,换了
3.2 获取有临床信息的样本
index <- clin$sample %in% colnames(dd) ###设立索引,找到同时有临床信息的样本
clin <- clin[index,]
dd <- dd[,c("sample",clin$sample)]
colnames(dd)[1] <- "gene_id" ###更换名称,便于后续操作
3.3 探针转换
library(rtracklayer)
gtf1 <- rtracklayer::import('gencode.v29.annotation.gtf') ###读取gtf数据
gtf_df <- as.data.frame(gtf1)
gtf2 <- gtf_df[ , c(10,12)] ###选取带有gene_id &gene_name两列
###此时gtf2和dd有共同的列gene_id
library(dplyr)
library(tidyr)
library(tibble)
exp_1 <- as.data.frame(dd) %>%
merge(gtf2 , by="gene_id") %>%
select(-gene_id)%>%
select(gene_name , everything())%>%
distinct(gene_name , .keep_all = T)%>%
column_to_rownames(colnames(.)[1])
(提示内存不够,换套代码
tcga_panmRNA_expr <- gtf_df %>%
filter(type=="gene",gene_type=="protein_coding") %>% #筛选gene,和编码指标
dplyr::select(c(gene_name,gene_id)) %>% #选取两列
inner_join(dd,by ="gene_id") %>% #和表达量数据合并
dplyr::select(-gene_id) %>% ## 去掉多余的列
filter(gene_name!="NA") %>% ## 去掉基因名称中可能的NA，可有可无
distinct(gene_name,.keep_all = T) %>% ## 去掉重复
column_to_rownames("gene_name") ## 列名转为行名
3 用ssGSEA来量化浸润水平
expr <- tcga_panmRNA_expr
expr <- as.matrix(expr)
library(GSVA)
gsva_data <- gsva(expr,cellMarker, method = "ssgsea") ###此步骤耗时较久
4 添加分组信息
tcga_gsva <- as.data.frame(t(gsva_data)) ##变回数据框
tcga_gsva1 <- cbind(clin,subtype=substring(rownames(tcga_gsva),14,15),tcga_gsva)###将临床数据和浸润水平结合
tcga_panmRNA_expr <- as.data.frame(t(tcga_panmRNA_expr))
tcga_panmRNA_expr <- cbind(clin,subtype=substring(rownames(tcga_panmRNA_expr),14,15),tcga_panmRNA_expr)
####将表达量的数据与临床数据结合
