如何用Python绘制Circos图

如何用Python绘制Circos图,第1张

用Python实现Circos图的绘制在线绘制的Circos有一定局限性,如对数据的要求、个性化的局限和处理速度等的问题,但如果你是一个Pythoneer或者喜欢用更加Pythonic的方式来个性化地绘制Circos图,那么今天就跟随我一起用代码实现这一目标吧!

安装Circos包

首先,登录Python的包索引网站Python Package Index(PyPI,正确读音是:Pie Pee Ai),找到Circos包的下载页:

https://pypi.python.org/pypi/Circos/1.3.5

该包/模块的作者是我的好友Eric Ma。你可以选择下载wheeler文件,然后本地安装。也可以在shell下直接通过pip进行安装:

pipinstallcircos

注意,所支持的Python版本必须是3.x,对2不支持。

选择数据

当安装了circos包后,我们就可以直接应用这个包来写代码了。为了演示方便,我需要应用一些数据。作为内科医师,就让我来展示一下老本行:处理药物与肝酶细胞色素P450的相互关系的可视化。由于是为了抛大销缺砖引玉,所以绘制出的Circos图相对简单。

我们先从美国FDA官网下载不同细胞色素相关的各种口服药物表。共202种常用的口服药物,涵盖内科学、肿瘤学、神经科和心理学等学科。数据文件如下:

可以看到这个数据的结构:是按肝细胞色素酶进行分类,共斗模分8个列。这8个细胞色素酶分别是:CYP1A2,CYP2B6, CYP2C8, CYP2C9, CYP2C19, CYP2D6, CYP2E1和CYP3A4。我们将要建立各个口服药与这些肝酶之间关系的Circos图,从而了解通过相同肝酶代谢或转化的药物之间是否存在相互作用。

导入各个模块和读入数据

导入各个模块:

fromcircosimportCircosPlot

importxlrd

importpandasaspd

importnumpyasnp

读入文件:

filename='.\\MedicationInteraction.xlsx'

book=xlrd.open_workbook(filename)

print('File loaded!')

提取数据:

nrows=book.sheet_by_name('Sheet1').nrows

header=book.sheet_by_name('Sheet1').row_values(0)

data=[book.sheet_by_name('Sheet1').row_values(i) fori inrange(1, nrows)]

df=pd.DataFrame(data, columns =header)

df[df==''] = np.nan

读取后,药物和酶的数据为pandas的DataFrame数据结构,细胞色素P450酶的名字为columns的名字。我们可以检查一下数据:

修数据,尤其是处理NA数据

df_dict={}

foriinrange(len(df.columns)):

df_dict[df.icol(i).name] =list(df.icol(i).dropna())

节点和连线

创建节点(nodes)数据,在我这个例子里就是各个药物和肝酶:

nodes=[]

forkeyindf_dict.keys():

nodes.extend(df_dict[key])

nodes=list(nodes)

headers=list(df.columns)

enzymes=['0'] * 5

forheaderinheaders:

enzymes.append(header)

enzymes.extend(['0']*5)

nodes.extend(enzymes)

创建连线(edges)数据,我们应用tuple(元组)这个数据结构来表示药物与特定肝酶之间的关系:

edges_origin=[]

forkeyindf_dict.keys():

forvalue indf_dict[key]:

edges_origin.append((key, value))

绘图

绘制Circos图:

c=CircosPlot(nodes, edges_origin, radius =10,

nodecolor="blue",

edgecolor="red",

)

c.draw()

得到了下面这张所滚辩有药物与肝酶之间的图:

左上方是8个肝脏细胞色素P450酶(CYP1A2、CYP2B6、CYP2C8、CYP2C9、CYP2C19、CYP2D6、CYP2E1和CYP3A4)。其它点即为202种口服药物。每种药物都与参与代谢和转化它的P450酶相连。与相同酶连接的不同药物,理论上应该都存在相互作用,但具体如何还要看与酶的作用机理。

个性化绘图

如果我们打算分别可视化出不同肝酶的关系图形,我们只需改变连线信息,即edges信息:

edges=[]

&#8205forvalueindf_dict['CYP2B6']:

edges.append(('CYP2B6', value))

c=CircosPlot(nodes, edges, radius =10,

nodecolor="orange",

edgecolor="orange",

)

c.draw()

从而我们得到了各种肝酶所代谢和转化药物的图形

用PS将它们合并:

相同肝酶所代谢和转化的药物用相同颜色的edges表示。

显示特定药物

最后,我们可以挑选其中一些感兴趣的药物来进行观察,例如,我从这202个药物中指定几个我感兴趣的药物:

propafenone(心律平), acetaminophen(对乙酰氨基酚), paclitaxel(紫杉醇), ibuprofen(布洛芬), losartan(洛沙坦), omeprazole(奥美拉唑), carvediolo(卡维地洛), codeine(可待因), theophylline(茶碱), quinidine(奎尼丁), verapamil(异搏定), lovastatin(洛伐他汀), nitrendipine(尼群地平)

然后重新建立edges:

medications=['propafenone', 'acetaminophen', 'paclitaxel', 'ibuprofen', 'losartan', 'omeprazole', 'carvedilol', 'codeine', 'theophylline', 'quinidine', 'verapamil', 'lovastatin', 'nitrendipine']

edges_candidate=set()

formedicationinmedications:

foredge inedges_origin:

if medication==edge[1]:

edges_candidate.add(edge)

edges_candidate=list(edges_candidate)

然后再绘图:

c=CircosPlot(nodes, edges_candidate, radius =10,

nodecolor="black",

edgecolor="black",

)

c.draw()

从而得到这张图。

相信大家都听说过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数据上传 、 二代测序数据解读 。


欢迎分享,转载请注明来源:内存溢出

原文地址: http://www.outofmemory.cn/tougao/12152169.html

(0)
打赏 微信扫一扫 微信扫一扫 支付宝扫一扫 支付宝扫一扫
上一篇 2023-05-21
下一篇 2023-05-21

发表评论

登录后才能评论

评论列表(0条)

保存