很多情况下，我们进行生信分析都会有一个感兴趣的点，可以是免疫、m6A甲基化、铁死亡、铜死亡等等，但有很多小伙伴不知道到哪里去寻找并下载目标基因集。本期推文，小碗不仅会给大家介绍如何从GSEA检索并下载感兴趣的基因集，而且还会教大家如何从新版TCGA中下载突变数据，并使用maftools包绘制目标基因集的突变瀑布图，以进一步研究这些基因对疾病发生发展的可能的机制。
1.从GSEA下载目标基因集
（1）进入GSEA（GSEA (gsea-msigdb.org)），选择MSigDB数据库
（2）点击目标物种下的“search”，这里以人类为例。
（3）输入关键词，比如铁死亡（ferroptosis）
（4）这里只检索出来一个包含64个基因的铁死亡相关基因集
（5）点击下载相应格式，推荐下载gmt格式，后续可以使用GSEABase包的getGMT()直接读入R环境中。
2.使用maftools包绘制瀑布图
# ①设置工作环境----
rm(list=ls())
gc()
setwd("../data")
# ②安装包，加载包----
if (!require("BiocManager"))
install.packages("BiocManager")
BiocManager::install("maftools")
BiocManager::install("GSEABase")
BiocManager::install("TCGAbiolinks")
library(maftools)
library(GSEABase)
library(TCGAbiolinks)
# ③导入目标基因集
gset <- getGmt("WP_FERROPTOSIS.v2023.1.Hs.gmt")
genes <- gset[["WP_FERROPTOSIS"]]@geneIds
# ④下载突变数据
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列的数据框。
每一行对应一个突变，每一列为该突变的详细信息，包括突变类型，发生该突变的样本等等。
# ⑤导入样本注释信息----
load("../rdata/样本注释信息.rdata")
class(anno)
dim(anno)
anno[1:4,]
样本注释信息也是一个数据框，对每个突变所在的样本进行注释，与之前的突变数据相对应，突变数据有45440行，注释信息就相应的有45440行。注释信息可以是性别、年龄、基因表达量、样本分组等等，可以同时包含多种注释信息。
# ⑥转换为maf格式文件
laml <- read.maf(maf=data,clinicalData = anno)
# ⑦绘制瀑布图
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）。图片最右侧的柱形图表示的是每个基因的突变频数。样本注释信息则显示在图片主体的正下方。
好了，今天的分享就到这，如果对你有帮助就给小碗一个免费的赞吧，我们下期再见！
网友评论