加了这个图，保底5+SCI！学会位点结构域绘制，让你笑傲江湖！（附代码）

基因突变氨基酸位点Lollipop图

大家好，我是阿琛。在上一期内容中，我们介绍了基因组突变中肿瘤突变负荷（TMB）的计算方法，及其结合转录组信息的简单应用方法（我用10分钟，拆解了肿瘤免疫高分SCI的万用分析，简单又高效！）。然而，TMB主要代表样品中整体的突变水平，其无法具体解释某个特定基因的突变发生在哪个重要的结构域，或者突变产生特定影响。
因此，我们需要针对特定基因，展示某个癌症类型里某个基因的突变位点，进而推测其中某位点的突变会产生的蛋白质功能。对于基因氨基酸位点的展示，maftools包中的lollipopPlot()函数提供了绘制功能，但是其中得到的图形相对比较简单。除此之外，trackViewer 包可以进行氨基酸位点变异位置精美图谱的展示。

下面，我们一起来看一下如何绘制基因突变后氨基酸位点的Lollipop图。

1. TCGA突变数据的下载

对于TCGA数据库中数据的下载，我们使用R的TCGAbiolinks包的GDCquery_Maf()函数来进行下载。

rm(list = ls()) #清空环境变量options(stringsAsFactors = F)library(TCGAbiolinks)luad.varscan2.maf <- GDCquery_Maf("LUAD", pipelines = "varscan2")save(luad.varscan2.maf, file = "LUAD.maf.rda") #保存突变结果

与之前的一样，在此我们还是选择varscan2类型的数据进行下载。

#load("LUAD.maf.rda")mut <- luad.varscan2.maf[luad.varscan2.maf$Hugo_Symbol == "TP53",]write.table(mut, "LUAD.mutation_TP53.txt",row.names = F,quote = F,sep = "\t")

提取突变数据中基因名为TP53的对应突变数据；结果显示，其中存在291个TP53突变信息。

2. 读取氨基酸突变数据

rt <- read.table("LUAD.mutation_TP53.txt",as.is = T,sep = "\t",header = T)rt <- rt[!is.na(rt$HGVSp),c("RefSeq", "HGVSp", "HGVSp_Short")]

读取前面保存的突变数据，通过is.na()函数分析其中基因突变但未引起氨基酸位点改变的样品信息，结合!取反，最终得到基因突变且引起氨基酸位点发生改变的样品，一共273个突变信息。

3. 提取氨基酸位点信息

3.1 读取参考结构域信息

library(maftools)gff = readRDS(file = system.file('extdata', 'protein_domains.RDs', package = 'maftools'))id <- strsplit(x = as.character(rt$RefSeq), split = '.', fixed = TRUE)[[1]][1]protein <- gff[HGNC %in% "TP53"][refseq.ID == id,]

接着，根据maftool包里面提供的蛋白质结构域数据，获取其中对应的TP53基因的结构域。结果显示，在TP53基因中，共存在3种不同的结构域类型，分别是P53 transactivation motif，P53 DNA-binding domain，以及P53 tetramerisation motif。

3.2 分析氨基酸位置

data <-rt$HGVSp_Shortprot.spl = strsplit(x = as.character(data),split = '.', fixed = TRUE)prot.conv = sapply(sapply(prot.spl,function(x) x[length(x)]), '[', 1)pos = gsub(pattern = 'Ter.*', replacement = '',x = prot.conv)pos = gsub(pattern = '[[:alpha:]]', replacement = '', x = pos)pos = gsub(pattern = '\\*$', replacement = '', x = pos)pos = gsub(pattern = '^\\*', replacement = '', x = pos)pos = gsub(pattern = '\\*.*', replacement = '', x = pos)pos = as.numeric(sapply(X = strsplit(x = pos,split = '_', fixed = TRUE),FUN = function(x)x[1])) aa = paste0(unlist(regmatches(data,gregexpr("p[.].[0-9]+", data))),"X")mutpos = data.frame(position = pos,mutation = data,aa = aa,stringsAsFactors = F)

根据参考结构域，结合突变信息中提供的氨基酸突变位点，我们需要将两者结合起来。

mutpos <- mutpos[order(mutpos$pos),]

最后，根据突变位置的大小顺序，对突变结果进行排序整理。

3.3 统计不同氨基酸突变类型得突变频率

pos2 <- mutpos[!duplicated(mutpos),]rownames(pos2) <- pos2$mutationpos2$rate <- 1pos2 <- pos2[names(table(mutpos$mutation)),]pos2$rate <-table(mutpos$mutation)head(pos2)

由于多个样品可能出现同一氨基酸位点的突变情况，因此我们需要对其进行去重复，合并，统计。

table(pos2$rate)

使用table()函数对突变频率进行统计描述，其中同一位点突变的最多存在7个样品。

4. 绘制Lollipop图

4.1 R包准备与读取

#if (!requireNamespace("BiocManager", quietly = TRUE))# install.packages("BiocManager")#BiocManager::install("trackViewer")library(trackViewer) #绘图包library(RColorBrewer) #颜色包

在此，我们使用trackViewer包来绘制氨基酸位点的Lollipop图，同时使用RColorBrewer包来获取图形对应的颜色。

4.2 构建基因特征

features <- GRanges("chr1", IRanges(start=protein$Start, end=protein$End, names=protein$Description))features$height <- 0.07features$fill<-brewer.pal(9,"Set1")[1:3]

对基因特征的起始位置，高度，以及颜色进行设置。

4.3 构建样品特征

sample.gr <- GRanges("chr1", IRanges(pos2$position, width = 1, names = pos2$mutation))sample.gr$label.parameter.rot <- 90#label角度sample.gr$label <- as.character(pos2$rate) #标记sample.gr$label.col <- "white"#标记的颜色sample.gr$color <-sample.gr$border <- brewer.pal(9,"Set1")[4] #点及连线颜色sample.gr$score<- log2(pos2$rate) #点所在高度

根据样品信息，分别设置标签对应的内容，角度，颜色等等信息。