基因组突变信息的circos图

写在前面

目前,将基因组多种突变信息如 SNV / INDEL 和 CNV 一起呈现在基因组上的可视化方式很多,比较受欢迎的就是以 CIRCOS 的形式来展示。有一个软件就叫 CIRCOS ,是perl语言写的,使用起来比较麻烦,然后在生信技能树也有介绍一个R包RCircos。

这里我们推荐用顾祖光老师的 R 包 circlize,circlize包在德国癌症中心的华人博士Zuguang Gu开发的。

安装R包

安装R包比较简单,但是如何使用这个R包,需要学习一下帮助文档。在文末给出的链接中,作者已经详细的描述了这个包的绘图思想和原理,以及各个参数的作用。在这里,我先用作者给出的测试数据,来学习一下如何将 SNV / INDEL 和 CNV 的数据与基因组展示在同一张 circos 图中。

首选是安装 R 包

options("repos" = c(CRAN="http://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
if (! require(circlize)) {
  install.packages('circlize')
}
library(circlize)

然后就是获取测试数据,在这个R包中,绘图的数据都是以 bed 格式,前三列一般是坐标,后面可以增加任意列作为注释。作者非常人性化的设计了一个函数 generateRandomBed 用于获取示例的 bed 数据

bed = generateRandomBed(nr = 200)
head(bed)
#    chr    start      end     value1
# 1 chr1  6255701 16825950 -0.7966037
# 2 chr1 23545777 26681342 -1.4405458
# 3 chr1 41443165 46719497  0.1430059
# 4 chr1 52909447 54931298 -0.1015327
# 5 chr1 65273586 74868674 -0.5058086
# 6 chr1 75096106 79815973  0.16868

然后第一步就是先绘制基因组,也是一个函数就可以实现,这里的基因组选择的是 hg38

circos.initializeWithIdeogram(species = "hg38")

第二步就是添加 SNV /INDEL 的信息,以点图展示,基于上面生成的 bed 示例数据:

circos.genomicTrack(bed, 
                    panel.fun = function(region, value, ...) {
                      circos.genomicPoints(region, value, pch = 16, cex = 0.5, col = 1...)

第三步自然就是 CNV 的信息啦,因为这里只是学习,所以就用同一个bed,最后要用一个函数结束绘图circos.clear()

circos.genomicTrack(bed, 
                    panel.fun = function(region, value, ...) {
                      circos.genomicRect(region, value, ytop.column = 1, ybottom = 0, 
                                         col = ifelse(value[[1]] > 0"red""green"), ...)
                      circos.lines(CELL_META$cell.xlim, c(00), lty = 2, col = "#00000040")
                    })
circos.clear()

数据准备

上面简单学习了一下绘图的方法,接下来我们就要对实际的数据进行绘图了(以 case1_biorep_A_techrep 样本为例)。首先要将 SNV 和 CNV 的数据处理成 bed 格式

somatic mutations

maf = read.table('./7.anotation/vep/case1_biorep_A_techrep_vep.maf', comment.char = "#",header = T,sep = 't',quote = "")
choose.col = c("Chromosome","Start_Position","End_Position","Hugo_Symbol","t_depth","t_ref_count","t_alt_count")
somatic = maf[ ,choose.col]
somatic$vaf = somatic$t_alt_count/somatic$t_depth
somatic.bed = somatic[,c(1,2,3,8)]

CNV

seg = read.table('./8.cnv/gatk/segment/case1_biorep_A_techrep.cr.igv.seg', header=T)
seg$Segment_Mean[ seg$Segment_Mean < -1 ] = -1
seg.bed = seg[,c(2,3,4,6)]

可视化

circos.initializeWithIdeogram(species = "hg38")
circos.genomicTrack(somatic.bed, 
                    panel.fun = function(region, value, ...) {
                      circos.genomicPoints(region, value, pch = 16, cex = 0.5, col = 1...)
                    })
circos.genomicTrack(seg.bed, 
                    panel.fun = function(region, value, ...) {
                      circos.genomicRect(region, value, ytop.column = 1, ybottom = 0, 
                                         col = ifelse(value[[1]] > 0"red""blue"), ...)
                      circos.lines(CELL_META$cell.xlim, c(00), lty = 2, col = "#00000040")
                    })
circos.clear()

基本上到这里就可以看出来SNV/INDEL和CNV的分布情况,可以进一步标注出体细胞突变位点所在的基因,也可以进一步美化,具体要看R包的帮助文档了。顾老师写的这个 circlize R包的功能非常强大,感兴趣的朋友可以深入了解。

参考链接

[1] https://academic.oup.com/bioinformatics/article/30/19/2811/2422259

[2] https://jokergoo.github.io/circlize_book/book/

生物信息学

如何安装github上的R包

2020-6-29 13:31:24

生物信息学

GWAS候选位点与候选基因的筛选思路

2020-6-29 18:02:19

0 条回复 A文章作者 M管理员
    暂无讨论,说说你的看法吧
个人中心
购物车
优惠劵
今日签到
有新私信 私信列表
搜索