|
作者：解螺旋.麦子 如需转载请注明来源：解螺旋·医生科研助手
上回我们简单介绍了R语言作Circos图的基本数据形式和方法，也谈到圈圈套圈圈意味着函数套函数、甚至是数据表的复杂组合。今天我们就从数据类形进一步探索怎么个套法。复习请进→_→《作Circos图不用学Perl了，我们有R！》 circlize包的一大特点就是可以用高级函数把低级函数组织起来，形成圈圈布局的图，每个圈、每个板块里都是我们曾经很熟悉的点图、线形图、柱形图等等。像circos. genomicInitialize()、circos.genomicTrackPlotRegion()这样用于组织布局的就是高级函数，像circos.genomicPoints()、circos.genomicLines()、circos.genomicRect()这样用来在每个圈、每个板块中填充内容的就是低级函数。嵌套的基本思路是这样的： circos.genomicInitialize(数据表名称) 数据表名称, ylim=c(n1,n2), panel.fun = function(region,value,…) { circos.genomicPoints( region, value, pch = 16, cex = 0.5) 其中，最关键的就是circos.genomicTrackPlotRegion和panel.fun，处理好它俩就能升一个境界。它们在不同的数据组织形式里的具体应用也有所不同，下面我们就按数据的组织形式来扒一扒。 circos.genomicTrackPlotRegion用来生成新的圈；panel.fun=function(region,value,…)则是一个自定义的函数，Region是板块的起始点和终止点，value则是描述Y轴上的值。在它后面的{}中便可嵌入各种低级函数。 ylim()括号中的两个数，就是Y轴的范围，方向是从圆心向外，即，假如我们数据的范围是3~57，那么设ylim=c(0,60)，则靠近圆心为y轴的0，外层为60，以此安排点的分布。X轴默认为顺时针方向。 Circos中，数据的基本形式是由3个列组成的表格，第1列决定一个圈分为哪几个板块，第2列是板块的起点，第3列是止点，这就叫BED格式，只能画一个原始的圈。我们上回演示的简单数据和原始圈，包括第一个使用circos.genomicTrackPlotRegion来生成的颜色标识圈： 要把数据搞大，有几种方式，一是增加行，一是增加列，还有一种是把好几个数据框（data frame）组成一个列表（list）。显然，能让图形变得更复杂更有趣的是后两种。 一、单个数据表 在原来的BED表中加入其他列，以数值形式呈现，比如找到的CNV、SNP数量等等，这样我们作图的素材还是单个数据表。 在R中运行circos.initializeWithIdeogram()会默认调用UCSC基因库中hg19的cytoband文件，在Excel里观察它的结构是这个样子的： 前三列就是基本格式，后面虽然还有两列，但显然不是“数值型”。生成的原始圈是以它为框架的，后面的染色标识则生成了第一圈，图形长这样： 要往里面添加内容，随机生成一个Bed来看看： bed=generateRandomBed(nc=1)，nc是指除了前三列之外的数值型列数。 它的第一列和原始圈那张表的第一列有对应的内容。后面的value1就是随具体课题而定的数值型变量。 在这样的单一数据框中，上面说的关键函数之一，panele.fun=function(region,value,…)里的region显然就是指start和end两列，value就是指多出来的那些数值型的列了，就是这里的value1。于是可以嵌个低级函数了，就用散点图演示吧： circos.genomicTrackPlotRegion( bed, track.height=0.1,ylim = c(-1,1), panel.fun = function(region, value, ...) { circos.genomicPoints( region, value, pch = 16, cex = 0.1,...) }) track.height是指圈的高度，以原始圈半径的比例来计算，此处0.1就是半径的10%；Pch和cex是点的样式，在RStudio的Help里可以搜到具体的样式代号。 二、数据列表 如果把一组数据表组成一个列表（list），panel.fun则迭代用于各个数据表。再随机生成一个list来观察它的结构： bedlist = list( bed01=generateRandomBed(nc=1,nr=200), bed02=generateRandomBed(nc=1,nr=200) ) 这里用nc和nr控制一下列和行的数量，不要像上面那样1万多行做出个密恐图。nc是指除了前三列之外的数值型列数；但nr就比较随缘，最后生成的行数可能并不是我们设的200，对此，circlize包的作者也只是提醒，并没有解释。就这样吧，反正除了做练习，应该也不大会用到随机生成数据的方法。 就这样，bed01和bed02是互有联系的两个数据框，它们的第1列也有相同的内容。还可以观察到两个表每行的起止点并不一样，这是因为起止点不是对第1列负责，而是对应后面的数值。然后再生成点图： circos.genomicTrackPlotRegion( bedlist,track.height=0.1,ylim = c(-1,1), panel.fun = function(region, value, ...) { i = getI(...) circos.genomicPoints( region, value,pch = 16, cex = 0.1,col = i, ...) }) 这段代码里，panel.fun的{}中，比单一数据表那里多了一个i=getI(...)，circos.genomicPoints里也相应地多了一个col = i，这俩是什么意思呢？任性的Help是这么说的： 所以就是“憋说话，加我”的意思嘛╮(╯▽╰)╭ 其实也可以理解为，一个list中有多个数据表，getI可以用来找到当前所用的表格中对应的数值型列。去掉这个不明觉厉的getI和col = i，生成下图中最里面那圈，图里就少了一组数据了。 如果表里有多个数值列，就可以用numeric.column来控制想要输出哪一列，像这段代码： circos.genomicTrackPlotRegion( bedlist,track.height=0.1,ylim = c(-1,1), numeric.column = c(4, 5), panel.fun = function(region, value, ...) { circos.genomicPoints(region, value,pch = 16, cex = 0.1, ...) }) 就比上面多了一个numeric.column = c(4, 5)，是指第一个表中的第4列和第二个表中的第5列。也就是说，每生成一个圈，都只能用到各个表中的1个列。 本例的低级函数只用了circos.genomicPoints()，但其他的circos.genomicLines()、circos.genomicRect()都大同小异，大家可以自行探索。还是多多练习circos.genomicTrackPlotRegion和panel.fun的用法吧。
|
|
|
