果叔今天给大家介绍一个CNV可视化与焦点分析的R包——GenomeTornadoPlot。它2022年3月28日发表于Bioinformatics期刊,快和果叔一起看看吧!
简介
在癌症研究中,焦点(focal)拷贝数变异(CNVs)非常重要,因为它们能够定位驱动基因。更具体地说,由于对癌细胞的选择性压力,致癌基因和肿瘤抑制基因更容易受到这些事件的影响,而邻近的乘客则相对较少受到影响。在基因组位点中存在多个候选基因共同居住的情况下,需要进行仔细的比较,以确定多基因最小缺失区域或真正的单一驱动基因。对于大型癌症基因组队列中的焦点CNVs的研究需要专门的可视化和统计分析工具。
GenomeTornadoPlot能够从队列测序数据生成CNV类型、位置和长度的基因中心可视化。此外,它还支持对近邻基因进行成对比较,以识别共突变模式或驱动基因-乘客基因层次结构。GenomeTornadoPlot提供的视觉检查进一步得到了适应性的局部和全局焦点评分的支持。集成到GenomeTornadoPlot R包中的是全面的PCAWG CNV数据库,包括来自Pan-cancer Whole Genome Analysis项目46个队列的2976个癌症基因组实体。GenomeTornadoPlot R包可用于基于PCAWG数据进行探索性或基于假设的分析,也可以与用户提供的数据相结合使用。
为什么要有这个包?
焦点拷贝数变异可以精确定位驱动基因。更具体地说,原癌基因和肿瘤抑制基因比邻近的乘客基因更经常受到这些事件的影响。此外,多基因最小缺失区域的检测使得能够检测协同共突变。
使用 GenomeTornadoPlot 包,我们可以可以:
•在队列水平上将交替选择的 CNVs 与选择的基因重叠可视化为tornado plot。
•观察同一染色体中两个选定基因的 CNVs,并比较这些事件
•用不同的方法计算焦点得分。
焦点度评分算法
标准焦点评分
一般来说,我们假设具有相对更多焦点事件的基因比具有更多广泛事件的基因具有更高的得分。在这里,我们通过以下方式定义标准焦点分数:
其中 M 是焦点变化事件的总数,Lmax是最长焦点变化事件的长度。
边得分
为了消除相邻基因的影响,我们实现了另一种算法,并称之为“边得分”。它被定义为:
其中近邻 1 和近邻 2 是目标基因的近邻基因,如果目标基因位于染色体的边缘,则仅有的近邻基因被计为近邻 1 和近邻 2。
简而言之,边得分是一个基因与其相邻基因的标准焦点得分之间的平均差异。
用户定义的分数
用户还可以根据 CNV 数据集定义自己的焦点评分方法。可以将分数导入到可视化。
方便用户选择是计算标准焦点值还是边缘焦点值 socre 作为“焦点值”。请注意,每个基因的焦点值是根据我们给出的数据计算出来的,因此在不同的数据集中并不是一个恒定值。
跑代码时卡顿、电脑不给力让人抓狂!找果叔试用稳定高速的服务器,让分析顺畅无比!
代码学不会?bug 频繁出现,束手无策?实操生信分析课程赶快学起来!滴滴果叔领取体验课程哦~
线上课程教学
课题设计、定制生信分析
云服务器租赁
加微信备注99领取使用
安装
在安装 GenomeTornadoPlot 之前,我们还要安装其他依赖项:
dependencies.packages = c('ggplot2', 'data.table', 'devtools','grid', 'gridExtra','tiff',"shiny","shinydashboard","entropy")
install.packages(dependencies.packages)
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install(c('GenomicRanges','quantsmooth','IRanges'))
工作流程
安装
1.在Git repository中,单击“Clone or Download”。
2.复制链接。
3.打开终端输入:
git clone https://github.com/chenhong-dkfz/GenomeTornadoPlot
1.打开文件夹 GenomeTornadoPlot 并在 RStudio 中打开“GenomeTornadoPlot.rproj ”文件。
2.在 RStudio 控制台中,键入:
devtools::install()
快速入门
第 0 步:准备数据
首先,你可以准备一个类似于bed的数据,并将其导入到 R Session 中。在 R 中,它应该是一个数据帧,如下所示:
我们可以从https://github.com/chenhong-dkfz/GenomeTornadoPlot-files)下载真实的 PCAWG 数据进行测试
load(file = "path to the file/chr17.RData")
chr17 <- cnv_chr
names(chr17)[4] <- "Score"
names(chr17)[6] <- "Cohort"
names(chr17)[7] <- "PID"
knitr::kable(head(chr17, 10))
该列得分记录每个 CNV 事件的拷贝数。
确保数据框的列名非常重要。请注意每列名称的第一个大写字母。
准备好数据后,我们可以将 PLOT 应用于数据
步骤 1:计算聚焦度分数
运行 MakeData() 函数:
library(GenomeTornadoPlot)
input_gene_1="TP53"
sdt <- MakeData(CNV=chr17,gene_name_1 = input_gene_1,score.type = "del")
这CNV是我们刚刚导入的一个类似于 Bed 的 Data.Frame.
Focality_score <- sdt@gene_score
这里将计算聚焦度分数。
其他参数定义如下:
•gene_name_1: 第一个基因的名称。
•gene_name_2: 第二个基因的名称(可选)。
•score.type: 如果值为"del",则计算删除的集中度得分。如果值为"dup",则计算复制的集中度得分。
•score.method: 如果值为"normal",则计算标准的集中度得分。如果值为"edge",则计算边缘得分。
•gene_score_1: 如果用户给定了该值,则将此输入值用作可视化中第一个基因的集中度得分。
•gene_score_2: 如果用户给定了该值,则将此输入值用作可视化中第二个基因的集中度得分(可选)。
这里SDT是包含所选基因的 CNV 信息的 R 对象。且它应该是步骤 2 输入。
步骤 2:画图
运行 TornadoPlots() 函数:
plotlist1 <- TornadoPlots(sdt,color.method="ploidy",sort.method="length",multi_panel=FALSE)
以下是每个参数的解释:
•data: 由MakeData()函数生成的R对象。
•legend: 可以设置为"pie"(默认)或"barplot"(可选)。
•color: CNV颜色的向量,可选。
•color.method: 如何给CNV着色。可以是"cohort"(默认)或"ploidy"(可选)。
•sort.method: 如何对CNV进行排序。可以是"length"(默认),"cohort"或"ploidy"(可选)。
•SaveAsObject: 如果为TRUE,则返回一个rastergrob对象。如果为FALSE,则仅保存绘图。
•format: 如果SaveAsObject为FALSE,包将在磁盘上的路径中保存绘图。如果此值为"tiff",则将绘图保存为tiff图像。如果此值为"eps",则将绘图保存为EPS矢量图像。
•path: 如果SaveAsObject为FALSE,包将在磁盘上的路径中保存绘图。图像将保存在路径中。
•multi_panel: 如果为TRUE,将显示多个面板图。
•zoomed: 值应为"global","region"或"gene"。它表示图形将如何放大。
•orient: 如果值为"v",则垂直排列的图将被显示。如果值为"h",则水平排列的图将被显示。
•drop.low.amp: 如果值为TRUE,则不会在图中显示CN<5的扩增。
•font.size.factor: 图中显示的字体大小的重新缩放。
这个query描述的是一个流程,其中plotlist1是一个包含绘图输出的列表。
在第一步中,如果你只给出gene_name_1,那么在完成第二步之后,你将得到一个标准tornado图和一个针对这个基因的“dup_del”图。否则,如果你也给出了gene_name_2,那么你将得到一个“twin”图和一个“mixed”图。
GenomeTornadoPlot Easy2Use(ShinyApp)
为了帮助用户方便地生成GenomeTornadoPlot , GenomeTornadoPlot 软件包中提供了一个shinyapp。
runExample()
用户可以从本地计算机上传基于文本的 CNV 文件,设置参数并在 Shiny 应用程序中下载生成的图。(我们可以从https://github.com/chenhong-dkfz/GenomeTornadoPlot-files/FOXP1_RYBP.txt)下载用于测试的示例输入文件
例子
在这里,你可以简单地使用以下代码来制作tornado图。虚拟数据附加在包中。第一个例子是单个基因。
data("cnv_STK38L", package = "GenomeTornadoPlot")
data_genea <- MakeData(CNV=cnv_STK38L,gene_name_1 = "STK38L",score.type="del")
plotlist1 <- TornadoPlots(data_genea,gene.name="STK38L",color.method="ploidy",
sort.method="length",SaveAsObject=TRUE,multi_panel=FALSE,orient="v")
如果你需要的只是焦点分数,只需使用以下命令:
data_genea@gene_score
下面的代码允许用户对基因列表的分数进行排序,并且他们可以很容易地生成图。我们使用 CHR21 中的几个基因作为例子。
g.list <- c("ERG","TMPRSS2","RUNX1")
x=data.frame(row.names=g.list,score=rep(1,length(g.list)))
for(gene in g.list){
print(gene)
sdt <- MakeData(CNV=chr21,gene_name_1 = gene,score.type = "del")
print(sdt@gene_score)
x[gene,]=sdt@gene_score
}
print(x)
然后,用户可以将前 N 个基因放入 for 循环中,以执行 TornadoPlots 函数。
我们也可以更进一步,试着打印一个标准的基因组龙卷风图:
grid.arrange(plotlist1[[1]])
彩色线条代表 CNV 事件。在图中,我们可以很容易地找到它们在染色体中的位置。饼图代表事件的组群贡献。本例中的颜色代表队列。但用户也可以更改参数,并使复制数量或长度的颜色。图表下方的分数是该基因的“焦点分数”。
在某些情况下,一个基因在不同的队列中扮演着不同的角色。deletion/duplication图有助于识别这种情况。
grid.arrange(plotlist1[[2]])
在这里,感兴趣的基因在大多数队列中是重复的,而在其他一些队列中缺失更频繁。
下面是另一个基因 PTEN 的例子。
load(file = "path to file/chr10.RData")
chr10 <- cnv_chr
names(chr10)[4] <- "Score"
names(chr10)[6] <- "Cohort"
names(chr10)[7] <- "PID"
input_gene_1="PTEN"
sdt <- MakeData(CNV=chr10,gene_name_1 = input_gene_1,score.type = "del", gene_score_1 = 30)
print(sdt@gene_score)
这里,我们将参数gene_score_1设置为 30,并且我们将看到聚焦性分数将为 30。
sdt <- MakeData(CNV=chr10,gene_name_1 = input_gene_1,score.type = "del")
print(sdt@gene_score)
如果我们不为gene_score_1设置任何值,则函数将使用默认的“正常”方法计算焦点分数。
然后,我们用不同的排序和颜色方法绘制龙卷风图。
plotlist1 <- TornadoPlots(sdt,color.method="cohort",sort.method="cohort",multi_panel=FALSE,
font.size.factor=1.5,orient="v")
grid.arrange(plotlist1[[1]])
按组群排序并按队列着色。
plotlist1 <- TornadoPlots(sdt,color.method="ploidy",sort.method="cohort",multi_panel=FALSE,
font.size.factor=1.5,orient="v")
grid.arrange(plotlist1[[1]])
按队列分类,按倍性分类。
plotlist1 <- TornadoPlots(sdt,color.method="ploidy",sort.method="length",multi_panel=FALSE,
font.size.factor=1.5,orient="v")
grid.arrange(plotlist1[[1]])
按长度排序,按倍性排序。
我们可以设置图形的方向为"h"(水平排列)并将缩放设置为TRUE,以便垂直地排列CNVs。
plotlist1 <- TornadoPlots(sdt,color.method="ploidy",sort.method="length",multi_panel=FALSE,
font.size.factor=1.5, orient="h",zoomed=TRUE)
grid.arrange(plotlist1[[1]])
grid.arrange(plotlist1[[1]])
在这种情况下,我们能够查看事件的详细信息。这里的虚线表示基因 PTEN 的起始和终止位置。
我们还可以通过将multi_panel设置为ture来绘制多面板图。详细的基因/变异统计数据、不同级别的龙卷风图和放大图将显示在一个复杂的图中。
TornadoPlots(sdt,color.method="ploidy",sort.method="length",multi_panel=TRUE,
font.size.factor=1.5, orient="v",zoomed=TRUE)
我们也可以在基因对上用 GenomeTornadoPlot 。
load(file = "path to file/chr21.RData")
chr21 <- cnv_chr
names(chr21)[4] <- "Score"
names(chr21)[6] <- "Cohort"
names(chr21)[7] <- "PID"
input_gene_1="ERG"
input_gene_2="TMPRSS2"
sdt <- MakeData(CNV=chr21,gene_name_1 = input_gene_1,gene_name_2=input_gene_2,score.type="del")
plot_twin <- TornadoPlots(sdt,sort.method="cohort",SaveAsObject=T,multi_panel=FALSE)
画两个图:
grid.arrange(plot_twin[[1]])
此外,混合图显示了单独基因 1、单独基因 2 或两个基因重叠的 CNVs 的比例。
grid.arrange(plot_twin[[2]])
总结
果叔今天给大家介绍了一个名为GenomeTornadoPlot的R包,用来可视化和对多个大型癌症队列中的焦点CNVs进行统计分析。该软件包集成了PCAWG数据库,允许用户在不需要额外数据或整合自生成的小数据集的情况下进行探索性或基于假设的分析,以解决虚拟生物学上的问题。为了便于应用,GenomeTornadoPlot还实现了ShinyApp。
小伙伴们快去试试吧,特别是癌症领域的小伙伴哟~
参考文献
[1] Hong C, Thiele R, Feuerbach L. GenomeTornadoPlot: a novel R package for CNV visualization and focality analysis. Bioinformatics. 2022 Mar 28;38(7):2036-2038. doi: 10.1093/bioinformatics/btac037. PMID: 35099519; PMCID: PMC8963283.
[2]https://github.com/chenhong-dkfz/GenomeTornadoPlot/blob/master/README.md
如果小伙伴有其他数据分析需求,可以尝试使用本公司新开发的生信分析小工具云平台,零代码完成分析,非常方便奥,云平台网址为:(http://www.biocloudservice.com/home.html),其中也包括了通路表达分析(http://www.biocloudservice.com/313/313.php),单细胞的基因共表达分析(http://www.biocloudservice.com/906/906.php)等各种小工具哦~,有兴趣的小伙伴可以登录网站进行了解。
不会分析还想用生信工具助力发文咋办?有这顾虑的朋友,想一步到位就带着想法来,不论是代码实操还是在线文章结果复现,果叔照样能提供,还有大家都想要的服务器,找果叔获取就对了!
往期回顾
01 |
02 |
03 |
04 熬夜=爆肝被实锤了!最新研究:失眠打鼾缺觉党,“困” 境重重小心脂肪肝风险 “爆” 表!UKB+两步MR+中介分析,深度剖析! |