2025.06.20【pacbio】|PB甲基化分析结果的统计与可视化介绍

VIP2025.06.20【pacbio】|PB甲基化分析结果的统计与可视化介绍
本文介绍了PacBio甲基化数据的统计分析与可视化方法。主要内容包括:1) 分析甲基化位点分布、比例及染色体分布;2) 修饰基序(Motif)的统计与热图展示;3) 全基因组和局部区域测序深度可视化;4) 使用IGV查看单个reads修饰情况。文章提供了R代码示例,涵盖柱状图、分布图、热图等多种可视化形式,并建议了差异甲基化分析等进阶方法,帮助科研人员深入挖掘甲基化数据。
开通可读全文 ·2k 人阅读
最低0.47元/天 解锁文章

引言

之前介绍了pacbio的甲基化基础分析流程,在完成了PacBio甲基化数据的自动化分析流程后,如何对结果进行统计和可视化,是科研工作中不可或缺的一环。高质量的图表不仅能帮助我们直观理解数据,还能为论文和报告增色不少。本文将以Snakemake流程输出的典型结果为例,详细介绍甲基化分析中的常见统计方法和可视化思路,并给出R/Python代码示例。


请添加图片描述

1. 甲基化分析结果文件简介

在前文的Snakemake流程中,主要输出了以下几类与统计和可视化相关的文件:

  • *.significant.bed.xls:高可信度甲基化位点表
  • motifs.csv:鉴定出的甲基化修饰基序
  • *.basemods.gff:GFF格式的修饰注释
  • *.kinetics.align.bam:比对BAM文件(可用于IGV可视化)
  • *.pdf/*.png:全基因组覆盖度分布图

2. 甲基化位点统计分析

2.1 统计不同类型修饰的数量和分布

目标:统计全基因组范围内各类甲基化修饰(如m6A、m4C)的数量、分布及其在不同染色体上的分布。

示例R代码:统计m6A/m4C位点数

*.significant.bed.xls 文件格式示例
该文件为制表符分隔文本,首行为表头,常见字段如下:

#seq_name	pos_left	pos	depth	strand	type	fraction
chr1	    12345	    12346	25	    +	    m6A	    0.85
chr1	    23456	    23457	18	    -	    m4C	    0.67
chr2	    34567	    34568	30	    +	    m6A	    0.92
...

seq_name:染色体或contig名称
pos_left:起始位置(0-based)
pos:终止位置(1-based)
depth:该位点的测序深度
strand:链信息(+/-)
type:修饰类型(如m6A、m4C)
fraction:甲基化比例(0~1)

library(data.table)
bed <- fread("sample1.basemods.significant.bed.xls", skip=1)
table(bed$type)
# 输出每种修饰类型的数量
可视化:不同修饰类型的柱状图
library(ggplot2)
ggplot(bed, aes(x=type, fill=type)) 
VIP 2025.06.20【pacbio】|使用Snakemake构建可重复的PacBio全基因组甲基化分析流程 本文介绍了一个基于Snakemake的自动化流程,用于分析PacBio SMRT测序数据的DNA甲基化。流程涵盖从原始数据比对到甲基化检测的全过程,包含参考基因组索引构建、序列比对、修饰位点识别等核心步骤。通过配置文件和样本信息文件实现参数自定义,确保结果可重复性。该流程整合了pbmm2、ipdSummary等PacBio官方工具,并通过自动化处理减少手动操作错误,适合批量分析多个样本。输出文件包括GFF格式修饰位点、CSV统计表以及可视化文件,为表观遗传学研究提供完整分析方案。 开通可读全文·1.6k 阅读·23 点赞·16 收藏 阅读全文

相关推荐

IGV——基因组可视化

例如,查看转录组数据不同区域的Reads丰度,不同样本的基因表达差异,查找SNP位点,查看甲基化,查看基因结构信息。(3)如果在本地可视化,需要将 bam 和 bai 或 bigwig 文件从 linux 系统上传到本地,且放置在同一个文件夹下,如果服务器搭建有 ftp 站点则可以直接使用 url 载入,也可以将数据上传至可以接受外部访问的服务器(如 https://de.cyverse.org/dashboard )上再使用 url 载入,具体可参考教程《通过 url 访问服务器文件》。

import__的博客 6647

甲基化分析流程

目的:找到基因组哪些区域甲基化,以及附近基因有哪些?

weixin_46795596的博客 1805

生物信息学——三代测序数据:Pacbio

受限于技术力问题,面对基因组长达数百万乃至数十亿个碱基,目前很难对完整的基因组进行测序。于是研究人员尝试利用小片段重构基因组序列信息。根据不同的测序仪,将基因组碎片化为合适的大小,然后才能进行实际的测序、重组和分析。基于它们的底层化学和所分析的DNA片段的长度,测序仪可以分为长读长和短读长两种类型。长读长测序检测长度为1,000到20,000个碱基或更长的DNA(或RNA)片段。无需经过扩增,可以直接对提取后的基因测序。短读长测序多数只能检测50-300个碱基长度的片段。

qq_45794091的博客 6387

ArcGIS 10.4 Desktop安装步骤。

步骤 1 :安装 LicenseManager 有了安装包以后,我们进入安装包文件目录,如图所示,首先需要安装ArcGIS license manager,即点击打开图中红色框中exe文件。 PS:安装 LicenseManager,关于 License 的安装路径一般默认即可,因为很小的,装C 盘也无所谓了,当然了自己改了也可以的,但是不建议! 即便修改路径,建议只改盘符,不修改相对路径方式! 如:C:\Program Files (x86)\ArcGIS\License10.6 修改为 D:\.

IT爱好者 2万+

【中软软件测试笔试题】中国软件与技术服务股份有限公司笔试一卷

请分别写出ISO7层架构,并分别说明TCP、UDP协议分别在处于哪一层,以及TCP、 UPD协议有什么区别。 答:7 应用层 6 表示层 5 会话层 4 传输层 3网络层 2 数据链路层 1 物理层 第四层;UDP是面向无连接、不可靠数据传输协议; TCP是面向连接、可靠数据传输协议......

Yanan990830的博客 1808

ST7102驱动移植与调试实战

本文详细介绍了在君正T40平台上移植和调试ST7102 LCD驱动芯片的实战经验。内容涵盖硬件接口配置、时序参数调整、寄存器初始化序列编写及常见问题解决,重点解析了ST7102单芯片集成显示与触控功能的特性,为嵌入式设备显示驱动开发提供实用指南。

k6l7m8n9的博客 1390

深度学习注意力机制新突破:YOLOv8结合SimAM无参数注意力模块,让目标检测更精准!

在计算机视觉领域,注意力机制已成为提升模型性能的关键技术。然而,传统的注意力模块如SENet、CBAM等往往引入大量额外参数,增加了模型的复杂度和计算负担。本文介绍一种创新的无参数注意力机制——SimAM(Simple Parameter-Free Attention Module),并将其成功应用于YOLOv8目标检测框架。SimAM基于神经科学理论,通过定义能量函数为特征图计算3D注意力权重,无需任何额外参数即可显著提升模型性能。本文详细阐述了SimAM的理论基础、实现原理,并提供了完整的YOLOv8集

YOLO 93

5.7 BGP属-ORIGIN

注意:当去往同一个目的地存在多条不同Origin属性的路由时,在其他条件都相同的情况下,BGP将按如Origin的下顺序优选路由:IGP > EGP > Incomplete。通过以上输出可以看到路由表选择了R3发过来的路由,因为R1的路由的起源为?通过以上输出我们可以看到针对于1.1.1.1这条路由有两条,从R1和R3传过来的起源都属于IGP。(4)把从R1传过来的路由的起源改为incomplete。通过以上输出可以看到,从R1传过来的路由是最优的。(2)在R2上查看1.1.1.1的路由详细信息。

lwljh134的博客 301

《从零学会 CST 仿真全流程》专栏简介

本专栏系统讲解CSTPCBStudio与MWSStudio在高速PCB仿真中的应用,涵盖SI/PI仿真全流程。重点介绍CSTPCBStudio用于板级信号完整性分析,CSTMWS用于3D高频电磁仿真,以及电缆、低频、多物理场等配套工具。课程以实战为导向,通过USB3.1、LPDDR4等案例,手把手教学软件操作与工程应用,适合硬件工程师、EMC工程师及电子专业学生学习。内容注重实操,弱化理论,帮助快速掌握仿真核心技能。

一只豌豆象的博客 590

【linux】解决报错:Network error: Connection refused

【linux】解决报错:Network error: Connection refused

weixin_44502754的博客 4929

温度PID控制PLC程序

温度PID稳定控制程序,适用于欧姆龙PLC,对其他PLC也具有很大参考价值

深度学习工具包 DeepLearningToolbox

深度学习工具包Deprecation notice.-----This toolbox is outdated and no longer maintained.There are much better tools available for deep learning than this toolbox, e.g. [Theano](http://deeplearning.net/software/theano/), [torch](http://torch.ch/) or [tensorflow](http://www.tensorflow.org/)I would suggest you use one of the tools mentioned above rather than use this toolbox.Best, Rasmus.DeepLearnToolbox================A Matlab toolbox for Deep Learning.Deep Learning is a new subfield of machine learning that focuses on learning deep hierarchical models of data.It is inspired by the human brain's apparent deep (layered, hierarchical) architecture.A good overview of the theory of Deep Learning theory is[Learning Deep Architectures for AI](http://www.iro.umontreal.ca/~bengioy/papers/ftml_book.pdf)For a more informal introduction, see the following videos by Geoffrey Hinton and Andrew Ng.* [The Next Generation of Neural Networks](http://www.youtube.com/watch?v=AyzOUbkUf3M) (Hinton, 2007)* [Recent Developments in Deep Learning](http://www.youtube.com/watch?v=VdIURAu1-aU) (Hinton, 2010)* [Unsupervised Feature Learning and Deep Learning](http://www.youtube.com/watch?v=ZmNOAtZIgIk) (Ng, 2011)If you use this toolbox in your research please cite [Prediction as a candidate for learning deep hierarchical models of data](http://www2.imm.dtu.dk/pubdb/views/publication_details.php?id=6284)```@MASTERSTHESIS\{IMM2012-06284, author = "R. B. Palm", title = "Prediction as a candidate for learning deep hierarchical models of data", year = "2012",}```Contact: rasmusbergpalm at gmail dot comDirectories included in the toolbox-----------------------------------`NN/` - A library for Feedforward Backpropagation Neural Networks`CNN/` - A library for Convolutional Neural Networks`DBN/` - A library for Deep Belief Networks`SAE/` - A library for Stacked Auto-Encoders`CAE/` - A library for Convolutional Auto-Encoders`util/` - Utility functions used by the libraries`data/` - Data used by the examples`tests/` - unit tests to verify toolbox is workingFor references on each library check REFS.mdSetup-----1. Download.2. addpath(genpath('DeepLearnToolbox'));Example: Deep Belief Network---------------------```matlabfunction test_example_DBNload mnist_uint8;train_x = double(train_x) / 255;test_x = double(test_x) / 255;train_y = double(train_y);test_y = double(test_y);%% ex1 train a 100 hidden unit RBM and visualize its weightsrand('state',0)dbn.sizes = [100];opts.numepochs = 1;opts.batchsize = 100;opts.momentum = 0;opts.alpha = 1;dbn = dbnsetup(dbn, train_x, opts);dbn = dbntrain(dbn, train_x, opts);figure; visualize(dbn.rbm{1}.W'); % Visualize the RBM weights%% ex2 train a 100-100 hidden unit DBN and use its weights to initialize a NNrand('state',0)%train dbndbn.sizes = [100 100];opts.numepochs = 1;opts.batchsize = 100;opts.momentum = 0;opts.alpha = 1;dbn = dbnsetup(dbn, train_x, opts);dbn = dbntrain(dbn, train_x, opts);%unfold dbn to nnnn = dbnunfoldtonn(dbn, 10);nn.activation_function = 'sigm';%train nnopts.numepochs = 1;opts.batchsize = 100;nn = nntrain(nn, train_x, train_y, opts);[er, bad] = nntest(nn, test_x, test_y);assert(er < 0.10, 'Too big error');```Example: Stacked Auto-Encoders---------------------```matlabfunction test_example_SAEload mnist_uint8;train_x = double(train_x)/255;test_x = double(test_x)/255;train_y = double(train_y);test_y = double(test_y);%% ex1 train a 100 hidden unit SDAE and use it to initialize a FFNN% Setup and train a stacked denoising autoencoder (SDAE)rand('state',0)sae = saesetup([784 100]);sae.ae{1}.activation_function = 'sigm';sae.ae{1}.learningRate = 1;sae.ae{1}.inputZeroMaskedFraction = 0.5;opts.numepochs = 1;opts.batchsize = 100;sae = saetrain(sae, train_x, opts);visualize(sae.ae{1}.W{1}(:,2:end)')% Use the SDAE to initialize a FFNNnn = nnsetup([784 100 10]);nn.activation_function = 'sigm';nn.learningRate = 1;nn.W{1} = sae.ae{1}.W{1};% Train the FFNNopts.numepochs = 1;opts.batchsize = 100;nn = nntrain(nn, train_x, train_y, opts);[er, bad] = nntest(nn, test_x, test_y);assert(er < 0.16, 'Too big error');```Example: Convolutional Neural Nets---------------------```matlabfunction test_example_CNNload mnist_uint8;train_x = double(reshape(train_x',28,28,60000))/255;test_x = double(reshape(test_x',28,28,10000))/255;train_y = double(train_y');test_y = double(test_y');%% ex1 Train a 6c-2s-12c-2s Convolutional neural network %will run 1 epoch in about 200 second and get around 11% error. %With 100 epochs you'll get around 1.2% errorrand('state',0)cnn.layers = { struct('type', 'i') %input layer struct('type', 'c', 'outputmaps', 6, 'kernelsize', 5) %convolution layer struct('type', 's', 'scale', 2) %sub sampling layer struct('type', 'c', 'outputmaps', 12, 'kernelsize', 5) %convolution layer struct('type', 's', 'scale', 2) %subsampling layer};cnn = cnnsetup(cnn, train_x, train_y);opts.alpha = 1;opts.batchsize = 50;opts.numepochs = 1;cnn = cnntrain(cnn, train_x, train_y, opts);[er, bad] = cnntest(cnn, test_x, test_y);%plot mean squared errorfigure; plot(cnn.rL);assert(er<0.12, 'Too big error');```Example: Neural Networks---------------------```matlabfunction test_example_NNload mnist_uint8;train_x = double(train_x) / 255;test_x = double(test_x) / 255;train_y = double(train_y);test_y = double(test_y);% normalize[train_x, mu, sigma] = zscore(train_x);test_x = normalize(test_x, mu, sigma);%% ex1 vanilla neural netrand('state',0)nn = nnsetup([784 100 10]);opts.numepochs = 1; % Number of full sweeps through dataopts.batchsize = 100; % Take a mean gradient step over this many samples[nn, L] = nntrain(nn, train_x, train_y, opts);[er, bad] = nntest(nn, test_x, test_y);assert(er < 0.08, 'Too big error');%% ex2 neural net with L2 weight decayrand('state',0)nn = nnsetup([784 100 10]);nn.weightPenaltyL2 = 1e-4; % L2 weight decayopts.numepochs = 1; % Number of full sweeps through dataopts.batchsize = 100; % Take a mean gradient step over this many samplesnn = nntrain(nn, train_x, train_y, opts);[er, bad] = nntest(nn, test_x, test_y);assert(er < 0.1, 'Too big error');%% ex3 neural net with dropoutrand('state',0)nn = nnsetup([784 100 10]);nn.dropoutFraction = 0.5; % Dropout fraction opts.numepochs = 1; % Number of full sweeps through dataopts.batchsize = 100; % Take a mean gradient step over this many samplesnn = nntrain(nn, train_x, train_y, opts);[er, bad] = nntest(nn, test_x, test_y);assert(er < 0.1, 'Too big error');%% ex4 neural net with sigmoid activation functionrand('state',0)nn = nnsetup([784 100 10]);nn.activation_function = 'sigm'; % Sigmoid activation functionnn.learningRate = 1; % Sigm require a lower learning rateopts.numepochs = 1; % Number of full sweeps through dataopts.batchsize = 100; % Take a mean gradient step over this many samplesnn = nntrain(nn, train_x, train_y, opts);[er, bad] = nntest(nn, test_x, test_y);assert(er < 0.1, 'Too big error');%% ex5 plotting functionalityrand('state',0)nn = nnsetup([784 20 10]);opts.numepochs = 5; % Number of full sweeps through datann.output = 'softmax'; % use softmax outputopts.batchsize = 1000; % Take a mean gradient step over this many samplesopts.plot = 1; % enable plottingnn = nntrain(nn, train_x, train_y, opts);[er, bad] = nntest(nn, test_x, test_y);assert(er < 0.1, 'Too big error');%% ex6 neural net with sigmoid activation and plotting of validation and training error% split training data into training and validation datavx = train_x(1:10000,:);tx = train_x(10001:end,:);vy = train_y(1:10000,:);ty = train_y(10001:end,:);rand('state',0)nn = nnsetup([784 20 10]); nn.output = 'softmax'; % use softmax outputopts.numepochs = 5; % Number of full sweeps through dataopts.batchsize = 1000; % Take a mean gradient step over this many samplesopts.plot = 1; % enable plottingnn = nntrain(nn, tx, ty, opts, vx, vy); % nntrain takes validation set as last two arguments (optionally)[er, bad] = nntest(nn, test_x, test_y);assert(er < 0.1, 'Too big error');```[![Bitdeli Badge](https://d2weczhvl823v0.cloudfront.net/rasmusbergpalm/deeplearntoolbox/trend.png)](https://bitdeli.com/free "Bitdeli Badge")

jre9下载 jdk9

(Java SE Runtime Environment 9.0) 9.0u175 官方最新版 【x86|x64】

SAP B1中文培訓教材_SDK

公司内部培训资料,对开发工具 ,DIAPI, UIAPI,Add-ON有详细解说

上一篇: 2025.06.20【pacbio】|使用Snakemake构建可重复的PacBio全基因组甲基化分析流程
下一篇: 2025.06.23【甲基化】methylKit:甲基化测序数据分析安装与详细使用教程

双11狂欢大促

盲盒赢 定制羽绒服

开启惊喜盲盒,抽取你的专属好礼

— 惊喜奖品抢先看 —
穆易青
博客等级 码龄11年 5706粉丝 397原创
评论
成就一亿技术人!
拼手气红包6.0元
还能输入1000个字符
 
 条评论被折叠 查看
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付:10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

穆易青

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值