【R画图】环形热图

Python022

【R画图】环形热图,第1张

热图(heatmap)在生信领域基本就是常规操作,基本技能,入门操作。能画热图的工具也有很多,我自己常用的R包是pheatmap。

最近经常看见环状的热图,所以就搜了一下资料学习一下,测试一下。

环状热图我也经常会在论文中看到,用法和热图相同,但更适合于需要展示较多基因(数据)时来使用。

===安装=====

install.packages("circlize")

install.packages("ComplexHeatmap")

library('ComplexHeatmap')

library('circlize')

===测试===

data<-read.table("expression.txt",header=T,row.names=1)

head(data)

//转化为矩阵并对其进行归一化

madt<-as.matrix(data)

madt2<-t(scale(t(madt)))

Heatmap(madt2)     //这就是默认参数,普通的热图

range(madt2)   //查看值的分布

mycol=colorRamp2(c(-1.7, 0.3, 2.3),c("blue", "white", "red")) //定义颜色范围

circos.heatmap(madt2,col=mycol)  //默认参数

注:

circos.clear()  //绘制完成后需要使用此函数完全清除布局

circos.par(gap.after=c(50))    //

调整圆环首尾间的距离,数值越大,距离越宽

circos.heatmap(madt2,col=mycol, dend.side="inside",rownames.side="outside",rownames.col="black",rownames.cex=0.9,rownames.font=1,cluster=TRUE)

注:dend.side:控制行聚类树的方向,inside为显示在圆环内圈,outside为显示在圆环外圈

#rownames.side:控制矩阵行名的方向,与dend.side相同;但注意二者不能在同一侧,必须一内一外

#cluster=TRUE为对行聚类,cluster=FALSE则不显示聚类

下面是换了inside和outside的结果:

===聚类树的调整和美化(需要用到两个别的包)===

install.packages("dendextend")     //改颜色

install.packages("dendsort")          //聚类树回调

library(dendextend)

library(dendsort)

circos.par(gap.after=c(50))

circos.heatmap(madt2,col=mycol,dend.side="inside",rownames.side="outside",track.height=0.38,rownames.col="black",rownames.cex=0.9,rownames.font=1,cluster=TRUE,dend.track.height=0.18,dend.callback=function(dend,m,si){color_branches(dend,k=15,col=1:15)})

注:track.height:轨道的高度,数值越大圆环越粗  dend.track.height:调整行聚类树的高度  dend.callback:用于聚类树的回调,当需要对聚类树进行重新排序,或者添加颜色时使用

包含的三个参数:dend:当前扇区的树状图;m:当前扇区对应的子矩阵;si:当前扇区的名称   color_branches():修改聚类树颜色

circos.clear()

//添加图例标签等

lg=Legend(title="Exp",col_fun=mycol,direction= c("vertical"))

grid.draw(lg)

//添加列名

circos.track(track.index=get.current.track.index(),panel.fun=function(x,y)

{

if(CELL_META$sector.numeric.index==1)

{

cn=colnames(madt2)

n=length(cn)

circos.text(rep(CELL_META$cell.xlim[2],n)+convert_x(0.3,"mm"), //x坐标

4.1+(1:n)*0.7, //y坐标和间距

cn,cex=0.8,adj=c(0,1),facing="inside")

}

},bg.border=NA)

mycol2=colorRamp2(c(-1.7, 0.5, 2.3),c("#57ab81", "white", "#ff9600"))

也可以更改上面的颜色,从而改变热图的配色:

===分组热图绘制========

#circos.heatmap()内只能是一个矩阵,但如果矩阵数据存在分组,可以用split参数来指定分类变量

split = sample(letters[1:2], 40, replace = TRUE)

split = factor(split, levels = letters[1:2])

circos.clear()

circos.par(gap.after=c(22))

circos.heatmap(madt2,col=mycol,split=split,dend.side="inside",rownames.side="outside",track.height = 0.38,

rownames.col="black",rownames.cex=0.9,

rownames.font=1,

cluster=TRUE,

dend.track.height=0.18,

dend.callback=function(dend,m,si) 

{

color_branches(dend,k=15,col=1:15)

}

)

//假设有两个热图的矩阵数据

madt2<-t(scale(t(madt)))

madt3<-t(scale(t(madt)))

split2 = sample(letters[1:2], 40, replace = TRUE)

split2 = factor(split2, levels = letters[1:2])

circos.clear()

circos.par(gap.after=c(8))

circos.heatmap(madt2,col=mycol2,split=split2,dend.side="outside",

cluster=TRUE,

dend.track.height=0.2,

dend.callback=function(dend,m,si) {

color_branches(dend,k=15,col=1:15)

}

)

circos.heatmap(madt3, col = mycol,rownames.side="inside",rownames.cex=0.8)

相信大家都听说过circos图,但是亲自画过的人可能就很少,这主要因为软件的安装和使用稍微有一点麻烦。其实,circos图也是可以在线绘制的,这样就简单多了!一起来了解一下吧!

在circos官网(http://circos.ca/)的最右方有个“CIRCOS ONLINE”选项,这里可以实现在线绘制部分circos图。

打开后界面如下:

以微生物多样性分析中样品与物种丰度circos图绘制为例,给大家讲解circos图的绘制功能。该图能够很直观的反映各样品中不同物种所占的比例,以及物种在不同分组或者样品中的分布关系。

绘制circos图

1.数据准备

首先我们要做的就是准备画图所用到的数据,所用数据为物种在各样品中的相对丰度,这里只选用丰度大于0.01的物种用于绘图,数据如下(列名A、B、C为样品,行名Acetobacteraceae等是科一水平的物种分类):

OTUABC

Acetobacteraceae0.5063653216696110.5968872412369940.457528142134733

Arcobacteraceae0.0003294904846044670.0179133872520980.000426249200782749

Bacteroidaceae0.01752092807693420.04558718113953470.352221339584988

Dysgonomonadaceae0.001842974249051360.02565003004872960.0330226880824598

Lachnospiraceae0.005691857602178260.01390206286339050.0173870923992018

Lactobacillaceae0.174952205775860.2379461150250890.0588340146862225

Pseudomonadaceae0.00213263621353880.02862956070929480.0127991010016856

Ruminococcaceae0.003124728441908290.005061219761203110.0274388235522058

Sphingomonadaceae0.2578607015612780.007113946230875610.00898610815104722

由于网站要求的数据格式为非负整数,故将所有的数据乘1000(系统会自动截掉小数点后的数据),输入数据则变为:

OTUABC

Acetobacteraceae506.365321669611596.887241236994457.528142134733

Arcobacteraceae0.32949048460446717.9133872520980.426249200782749

Bacteroidaceae17.520928076934245.5871811395347352.221339584988

Dysgonomonadaceae1.8429742490513625.650030048729633.0226880824598

Lachnospiraceae5.6918576021782613.902062863390517.3870923992018

Lactobacillaceae174.95220577586237.94611502508958.8340146862225

Pseudomonadaceae2.132636213538828.629560709294812.7991010016856

Ruminococcaceae3.124728441908295.0612197612031127.4388235522058

Sphingomonadaceae257.8607015612787.113946230875618.98610815104722

2.绘图

数据准备好就可以来绘制circos图了,只需要导入数据就可以。

生成的图片如下:

可以看到,图中的物种和样品完全是按照字母顺序排列的,我们希望物种和样品分别位列两边,这里可以人为的对其指定顺序。方法也很简单,就是在数据的第一行和第一列用数字来指定顺序。如下:

OTUOTU1   2   3

OTUOTUABC

12Acetobacteraceae506.365321669611596.887241236994457.528142134733

10Bacteroidaceae17.520928076934245.5871811395347352.221339584988

8Dysgonomonadaceae1.8429742490513625.650030048729633.0226880824598

6Lachnospiraceae5.6918576021782613.902062863390517.3870923992018

11Lactobacillaceae174.95220577586237.94611502508958.8340146862225

7Pseudomonadaceae2.132636213538828.629560709294812.7991010016856

5Ruminococcaceae3.124728441908295.0612197612031127.4388235522058

9Sphingomonadaceae257.8607015612787.113946230875618.98610815104722

4Arcobacteraceae0.32949048460446717.9133872520980.426249200782749

第一行指定了样品的顺序,而第一列按丰度指定物种的顺序。生成图片时要勾选下图红框中的选项(排序所用),不然会报错哦!

新图如下:

图中由于部分物种丰度较低,导致物种名重叠,解决这个问题可以改变文字的布局。这时就需要进行设置了。

3.图片设置

点击"settings"进入设置界面,会有很多的设置选项,可以对图片进行细调。

这里只需要修改两个地方即可,将下图第一个红框改为“no”,可以调整文字为垂直布局,避免重叠;但是如果物种名太长,又可能会超出图片范围,所以要缩小圆圈的半径,即将第二个红框改为small。

修改并保存设置后,重新生成图片:

好了,今天我就先给大家介绍到这里,希望对您的科研能有所帮助!祝您工作生活顺心快乐!

更多生物信息课程:

1. 文章越来越难发?是你没发现新思路,基因家族分析发2-4分文章简单快速,学习链接: 基因家族分析实操课程 、 基因家族文献思路解读

2. 转录组数据理解不深入?图表看不懂?点击链接学习深入解读数据结果文件,学习链接: 转录组(有参)结果解读 ; 转录组(无参)结果解读

3. 转录组数据深入挖掘技能-WGCNA,提升你的文章档次,学习链接: WGCNA-加权基因共表达网络分析

4. 转录组数据怎么挖掘?学习链接: 转录组标准分析后的数据挖掘 、 转录组文献解读

5.  微生物16S/ITS/18S分析原理及结果解读 、 OTU网络图绘制 、 cytoscape与网络图绘制课程

6. 生物信息入门到精通必修基础课,学习链接: linux系统使用 、 perl入门到精通 、 perl语言高级 、 R语言画图

7. 医学相关数据挖掘课程,不用做实验也能发文章,学习链接: TCGA-差异基因分析 、 GEO芯片数据挖掘 、 GSEA富集分析课程 、 TCGA临床数据生存分析 、 TCGA-转录因子分析 、 TCGA-ceRNA调控网络分析

8.其他课程链接: 二代测序转录组数据自主分析 、 NCBI数据上传 、 二代测序数据解读 。