当前位置: 首页 > news >正文

pycirclize python包画circos环形图

pycirclize python包画circos环形图

很多小伙伴都有画环形图的需求,网上也有很多画环形图的教程,讲解circos软件和circlize R包的比较多,本文介绍一款python包:pyCirclize。适合喜欢python且希望更灵活作图的小伙伴。

pyCirclize包实际上也是以matplotlib模块为基础进行开发,个人使用体验感觉比circos软件灵活很多,例如circos软件没办法为各圈写标题,图注也比较单一,相比之下pycirclize的图注和对扇区的调整更加灵活。详见官方的教程文档:https://moshi4.github.io/pyCirclize/。

1. 安装pycirclize

直接使用pip工具安装即可,要求Python 3.9以上版本

pip install pycirclize

2. 实例图

多说无益,直接上一个样图。
下图是一个甲基化相关环形图,包含了常用环形图的诸多要素,根据实例图代码修改应该可以满足大部分作图需求了。

  • 第一圈为染色体:需要显示染色体ID和刻度(大刻度标出刻度值,但起始的大刻度不显示,免得首位的刻度值重叠显得很乱),标出小刻度。
  • 第二圈为case组相对于control组的高甲基化位点:CG、CHG、CHH三种颜色均显示,图中由于CHH类型位点太多导致其他两个看不太清了。标注出甲基化水平刻度线,从0到0.9共10条浅灰色的线。
  • 第三圈为基因密度热图:用黑白渐变展示。
  • 第四圈为case组相对于control组的低甲基化位点:同第二圈,但方向相反。
    在这里插入图片描述

3. 作图

俗话说得好:Talk is cheap. Show me the code.

3.1 数据准备

  1. 第一圈的chromosome.bed 文件:
    共三列,染色体名,start,end。
CHR1 0       27139696
CHR2 0       27139696
  1. 第二圈和第四圈的gDMR.hyper.txt
    共4列,染色体ID,Start,End,值。
CHR1 1482   1487   0.167943
  1. 基因密度
    根据gff文件提取
cut -f 1,3 chromosome.bed > test
bedtools makewindows -g test -w 1000000 > 1M
grep -w "gene" my.gff3 |awk '{print $1"\t"$4"\t"$5}'|uniq > gene.pos
bedtools intersect -a 1M -b gene.pos -c >gene.density

3.2 实例代码

python代码见下,细节部分注释了内容。希望能达到抛砖引玉的作用吧,祝大家科研顺利!

from pycirclize import Circos
import matplotlib.pyplot as plt
import matplotlib as mpl
import numpy as np
import pandas as pd
np.random.seed(0)
from matplotlib.patches import Patch
from matplotlib.lines import Line2D# Initialize Circos from BED chromosomes
circos = Circos.initialize_from_bed("chromosome.bed", space=1,start=5,end=355,endspace=False) # 定义染色体,space设置间距circos.text("图中标题,可以设置为组名vs组名", size=12, r=25,weight='bold') # 标题文字,默认位于在图中央
circos.text("Gene density", size=10, r=16,weight='bold')# Plot chromosome name
for sector in circos.sectors:sector.text(sector.name, size=5)# Plot outer trackouter_track = sector.add_track((98, 100))outer_track.axis(fc="lightgrey")major_interval = 10000000 # 设置大刻度minor_interval = 1000000 # 设置小刻度if sector.size > minor_interval:outer_track.xticks_by_interval(major_interval, label_formatter=lambda v: f"{v / 1000000:.0f} Mb" if v != 0 else None,label_size=4)outer_track.xticks_by_interval(minor_interval, tick_length=1, show_label=False)Hyper_CG=pd.read_table("CG_gDMR.hyper.txt",header=None)x = np.array(Hyper_CG[Hyper_CG[0]==sector.name][1]+(Hyper_CG[Hyper_CG[0]==sector.name][2]-Hyper_CG[Hyper_CG[0]==sector.name][1])/2)y = np.array(Hyper_CG[Hyper_CG[0]==sector.name][3])track1 = sector.add_track((75, 95), r_pad_ratio=0.1)# 注释圈名if sector.name == circos.sectors[0].name:circos.text("Hyper", r=track1.r_center, size=8)for r in range(77, 95, 2):sector.line(r=r, ec="lightgrey")track1.scatter(x, y,c="#EECA40",vmin=0,vmax=1,ec="black",lw=0.1,alpha=0.5)Hyper_CHG=pd.read_table("CHG_gDMR.hyper.txt",header=None)x = np.array(Hyper_CHG[Hyper_CHG[0]==sector.name][1]+(Hyper_CHG[Hyper_CHG[0]==sector.name][2]-Hyper_CHG[Hyper_CHG[0]==sector.name][1])/2)y = np.array(Hyper_CHG[Hyper_CHG[0]==sector.name][3])track1.scatter(x, y,c="#FD763F",vmin=0,vmax=1,ec="black",lw=0.1,alpha=0.5)Hyper_CHH=pd.read_table("CHH_gDMR.hyper.txt",header=None)x = np.array(Hyper_CHH[Hyper_CHH[0]==sector.name][1]+(Hyper_CHH[Hyper_CHH[0]==sector.name][2]-Hyper_CHH[Hyper_CHH[0]==sector.name][1])/2)y = np.array(Hyper_CHH[Hyper_CHH[0]==sector.name][3])track1.scatter(x, y,c="#23BAC5",vmin=0,vmax=1,ec="black",lw=0.1,alpha=0.5)# Add cytoband tracks from cytoband filegene_density=pd.read_table("gene.density",header=None)track2 = sector.add_track((70, 75), r_pad_ratio=0.1)# 注释圈名if sector.name == circos.sectors[0].name:circos.text("Gene", r=track2.r_center, size=8)track2.heatmap(gene_density[gene_density[0]==sector.name][3],vmin=0,vmax=max(gene_density[3]),cmap="Greys")# 内圈的负值Hyper_CG=pd.read_table("CG_gDMR.hypo.txt",header=None)x = np.array(Hyper_CG[Hyper_CG[0]==sector.name][1]+(Hyper_CG[Hyper_CG[0]==sector.name][2]-Hyper_CG[Hyper_CG[0]==sector.name][1])/2)y = np.array(Hyper_CG[Hyper_CG[0]==sector.name][3])track3 = sector.add_track((50, 70), r_pad_ratio=0.1)# 注释圈名if sector.name == circos.sectors[0].name:circos.text("Hypo", r=track3.r_center, size=8)for r in range(52, 70, 2):sector.line(r=r, ec="lightgrey")track3.scatter(x, y,c="#EECA40",vmin=-1,vmax=0,ec="black",lw=0.1,alpha=0.5)Hyper_CHG=pd.read_table("CHG_gDMR.hypo.txt",header=None)x = np.array(Hyper_CHG[Hyper_CHG[0]==sector.name][1]+(Hyper_CHG[Hyper_CHG[0]==sector.name][2]-Hyper_CHG[Hyper_CHG[0]==sector.name][1])/2)y = np.array(Hyper_CHG[Hyper_CHG[0]==sector.name][3])track3.scatter(x, y,c="#FD763F",vmin=-1,vmax=0,ec="black",lw=0.1,alpha=0.5)Hyper_CHH=pd.read_table("CHH_gDMR.hypo.txt",header=None)x = np.array(Hyper_CHH[Hyper_CHH[0]==sector.name][1]+(Hyper_CHH[Hyper_CHH[0]==sector.name][2]-Hyper_CHH[Hyper_CHH[0]==sector.name][1])/2)y = np.array(Hyper_CHH[Hyper_CHH[0]==sector.name][3])track3.scatter(x, y,c="#23BAC5",vmin=-1,vmax=0,ec="black",lw=0.1,alpha=0.5)circos.colorbar(bounds=(0.35, 0.55, 0.3, 0.01),vmin=0,vmax=max(gene_density[3]),orientation="horizontal",cmap="Greys")
fig = circos.plotfig()# 图注画在圈中间
_ = circos.ax.legend(handles=[Line2D([], [], color="#EECA40", marker="o", label="CG", ms=6, ls="None"),Line2D([], [], color="#FD763F", marker="o", label="CHG", ms=6, ls="None"),Line2D([], [], color="#23BAC5", marker="o", label="CHH", ms=6, ls="None"),],bbox_to_anchor=(0.5, 0.45),loc="center",ncols=1,
)
fig.savefig("Circos.pdf") # 保存
fig.savefig("Circos.png",dpi=300)

更多内容敬请关注微信公众号:
在这里插入图片描述


http://www.mrgr.cn/news/36851.html

相关文章:

  • liunxcentos7下 跟目录空间不足docker load镜像报错空间不足
  • 记录PC端处理图片颜色翻转的java、JNA、JNI处理耗时情况
  • 深度拆解:如何在Facebook上做跨境电商?
  • TinyWebSever项目面试题整理
  • 曦烽盛启——丰宁坝上草原大拍暨曦烽摄影学院十五周年庆典活动
  • 2025年第十三届生物信息学与计算生物学国际会议(ICBCB 2025)即将召开!
  • 文件防泄密软件有哪些|5款简单易操作软件有效防止文件泄密!
  • 美客多在巴西疯狂扩张,支持时尚和超市品类增长,美客多怎么上传产品?
  • vue 中获取数值但是只获取到了 Promise 属性,获取不到其中的值
  • 传递正能量
  • AWS Network Firewall - 配置只应许白名单域名出入站
  • MySQL数据库进阶知识(五)《锁》
  • JVM 基本组成
  • 【工具分享】Onyx2勒索病毒解密工具
  • GO Message Bus
  • 希捷电脑硬盘好恢复数据吗?探讨可能性、方法以及注意事项
  • 生物技术新突破:DAP-seq技术在非模式物种中的应用
  • linux-自动化构建工具与进度条程序的实现
  • 常用CSS选择器
  • Altium Designer创建接口类器件原理图库