基因组突变信息的circos图

生信菜鸟团

发布于 2020-07-02 11:32:56

3.8K1

发布于 2020-07-02 11:32:56

写在前面

目前，将基因组多种突变信息如 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(0, 0), 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(0, 0), lty = 2, col = "#00000040")
                    })

circos.clear()