很多情况下，我们进行生信分析都会有一个感兴趣的点，可以是免疫、m6A甲基化、铁死亡、铜死亡等等，但有很多小伙伴不知道到哪里去寻找并下载目标基因集。本期推文，小碗不仅会给大家介绍如何从GSEA检索并下载感兴趣的基因集，而且还会教大家如何从新版TCGA中下载突变数据，并使用maftools包绘制目标基因集的突变瀑布图，以进一步研究这些基因对疾病发生发展的可能的机制。
01 从GSEA下载目标基因集
(1) 进入GSEA，选择MSigDB数据库
长按扫描左侧二维码即可进入GSEA
网址：https://www.gsea-msigdb.org/
(2) 点击目标物种下的“search”，这里以人类为例
(3) 输入关键词，比如铁死亡（ferroptosis）
(4) 这里只检索出来一个包含64个基因的铁死亡相关基因集
(5) 点击下载相应格式，推荐下载gmt格式，后续可以使用GSEABase包的getGMT()直接读入R环境中
02 使用maftools包绘制瀑布图
(1) 设置工作环境
rm(list=ls())
gc()
setwd("../data")
(2) 安装包，加载包
if (!require("BiocManager"))
install.packages("BiocManager")
BiocManager::install("maftools")
BiocManager::install("GSEABase")
BiocManager::install("TCGAbiolinks")
library(maftools)
library(GSEABase)
library(TCGAbiolinks)
(3) 导入目标基因集
gset <-getGmt("WP_FERROPTOSIS.v2023.1.Hs.gmt")
genes <- gset[["WP_FERROPTOSIS"]]@geneIds
(4) 下载突变数据
query <- GDCquery(project = "TCGA-LIHC", #以TCGA肝癌为例
data.category = "Simple Nucleotide Variation",
data.type = "Masked Somatic Mutation",
access = "open")
GDCdownload(query)
GDCprepare(query,
save = TRUE,
save.filename = "../rdata/TCGA-LIHC_SNP.rdata")
load("../rdata/TCGA-LIHC_SNP.rdata")
class(data)
dim(data)
可以看到，下载的突变数据是一个有45440行，141列的数据框。
每一行对应一个突变，每一列为该突变的详细信息，包括突变类型，发生该突变的样本等等。
(5) 导入样本注释信息：
oad("../rdata/样本注释信息.rdata")
class(anno)
dim(anno)
anno[1:4,]
样本注释信息也是一个数据框，对每个突变所在的样本进行注释，与之前的突变数据相对应，突变数据有45440行，注释信息就相应的有45440行。注释信息可以是性别、年龄、基因表达量、样本分组等等，可以同时包含多种注释信息。
(6) 转换为maf格式文件
laml <- read.maf(maf=data,clinicalData = anno)
(7) 绘制瀑布图
pdf("../figure/瀑布图.pdf",width = 8,height = 8)
par(oma=c(2,3,3,3))
oncoplot(maf = laml,
clinicalFeatures = "group", #若无样本注释信息则删去此行，若有多种注释信息则使用c()连接
sortByAnnotation = TRUE, #若无样本注释信息则删去此行
genes = genes,
fontSize = 0.4)
dev.off()
到此，一张发表级的瀑布图就画好了，图片中间以灰色为底色的部分就是瀑布图的主体了，每一行代表目标基因集的一个基因，每一列代表一个样本，每个有色方块表示其对应的基因在其对应的样本发生了突变，突变类型以方块颜色不同而异分别对应图片最下方的图例，例如绿色对应错义突变、红色对应无义突变。图片最上方的柱形图表示的是每个样本的突变负荷（tumor mutation burden, TMB）。图片最右侧的柱形图表示的是每个基因的突变频数。样本注释信息则显示在图片主体的正下方。
