当前位置: 首页 > 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.lryc.cn/news/446792.html

相关文章:

  • Redis Sorted Set 跳表的实现原理和分析
  • 新手教学系列——在MySQL分表中批量调整表结构的实践与优化
  • 解决事务提交延迟问题:Spring中的事务绑定事件监听机制解析
  • Python 异步编程的秘密武器:Asyncio
  • 10年计算机考研408-计算机网络
  • 深信服校招面试总结
  • 【LeetCode热题100】模拟
  • 如何在Chrome最新浏览器中调用ActiveX控件?
  • 一款好用的远程连接工具:MobaXterm
  • Spring Boot使用配置方式整合MyBatis
  • HarmonyOS第一课-应用程序框架基础习题答案
  • 滚雪球学SpringCloud[10.2讲]:微服务项目的性能优化与调优
  • EasyExcel将数据库里面的数据生成excel文件
  • 【YOLO学习】YOLOv1详解
  • HarmonyOS应用开发(组件库)--组件模块化开发、工具包、设计模式(持续更新)
  • python测试开发---前后端交互Axios
  • 删除视频最后几帧 剪切视频
  • SSM框架学习(四、SpringMVC实战:构建高效表述层框架)
  • 戴尔笔记本电脑——重装系统
  • 领夹麦克风哪个品牌音质最好,主播一般用什么麦克风
  • 华为静态路由(route-static)
  • Focalboard开源项目管理系统本地Windows部署与远程访问协同办公
  • Java如何操作Elasticsearch
  • cpu路、核、线程、主频、缓存
  • 【AI算法岗面试八股面经【超全整理】——深度学习】
  • STL——map和set【map和set的介绍和使用】【multimap和multiset】
  • 【笔记】神领物流配置本地hosts无法访问域名(排除DNS 排除文件编码问题)已解决
  • Java | Leetcode Java题解之第424题替换后的最长重复字符
  • Xcode 16 Pod init 报错
  • 【数据结构】Java的HashMap 和 HashSet 大全笔记,写算法用到的时候翻一下,百度都省了!(实践篇)