学习笔记总结于『生信技能树』马拉松课程
本文学习复现《Oxidative Stress Response Biomarkers of Ovarian Cancer Based on Single-Cell and Bulk RNA Sequencing》（基于单细胞和Bulk转录组的卵巢癌氧化应激反应生物标志物）一文中的图，其中Oxidative Stress Response缩写为ROS。文章包含了芯片、转录组、单细胞等技术，适合用来复现及学习
八、突变频谱、免疫浸润、gsva
1.maftools
1.1突变频谱图
rm(list = ls())
load("../6.model/lasso_model.Rdata")
load("../6.model/rsurv.Rdata")
# TCGAmutations包整合了TCGA中全部样本的maf文件，所以数据就来源这个R包，tcga_load()会给出一个maf文件
# devtools::install_github(repo = "PoisonAlien/TCGAmutations")
library(TCGAmutations)
library(tidyverse)
maf = tcga_load("OV")
简单查看maf，其中data为主要，clinical.data是临床信息，其余为补充信息
再来看一下主体data的内容，如图84所示
dim(maf@data)
# [1] 35985 21
# 35985行代表着35985条突变，至于是多少人、多少基因，需要另做统计
length(unique(maf@data$Tumor_Sample_Barcode))
# [1] 411 # 病人数
length(unique(maf@data$Hugo_Symbol))
# [1] 12861 # 基因数
# lasso回归模型里面有367个病人，这些病人的ID是16位的，而图84中的突变病人ID是完整的
# 为了能将两个数据连接联动起来，需要把后者也截成16位
samp = data.frame(ID = str_sub(unique(maf@data$Tumor_Sample_Barcode),1,16),
long = unique(maf@data$Tumor_Sample_Barcode))
samp= merge(rsurv,samp,by = "ID") # 一一对应后只剩249人
突变频谱图（瀑布图）的含义
1.2 TMB肿瘤突变负荷
TMB：衡量每个病人的突变情况是否严重，数字越大，突变越多
identical(tmb$Tumor_Sample_Barcode,samp$long)
# [1] FALSE
# 发现顺序不同
samp = merge(samp,tmb,by.x = "long",by.y = "Tumor_Sample_Barcode")
samp$tmb_group = ifelse(samp$total_perMB_log>median(samp$total_perMB_log),"high","low")
library(ggplot2)
library(ggpubr)
ggplot(samp,aes(x = total_perMB_log,y = RS))+
geom_point(aes(color = group))+
geom_smooth(method = "lm")+
theme_bw()+
stat_cor(method = "pearson")
相关系数接近于0，差不多不相关了
load("../6.model/TCGA-OV_sur_model.Rdata")
library(survival)
library(survminer)
sfit <- survfit(Surv(time, event)~tmb_group, data=samp)
ggsurvplot(sfit,
palette = "jco",
risk.table =TRUE,
pval =TRUE,
conf.int =TRUE)
免疫浸润：根据基因表达量推断每一种细胞的组成比例和相对丰度是多少
目前主流的免疫浸润计算方法是CIBERSORT和ssgsea，这里介绍CIBERSORT
2.1对输入数据的要求
下面这段话摘自CIBERSORT的介绍
> 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》
非常清楚的写出了输入数据的要求：
1.不可以有负值和缺失值
2.不要取log
3.如果是芯片数据，昂飞芯片使用RMA标准化，Illumina 的Beadchip 和Agilent的单色芯片，用limma处理。
4.如果是RNA-Seq表达量，使用FPKM和TPM都很合适。
GEO常规的表达矩阵都是这样得到的，直接下载使用即可。注意有的表达矩阵下载下来就已经取过log，需要逆转回去。有的经过了标准化或者有负值，需要处理原始数据，前面写过介绍文了
2.2做成cibersort要求的输入文件
当年这个算法并没有被写成R包，而是只有一个放着函数的脚本CIBERSORT.R，把它下载下来放在工作目录即可。现在这个R包写出来了，但还是要求提供一个txt文件，详见带代码
rm(list = ls())
library(tidyverse)
load("../6.model/TCGA-OV_sur_model.Rdata")#生存分析的数据
load("../6.model/lasso_model.Rdata")#模型
load("../6.model/rsurv.Rdata")
# 需要两个输入文件:
# 一个是表达矩阵文件
# 一个是官网提供的LM22.txt，记录了22种免疫细胞的基因表达特征数据，现在这个文件已经内置在R包里面了
# 要求提供一个txt文件
library(tidyverse)
exp = exprSet#完整的表达矩阵，而不是lasso回归后的基因
exp2 = as.data.frame(exp)
exp2 = rownames_to_column(exp2)#由于CIBERSORT.R读取文件的代码比较粗暴，为了适应它，导出文件之前需要把行名变成一列。不然后面就会有报错。
write.table(exp2,file = "exp.txt",row.names = F,quote = F,sep = "\t")
2.3运行CIBERSORT
# 这段代码大致需要半小时才能运行完毕，不过已经设置了限速步骤了
f = "ciber_OV.Rdata"
if(!file.exists(f)){
#devtools:: install_github ("Moonerss/CIBERSORT")
library(CIBERSORT)
lm22f = system.file("extdata", "LM22.txt", package = "CIBERSORT") #调取R包的内置数据LM22.txt
TME.results = cibersort(lm22f,
"exp.txt" ,
perm = 1000,
QN = T)
save(TME.results,file = f)
}
load(f)
TME.results[1:4,1:4]
# TME.results即运行的结果，一行是一个样本，一列是一种细胞，这些数字即这个细胞在该样本中的相对丰度
# 计算方法：每一种细胞是有自己的marker gene的，根据这些基因的表达量推断，所以才要求我们的输入数据中有表达矩阵
re <- TME.results[,-(23:25)]# 前22列是计算结果，后三列是总结，这三列一般用不上
2.4经典的免疫细胞丰度热图
library(pheatmap)
k <- apply(re,2,function(x) {sum(x == 0) < nrow(TME.results)/2})#在一半以上样本里丰度为0的免疫细胞，就不展示在热图里了
table(k)
# k
# FALSE TRUE
# 10 12
re2 <- as.data.frame(t(re[,k]))
Group = rsurv$group
table(Group)
# Group
# low high
# 184 183
an = data.frame(group = Group,
row.names = colnames(exp))
pheatmap(re2,scale = "row",
show_colnames = F,
cluster_cols = F,
annotation_col = an,
color = colorRampPalette(c("navy", "white", "firebrick3"))(50))
得到如图90所示，不过这个就比较凌乱了，聚类与否都是乱的。可以尝试用分组聚类的方法把图变好看点，不过没必要，既然这么乱就不用放进文章了
2.5直方图
# 可以展示出每个样本的免疫细胞比例
library(RColorBrewer)
# mypalette <- colorRampPalette(brewer.pal(8,"Set1"))#设置配色，不过颜色不好看，换一个
?draw_boxplot()
grcolor = c("#2874C5", "#f87669", "#e6b707", "#868686", "#66C2A5", "#FC8D62", "#8DA0CB",
"#E78AC3", "#A6D854", "#FFD92F", "#E5C494", "#B3B3B3")
# 调整格式：宽变长
dat <- re %>%
as.data.frame() %>%
rownames_to_column("Sample") %>%
mutate(group = Group) %>%
gather(key = Cell_type,value = Proportion,-Sample,-group) %>%
arrange(group)
dat$Sample = factor(dat$Sample,ordered = T,levels = unique(dat$Sample)) #为了固定横坐标的顺序，以因子的时候把细胞分出顺序，这样在画图的时候能一一对应
# 先把group排序，然后将sample设为了因子，确定排序后的顺序为水平，所以两图的顺序是对应的
dat2 = data.frame(a = 1:ncol(exp),
b = 1,
group = sort(Group)) #这个dat2是为了画注释条
# p2才是核心
p1 = ggplot(dat2,aes(x = a, y = b)) +
geom_tile(aes(fill = group)) +
scale_fill_manual(values = grcolor) + # mypalette(22)[1:length(unique(Group))]
theme(panel.grid = element_blank(),
panel.background = element_blank(),
axis.line = element_blank(),
axis.ticks = element_blank(),
axis.text = element_blank(),
axis.title = element_blank()) +
scale_x_continuous(expand = c(0, 0)) +
labs(fill = "Group")
p2 = 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()
) +
scale_y_continuous(expand = c(0.01,0)) +
scale_fill_manual(values = mypalette(22))
library(patchwork)
p1 / p2 + plot_layout(heights = c(1,10),guides = "collect" ) &
theme(legend.position = "bottom")
得到图91，为了看高风险组和低风险组的病人的细胞组成是否存在差别，但肉眼感觉看不出来
2.6箱线图
# 展示免疫细胞之间的比较
ggplot(dat,aes(Cell_type,Proportion,fill = Cell_type)) +
geom_boxplot() +
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))
得到图92，不同的免疫细胞的丰度比较，有些完全是0
下面把每种免疫细胞不同分组的结果做个比较
# 画出全部在组间差异显著的细胞
# 全是0的行去掉
k = colSums(re)>0;table(k)
re = re[,k]
library(tinyarray)
draw_boxplot(t(re),factor(Group),
drop = F,#color = mypalette(length(unique(Group)))
)+
labs(x = "Cell Type", y = "Estimated Proportion")
得到图93
单画某一个感兴趣的免疫细胞
dat2 = dat[dat$Cell_type=="Dendritic cells activated",]
library(ggpubr)
ggplot(dat2,aes(group,Proportion,fill = group)) +
geom_boxplot() +
theme_bw() +
labs(x = "Group", y = "Estimated Proportion") +
theme(legend.position = "top") +
scale_fill_manual(values = grcolor)+ #mypalette(length(unique(Group)))
stat_compare_means(aes(group = group,label = ..p.signif..),method = "kruskal.test")
文章中的图如图95，是为了看免疫检查点基因和构建模型使用的基因的相关性
3.1免疫检查点基因
# 免疫检查点基因参考自文献 https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7786136
IF: 4.7 Q2 /
rm(list = ls())
library(tidyverse)
tmp = readLines("immune_checkpoint.txt")
genes = str_split(tmp,", ")[[1]]#将tmp的内容拆分，这些基因如果之间从表达矩阵里面找，有7个是找不到的。有可能是这7个有别名，只能手动找了，别名也是可以用的
load("../6.model/TCGA-OV_sur_model.Rdata")
k = genes %in% rownames(exprSet) ;table(k)#不过也有41个了，偷点懒不找了
genes = genes[k]
genes
load("../6.model/lassogene.Rdata")
library(Hmisc)
exp = exprSet
da = t(exp[c(genes,lassoGene),])#把目标基因和lasso回归基因拼到一起，并转置。计算基因之间的相关性要把基因放到列名上。如果不放，列名是样本，那就变成样本的相关性了
# Hmisc包的rcorr函数能算R值和P值
m = rcorr(da)$r[1:length(genes),(ncol(da)-length(lassoGene)+1):ncol(da)]#前41行，后26列，右上角
p = rcorr(da)$P[1:length(genes),(ncol(da)-length(lassoGene)+1):ncol(da)]
p[1:3,1:4]
# 倒数第26列：ncol(da)基因的总数67个 - length(lassoGene)26个 + 1
library(dplyr)
tmp = matrix(case_when(as.logical(p<0.01)~"**",
as.logical(p<0.05)~"*",
T~""),nrow = nrow(p))
library(pheatmap)
if (as.numeric(grDevices::dev.cur()) != 1) grDevices::graphics.off() #为了关闭画板
png("ff.png",width = 1200,height = 900)
pheatmap(t(m),
display_numbers =t(tmp),
angle_col =45,
color = colorRampPalette(c("#2fa1dd", "white", "#f87669"))(100),
border_color = "white",
width = 7,
height=9.1,
treeheight_col = 0,
treeheight_row = 0)
dev.off()
得到如图96，文章是竖着画的，我们是横着画的
不过注意色带，虽然很蓝，但其实数值也只到0.2左右。可以试着对pheatmap()函数调整颜色范围
3.2 gsva
gsva是为了做样本（或分组）之间的富集程度比较
rm(list = ls())
load("../6.model/TCGA-OV_sur_model.Rdata")
load("../6.model/lasso_model.Rdata")
library(GSVA)# 和gsea一样，需要基因和通路之间的对应关系
library(msigdbr)# 这个对应关系可以从msigdbr里面找，也可以从网站上下载
KEGG_df = msigdbr(species = "Homo sapiens",category = "C2",subcategory = "CP:KEGG") %>%
dplyr::select(gs_exact_source,gs_name,gene_symbol)
head(KEGG_df)
H_df = msigdbr(species = "Homo sapiens",category = "H") %>%
dplyr::select(gene_symbol,gs_exact_source,gs_subcat,gs_name)
dim(H_df)
#把数据框拆分成列表
kegg_list = split(KEGG_df$gene_symbol,KEGG_df$gs_exact_source)
lapply(kegg_list[1:3], head)
H_list = split(H_df$gene_symbol,H_df$gs_name)
lapply(H_list[1:3], head)
f = "gsva.Rdata"
if(!file.exists(f)){
KEGG_ES <- gsva(expr=exprSet,
gset.idx.list=kegg_list,
parallel.sz=10)# 注意设置10线程对内存有要求，实在不够，设置为1，让它慢慢跑吧//gsva只需要表达矩阵和基因通路对应关系列表就能完成分析
H_ES <- gsva(expr=exprSet,
gset.idx.list=H_list,
parallel.sz=10)
save(KEGG_ES,H_ES,file = f)
}
load(f)#载入得到的结果H_ES，行名是通路，列名是样本，数字是富集分数
# 有了这个结果，可以用limma做差异分析
load("../6.model/rsurv.Rdata")
library(limma)
design = model.matrix(~rsurv$group)
fit = lmFit(KEGG_ES, design)
fit = eBayes(fit)
DEG_KEGG = topTable(fit, coef = 2, number = Inf)
head(DEG_KEGG)
fit2 = lmFit(H_ES, design)
fit2 = eBayes(fit2)
DEG_H = topTable(fit2, coef = 2, number = Inf)
head(DEG_H)
# 做热图，那图个方便，就拿前20个来做
top20_kegg = head(rownames(DEG_KEGG)[order(DEG_KEGG$logFC,decreasing = T)],20)
top20_go = head(rownames(DEG_H)[order(DEG_H$logFC,decreasing = T)],20)
n = KEGG_ES[top20_kegg,]
n = n[,order(rsurv$RS)]
rownames(n) = KEGG_df$gs_name[match(rownames(n),KEGG_df$gs_exact_source)]
g = sort(rsurv$group)
library(tinyarray)
draw_heatmap(n,g,show_rownames = T,cluster_cols = F,scale = F,legend = T,annotation_legend = T)
n2 = H_ES[top20_go,]
n2 = n2[,order(rsurv$RS)]
rownames(n2) = H_df$gs_name[match(rownames(n2),H_df$gs_name)]
draw_heatmap(n2,g,show_rownames = T,cluster_cols = F,scale = F,legend = T,annotation_legend = T)
文章还有一张热图，为了显示差异富集的 HALLMARK 通路的富集评分与 RS 之间的相关性，我们也来做一下
library(corrplot)
M = cor(cbind(t(n),riskscore = rsurv$RS))
p.mat <- cor.mtest(cbind(t(n),riskscore = rsurv$RS))$p
library(paletteer)
my_color = rev(paletteer_d("RColorBrewer::RdYlBu")[-1])
my_color = colorRampPalette(my_color)(10)
png("corrplot.png",height = 1000,width = 1000)
corrplot(M, type="lower",
order="hclust",
col = my_color,
p.mat = p.mat,
sig.level = 0.01,
insig = "blank",
tl.col = "black",
tl.pos = 'l',
cl.pos = 'r')
dev.off()
得到如图99，不过看起来都不怎么相关
谢谢观看
