生信图可视化大全基因组数据绘图科研信息学

R语言circos绘制

2019-09-28  本文已影响0人  nnlrl

介绍

circos可以作为数据展示的一种手段，主要用于基因组数据的可视化，比如甲基化，组蛋白修饰，snp，indel，sv，以及常见的基因表达水平的展示，一般的circos都是基于perl语言的，数量繁多的配置文件总是很难按照我们想要的方式展示最终结果。而在R中也存在一些绘制circos的包，本例中主要介绍一下OmicCircos这个包的使用

用法

在使用OmicCircos绘制circos图之前，我们首先也需要准备配置文件，比如染色体长度信息，染色体定位信息等等，如果还有snp，甲基化等信息也可以一并展示，在OmicCircos中主要使用circos函数来绘制circos图形。

rm(list=ls());
library(OmicCircos);
options(stringsAsFactors = FALSE);

#首先导入自带数据进行测试，包括基因表达，cnv等数据
data("TCGA.PAM50_genefu_hg18");
data("TCGA.BC.fus");
data("TCGA.BC.cnv.2k.60");
data("TCGA.BC.gene.exp.2k.60");
data("TCGA.BC.sample60");
data("TCGA.BC_Her2_cnv_exp");

#cnv的p值进行log10转化
pvalue <- -1 * log10(TCGA.BC_Her2_cnv_exp[,5]);
pvalue <- cbind(TCGA.BC_Her2_cnv_exp[,c(1:3)], pvalue);

#筛选Her2亚型
Her2.i <- which(TCGA.BC.sample60[,2] == "Her2");
Her2.n <- TCGA.BC.sample60[Her2.i,1];
Her2.j <- which(colnames(TCGA.BC.cnv.2k.60) %in% Her2.n);

#筛选Her2亚型患者cnv信息
cnv    <- TCGA.BC.cnv.2k.60[,c(1:3,Her2.j)]; 
cnv.m  <- cnv[,c(4:ncol(cnv))];
cnv.m[cnv.m >  2] <- 2;
cnv.m[cnv.m < -2] <- -2;
cnv <- cbind(cnv[,1:3], cnv.m);

#筛选Her2亚型患者的基因表达信息
Her2.j   <- which(colnames(TCGA.BC.gene.exp.2k.60) %in% Her2.n);
gene.exp <- TCGA.BC.gene.exp.2k.60[,c(1:3,Her2.j)]; 

#设定颜色
colors <- rainbow(10, alpha=0.5);

#画图，首先建立一个空白图层
pdf("OmicCircos4vignette10.pdf", 8,8);
par(mar=c(2, 2, 2, 2));
plot(c(1,800), c(1,800), type="n", axes=FALSE, xlab="", ylab="", main="");


#选择需要放大的坐标
zoom <- c(1, 22, 939245.5, 154143883, 0, 180);

#开始画图
circos(R=400, cir="hg18", W=4,   type="chr", print.chr.lab=TRUE, scale=TRUE, zoom=zoom);#最外层染色体信息
circos(R=300, cir="hg18", W=100, mapping=gene.exp, col.v=4,  type="heatmap2", cluster=TRUE, col.bar=TRUE, lwd=0.01, zoom=zoom);#次外层基因表达信息
circos(R=220, cir="hg18", W=80,  mapping=cnv,      col.v=4,   type="ml3", B=FALSE, lwd=1, cutoff=0, zoom=zoom);#第三层cnv信息
circos(R=140, cir="hg18", W=80,  mapping=pvalue,   col.v=4,    type="l",   B=TRUE, lwd=1, col=colors[1], zoom=zoom);#第四层p值曲线
circos(R=130, cir="hg18", W=10,  mapping=TCGA.BC.fus, type="link", lwd=2, zoom=zoom);#基因共表达连线


## 局部放大操作
the.col1=rainbow(10, alpha=0.5)[1];
highlight <- c(140, 400, 11, 282412.5, 11, 133770314.5, the.col1, the.col1);
circos(R=110, cir="hg18", W=40, mapping=highlight, type="hl", lwd=2, zoom=zoom);
the.col2=rainbow(10, alpha=0.5)[6];
highlight <- c(140, 400, 17, 739525, 17, 78385909, the.col2, the.col2);
circos(R=110, cir="hg18", W=40, mapping=highlight, type="hl", lwd=2, zoom=zoom);

## highlight link
highlight.link1 <- c(400, 400, 140, 376.8544, 384.0021, 450, 540.5);
circos(cir="hg18", mapping=highlight.link1, type="highlight.link", col=the.col1, lwd=1);
highlight.link2 <- c(400, 400, 140, 419.1154, 423.3032, 543, 627);
circos(cir="hg18", mapping=highlight.link2, type="highlight.link", col=the.col2, lwd=1);

## zoom in chromosome 11
zoom <- c(11, 11, 282412.5, 133770314.5, 180, 270);
circos(R=400, cir="hg18", W=4,   type="chr", print.chr.lab=TRUE, scale=TRUE, zoom=zoom);
circos(R=300, cir="hg18", W=100, mapping=gene.exp, col.v=4,  type="heatmap2", cluster=TRUE, lwd=0.01, zoom=zoom);
circos(R=220, cir="hg18", W=80,  mapping=cnv,      col.v=4,   type="ml3", B=FALSE, lwd=1, cutoff=0, zoom=zoom);
circos(R=140, cir="hg18", W=80,  mapping=pvalue,   col.v=4,    type="l",   B=TRUE, lwd=1, col=colors[1], zoom=zoom);

## zoom in chromosome 17
zoom <- c(17, 17, 739525, 78385909, 274, 356);
circos(R=400, cir="hg18", W=4,   type="chr", print.chr.lab=TRUE, scale=TRUE, zoom=zoom);
circos(R=300, cir="hg18", W=100, mapping=gene.exp, col.v=4,  type="heatmap2", cluster=TRUE, lwd=0.01, zoom=zoom);
circos(R=220, cir="hg18", W=80,  mapping=cnv,      col.v=4,   type="ml3", B=FALSE, lwd=1, cutoff=0, zoom=zoom);
circos(R=140, cir="hg18", W=80,  mapping=pvalue,   col.v=4,    type="l",   B=TRUE, lwd=1, col=colors[1], zoom=zoom);
dev.off()

需要注意的是，omicCircos是使用type函数来指定不同图形的展示方式的，而type可以指定的图片形式包括以下几种：

上一篇 下一篇

猜你喜欢

热点阅读