(转帖)-circos-RCircos

原文地址:https://mp.weixin.qq.com/s/QMYRoYktlu0bUTsUCsM4Fg


install.packages('RCircos')

library(RCircos)

如果你安装好了,就直接加载它们即可

library(RCircos)

最基本的circos图

最基本的图,就是直接利用物种的染色体信息,展示出来即可

#导入内建人类染色体数据

data(UCSC.HG38.Human.CytoBandIdeogram)

cyto.info<-UCSC.HG38.Human.CytoBandIdeogram

RCircos.Set.Core.Components(cyto.info,chr.exclude=NULL,tracks.inside=10,tracks.outside=0)

##

## RCircos.Core.Components initialized.

## Type ?RCircos.Reset.Plot.Parameters to see how to modify the core components.

#chr.exclude <- NULL;           设置不显示的染色体,如 c(1,3)

#cyto.info <- UCSC.HG19.Human.CytoBandIdeogram; 设置染色体数据

#tracks.inside <- 10;    设置内部track 个数

#tracks.outside <- 0;    设置外部track 个数

RCircos.Set.Plot.Area()RCircos.Chromosome.Ideogram.Plot()

#绘制染色体表意文字,染色体的默认方法亮点和染色体名称。

复杂化

在前面最基本的图的基础上面把circos复杂化,一步步添加到我们想要的信息。 ### 基础染色体信息 同上

#导入内建人类染色体数据a<-function(...){

data(UCSC.HG19.Human.CytoBandIdeogram);

head(UCSC.HG19.Human.CytoBandIdeogram);

## 这里换了个参考基因组版本,请注意chr.exclude<-NULL;#设置不显示的染色体,如 c(1,3)          cyto.info<-UCSC.HG19.Human.CytoBandIdeogram;#设置染色体数据tracks.inside<-10;#设置内部track 个数tracks.outside<-0;#设置外部track 个数#导入上面四个基本参数RCircos.Set.Core.Components(cyto.info,chr.exclude,tracks.inside,tracks.outside);

RCircos.Set.Plot.Area()

RCircos.Chromosome.Ideogram.Plot()}a()

#绘制染色体表意文字,染色体的默认方法亮点和染色体名称。

在基因组上面根据坐标标记基因

前面已经介绍过这个RCircos.Gene.Label.Data数据是什么了

b<-function(...){a()#Gene Labels and connectors on RCircos Plot#RCircos.Gene.Connector.Plot   绘制染色体表意文字和基因名称之间的连接#RCircos.Gene.Name.Plot 在数据轨道上绘制基因名称data(RCircos.Gene.Label.Data);

name.col<-4;

side<-"in";

track.num<-1;RCircos.Gene.Connector.Plot(RCircos.Gene.Label.Data,track.num,side);

track.num<-2;

RCircos.Gene.Name.Plot(RCircos.Gene.Label.Data,name.col,track.num,side);}

b()


用热图来标记基因组每个区域的数据

一般是拷贝数变异等数据,用RCircos.Heatmap.Plot函数绘制一个数据轨道的热图

data(RCircos.Heatmap.Data);head(RCircos.Heatmap.Data);

##   Chromosome chromStart chromEnd GeneName  X786.O     A498 A549.ATCC

## 1       chr1     934341   935552     HES4 6.75781  7.38773   6.47890

## 2       chr1     948846   949919    ISG15 7.56297 10.49590   5.89893

## 3       chr1    1138887  1142089 TNFRSF18 4.69775  4.55593   4.38970

## 4       chr1    1270657  1284492     DVL1 7.76886  7.52194   6.87125

## 5       chr1    1288070  1293915    MXRA8 4.49805  4.72032   4.62207

## 6       chr1    1592938  1624243 SLC35E2B 8.73104  8.10229   8.36599

##      ACHN   BT.549  CAKI.1

## 1 6.05517  8.85062 7.00307

## 2 7.58095 12.08470 7.81459

## 3 4.50064  4.47525 4.47721

## 4 7.03517  7.65386 7.69733

## 5 4.58575  5.66389 4.93499

## 6 9.04116  9.24175 9.89727

c<-function(...){b()data.col<-6;track.num<-5;side<-"in";RCircos.Heatmap.Plot(RCircos.Heatmap.Data,data.col,track.num,side);}c()


加上散点图

用RCircos.Scatter.Plot函数来添加一个数据轨迹的扫描图

data(RCircos.Scatter.Data);head(RCircos.Scatter.Data);

##   chromosome   start     stop num.mark seg.mean

## 1          1   61735   228706       18  -0.4459

## 2          1  228729   356443       10   0.5624

## 3          1  356542   564621        4  -0.9035

## 4          1  603590  1704138      227   0.3545

## 5          1 1709023  1711414        6   1.2565

## 6          1 1714558 12862252     6276   0.4027

d<-function(...){c()RCircos.Scatter.Data$chromosome=paste0('chr',RCircos.Scatter.Data$chromosome)head(RCircos.Scatter.Data);data.col<-5;track.num<-6;side<-"in";by.fold<-1;RCircos.Scatter.Plot(RCircos.Scatter.Data,data.col,track.num,side,by.fold);}d()

> 这个数据里面的染色体号码起初并没有加上’chr’,我做了个修改


加上直线

用RCircos.Line.Plot函数绘制线为一个数据轨道

data(RCircos.Line.Data);head(RCircos.Line.Data);

##   chromosome    start     stop num.mark seg.mean

## 1          1    61735 16895627     8732   0.1797

## 2          1 16896821 17212714      105  -0.2117

## 3          1 17214822 25574471     5321   0.1751

## 4          1 25574707 25662212       37   0.5064

## 5          1 25663310 30741496     2400   0.1384

## 6          1 30741656 30745210        3  -1.4742

e<-function(...){d()RCircos.Line.Data$chromosome=paste0('chr',RCircos.Line.Data$chromosome)head(RCircos.Line.Data);data.col<-5;track.num<-7;side<-"in";RCircos.Line.Plot(RCircos.Line.Data,data.col,track.num,side);}e()


加上直方图

这里用RCircos.Histogram.Plot 函数绘制一个数据轨迹的直方图

data(RCircos.Histogram.Data);head(RCircos.Histogram.Data);

##   Chromosome chromStart chromEnd     Data

## 1       chr1   45000000 49999999 0.070859

## 2       chr1   55000000 59999999 0.300460

## 3       chr1   60000000 64999999 0.125421

## 4       chr1   70000000 74999999 0.158156

## 5       chr1   75000000 79999999 0.163540

## 6       chr1   80000000 84999999 0.342921

f<-function(...){e()data.col<-4;track.num<-8;side<-"in";RCircos.Histogram.Plot(RCircos.Histogram.Data,data.col,track.num,side);}f()


给每个区间加上title

老实说,我觉得这一个圈圈没有必要绘

data(RCircos.Tile.Data);head(RCircos.Tile.Data)

##   Chromosome chromStart  chromEnd

## 1       chr1          0  23900000

## 2       chr1   12700000  44100000

## 3       chr1   28000000  68900000

## 4       chr1   59000000  94700000

## 5       chr1   99700000 120600000

## 6       chr1  147000000 234700000

g<-function(...){f()track.num<-9;side<-"in";RCircos.Tile.Plot(RCircos.Tile.Data,track.num,side);}g()


最后把有连接关系的基因连线

用RCircos.Link.Plot函数绘制两个或多个基因组位置之间的链接线,请自行查看RCircos.Link.Data数据是什么,如何映射到这个circos图上的。

data(RCircos.Link.Data);head(RCircos.Link.Data);

##   Chromosome chromStart  chromEnd Chromosome.1 chromStart.1 chromEnd.1

## 1       chr1    8284703   8285399         chr1      8285752    8286389

## 2       chr1   85980143  85980624         chr7    123161313  123161687

## 3       chr1  118069850 118070319         chr1    118070329  118070689

## 4       chr1  167077258 167077658         chr1    169764630  169764965

## 5       chr1  171671272 171671550         chr1    179790879  179791292

## 6       chr1  174333479 174333875         chr6    101861516  101861840

h<-function(...){g()track.num<-11;RCircos.Link.Plot(RCircos.Link.Data,track.num,TRUE);data(RCircos.Ribbon.Data);RCircos.Ribbon.Plot(ribbon.data=RCircos.Ribbon.Data,track.num=11,by.chromosome=FALSE,twist=FALSE)}h()


保存为PDF

这个需要具体设置,因为circos图占画布跟其它图不太一样

out.file<-"RCircosDemoHumanGenome.pdf";pdf(file=out.file,height=8,width=8,compress=TRUE);RCircos.Set.Plot.Area();par(mai=c(0.25,0.25,0.25,0.25));plot.new();plot.window(c(-2.5,2.5),c(-2.5,2.5));## 把circos图画在这里即可dev.off();

最后给一个sessionInfo()吧

sessionInfo()

## R version 3.3.3 (2017-03-06)

## Platform: x86_64-apple-darwin13.4.0 (64-bit)

## Running under: macOS Sierra 10.12.4

##

## locale:

## [1] C

##

## attached base packages:

## [1] stats     graphics  grDevices utils     datasets  methods   base    

##

## other attached packages:

## [1] RCircos_1.2.0

##

## loaded via a namespace (and not attached):

##  [1] backports_1.0.5 magrittr_1.5    rprojroot_1.2   tools_3.3.3    

##  [5] htmltools_0.3.5 yaml_2.1.14     Rcpp_0.12.10    stringi_1.1.3  

##  [9] rmarkdown_1.4   knitr_1.15.1    stringr_1.2.0   digest_0.6.12  

## [13] evaluate_0.10

看看内置的测试数据

circos

本来是一个perl程序,专门用来genomic features在全基因组的分布,连接的。 所以输入数据必然要包含染色体,起始终止坐标的。

然后可以是一列value用来调色或者高度信息,可以加入基因名字等各种信息。 如果有两个genomic features的坐标信息,就是连接图啦。

下面是几个例子:

data(RCircos.Gene.Label.Data);head(RCircos.Gene.Label.Data);

##   Chromosome chromStart chromEnd   Gene

## 1       chr1    8921418  8934967   ENO1

## 2       chr1   17345375 17380514   SDHB

## 3       chr1   27022894 27107247 ARID1A

## 4       chr1   41976121 42501596 HIVEP3

## 5       chr1   43803519 43818443    MPL

## 6       chr1   45794977 45805926  MUTYH

data(RCircos.Histogram.Data)head(RCircos.Histogram.Data)

##   Chromosome chromStart chromEnd     Data

## 1       chr1   45000000 49999999 0.070859

## 2       chr1   55000000 59999999 0.300460

## 3       chr1   60000000 64999999 0.125421

## 4       chr1   70000000 74999999 0.158156

## 5       chr1   75000000 79999999 0.163540

## 6       chr1   80000000 84999999 0.342921

data(RCircos.Heatmap.Data)head(RCircos.Heatmap.Data)

##   Chromosome chromStart chromEnd GeneName  X786.O     A498 A549.ATCC

## 1       chr1     934341   935552     HES4 6.75781  7.38773   6.47890

## 2       chr1     948846   949919    ISG15 7.56297 10.49590   5.89893

## 3       chr1    1138887  1142089 TNFRSF18 4.69775  4.55593   4.38970

## 4       chr1    1270657  1284492     DVL1 7.76886  7.52194   6.87125

## 5       chr1    1288070  1293915    MXRA8 4.49805  4.72032   4.62207

## 6       chr1    1592938  1624243 SLC35E2B 8.73104  8.10229   8.36599

##      ACHN   BT.549  CAKI.1

## 1 6.05517  8.85062 7.00307

## 2 7.58095 12.08470 7.81459

## 3 4.50064  4.47525 4.47721

## 4 7.03517  7.65386 7.69733

## 5 4.58575  5.66389 4.93499

## 6 9.04116  9.24175 9.89727

data(RCircos.Link.Data)head(RCircos.Link.Data)

##   Chromosome chromStart  chromEnd Chromosome.1 chromStart.1 chromEnd.1

## 1       chr1    8284703   8285399         chr1      8285752    8286389

## 2       chr1   85980143  85980624         chr7    123161313  123161687

## 3       chr1  118069850 118070319         chr1    118070329  118070689

## 4       chr1  167077258 167077658         chr1    169764630  169764965

## 5       chr1  171671272 171671550         chr1    179790879  179791292

## 6       chr1  174333479 174333875         chr6    101861516  101861840

熟悉物种的染色体信息

就是很简单染色体的每个片段的起始终止坐标信息,取决于物种还有它对应的参考基因组版本

data(UCSC.HG19.Human.CytoBandIdeogram);head(UCSC.HG19.Human.CytoBandIdeogram);

##   Chromosome ChromStart ChromEnd   Band  Stain

## 1       chr1          0  2300000 p36.33   gneg

## 2       chr1    2300000  5400000 p36.32 gpos25

## 3       chr1    5400000  7200000 p36.31   gneg

## 4       chr1    7200000  9200000 p36.23 gpos25

## 5       chr1    9200000 12700000 p36.22   gneg

## 6       chr1   12700000 16200000 p36.21 gpos50

data(UCSC.Mouse.GRCm38.CytoBandIdeogram);

head(UCSC.Mouse.GRCm38.CytoBandIdeogram);

##   Chromosome ChromStart ChromEnd Band   Stain

## 1       chr1          0  8840440  qA1 gpos100

## 2       chr1    8840440 12278390  qA2    gneg

## 3       chr1   12278390 20136559  qA3  gpos33

## 4       chr1   20136559 22101102  qA4    gneg

## 5       chr1   22101102 30941543  qA5 gpos100

## 6       chr1   30941543 43219933   qB    gneg

data(UCSC.Baylor.3.4.Rat.cytoBandIdeogram);

head(UCSC.Baylor.3.4.Rat.cytoBandIdeogram);

##   Chromosome ChromStart ChromEnd Band Stain

## 1       chr1          0 10142096  p13  gneg

## 2       chr1   10142096 24272657  p12  gvar

## 3       chr1   24272657 38517175  p11  gneg

## 4       chr1   38517175 48659271  q11  gpos

## 5       chr1   48659271 69741157  q12  gneg

## 6       chr1   69741157 90025350  q21  gpos

参数列表如下;

rcircos.params<-RCircos.Get.Plot.Parameters();

rcircos.cyto<-RCircos.Get.Plot.Ideogram();

rcircos.position<-RCircos.Get.Plot.Positions();

RCircos.List.Plot.Parameters()

## Parameters for current RCircos session.

##

## Parameters in inch:

## ==============================

## radius.len:      1.84

## chr.ideo.pos:        1.94

## highlight.pos:       2.09

## chr.name.pos:        2.14

## plot.radius:     2.64

## track.in.start:      1.89

## track.out.start: 2.49

## chrom.width:     0.1

## track.padding:       0.02

## track.height:        0.1

##

## Parameters in chromosome unit:

## ==============================

## base.per.unit:       30000

## chrom.paddings:      300

## heatmap.width:       100

## hist.width:      100

## gene name char. width:   500

##

## General R graphic parameters:

## ==============================

## text.size:       0.4

## highlight.width: 2

## point.type:      .

## point.size:      1

## text.color:      black

## heatmap.color:       BlueWhiteRed

## hist.color:      red

## line.color:      black

## scatter.color:       black

## tile.color:      black

## track.background:    wheat

## grid.line.color: gray

## Bezier.point:        1000

## max.layers:      5

## sub.tracks:      5

##

## Data track numbers:

## ==============================

## tracks.inside:       10

## tracks.outside:      0

## Following are procedures to change RCircos plot parameters:

## params <- RCircos.Get.Plot.Parameters();

## params$radius.len <- 2.0;

## params$base.per.unit <- 5000;

## RCircos.Reset.Plot.Parameters(params)

##

## Chromosome ideogram data were automatically modified.

可以看出一个circos图由上面3个部分组成,如下:

RCircos cytoband data

RCircos plot positions

RCircos plot parameters(最复杂啦)

每个部分都可以单独调整细节:

params<-RCircos.Get.Plot.Parameters();

#获得参数列表

params$radius.len<-2.0;

#更改其中的参数

params$base.per.unit<-5000;

#更改其中的参数

RCircos.Reset.Plot.Parameters(params)#参数重置

快分享给你的小伙伴吧,免得他/她妈妈担心

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 204,684评论 6 478
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 87,143评论 2 381
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 151,214评论 0 337
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 54,788评论 1 277
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 63,796评论 5 368
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,665评论 1 281
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 38,027评论 3 399
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,679评论 0 258
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 41,346评论 1 299
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,664评论 2 321
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,766评论 1 331
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,412评论 4 321
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 39,015评论 3 307
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,974评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,203评论 1 260
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 45,073评论 2 350
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,501评论 2 343

推荐阅读更多精彩内容