modPhred

最新推荐文章于 2025-06-05 19:57:55 发布

原创

最新推荐文章于 2025-06-05 19:57:55 发布

· 663 阅读

21 ·

版权

文章标签：

#服务器 #网络 #运维 #青少年编程 #考研

比对

比对使用minimap2进行，排序后的BAM文件存储在out_dir/minimap2目录中。

我们使用推荐的设置来映射ONT读取：-ax map-ont。对于RNA样本，自动执行剪接感知映射：-ax splice。

由于修饰状态在FastQ中编码，它将通过下游分析（如比对）传播，并存储在BAM文件中。在BAM中编码修饰状态允许直接在基因组浏览器中可视化修饰概率。

你可以在这里 <output.html#data-formats>_找到更多关于数据格式的信息。

碱基识别

碱基识别使用guppy碱基识别器进行。你可以选择：

• 在之前进行Fast5的碱基识别（支持guppy版本>3.2）
• 或运行实时碱基识别（支持guppy版本>3.6）

在实时模式下，碱基识别使用默认设置进行。默认情况下，使用CpG/dam/dcm模型（dna_r9.4.1_450bps_modbases_dam-dcm-cpg_hac.cfg）。这可以使用-c / --config参数更改。

注意，如果你选择运行实时碱基识别，你需要安装匹配版本的pyguppyclient <install.html#which-pyguppyclient-version-should-i-install>_。

修饰感知模型

碱基识别应使用修饰感知模型。目前，没有公开可用的RNA修饰模型。本仓库中提供的RNA模型是在体外转录的分子上训练的，其中所有碱基都被修饰或全部未修饰，因此它不适用于仅部分碱基预期被修饰的生物样本。你可以在手稿的补充信息中找到更多细节。

读取聚类

你也可以根据修饰状态对读取进行聚类。注意，这目前聚类整个参考序列，因此仅在单个转录本上运行是合理的。绝对避免聚类来自整个染色体的读取。

~/src/modPhred/src/mod_cluster.py --minfreq 0.01 -i mod.gz -e pdf

这将生成类似于以下的图：

你可以使用-r指定感兴趣的区域，例如，如果你想聚类在chr1上100到200位置的正链上对齐的读取，你可以这样运行：

~/src/modPhred/src/mod_cluster.py --minfreq 0.01 -i mod.gz -e pdf -r chr1:100-200+

聚类细节

我们按以下步骤进行读取聚类：

• 用户选择读取聚类的区域，即某些转录本或链式基因组区域
• 识别该区域中被修饰的一组位置，即超过20%的读取被修饰（--minfreq 0.2），并且在至少一个样本中至少有25个读取对齐（-d / --mindepth 25）。重要的是，我们默认包括所有类型的修饰。可以使用--mod将分析限制为单一修饰。
• 从所有样本的所有读取中提取一组修饰位置的修饰概率。注意，读取通过映射质量-m / --mapq 15和它们与选定区域的重叠进行过滤：默认情况下，读取必须覆盖至少80%的选定区域--minAlgFrac 0.8。
• 执行层次聚类并保存图。仅聚类读取（修饰位置不聚类）。我们使用seaborn.clustermap <https://seaborn.pydata.org/generated/seaborn.clustermap.html>_的实现。

目前我们不报告读取聚类。

位置之间的相关性

你可以绘制某些区域内修饰位置之间的相关性。这很简单：

~/src/modPhred/src/mod_correlation.py -i modPhred/PRJEB22772/mod.gz -r NC_000913.3:1-5000

这将生成类似以下的图：

此外，使用以下命令可以将这些图缩小到仅：
• 来自特定链的修饰
• 和来自特定链的一个修饰

~/src/modPhred/src/mod_correlation.py -i modPhred/PRJEB22772/mod.gz -r NC_000913.3:1-5000+
~/src/modPhred/src/mod_correlation.py -i modPhred/PRJEB22772/mod.gz -r NC_000913.3:1-5000+ --mod 6mA

修饰的编码

由于碱基识别质量对于纳米孔测序来说不是非常有用（通常在7-12之间），我们决定将修饰概率存储为FastQ碱基质量。这是通过以下方式实现的。

Guppy碱基识别器在Fast5文件中报告碱基被修饰的概率，位于/Analyses/Basecall_1D_000/BaseCalled_template/ModBaseProbs下。这些概率存储为从0（无修饰）到255（修饰）的8位整数，并分别报告给定模型训练的所有修饰。如果我们只想存储每个碱基的一个修饰，我们可以简单地将这些值重新缩放到ASCII范围（0-93）并存储（FastQ文件中的碱基质量存储为ASCII字符）。然而，特别是对于RNA修饰，存储给定碱基的多个可能修饰的信息将是有益的（例如，m5C和5hmC对于C来说非常常见），因为RNA中有超过170种已知的修饰。如果你想存储每个碱基的3个修饰信息，这总共最多12个修饰，可以将修饰概率重新缩放到31个值（而不是255或93），并存储给定碱基具有最高概率的修饰的概率。例如，如果你想存储关于C的3个修饰（m5C，5hmC和m3C）的信息，你可以分配以下值：

0-30到第一个（m5C，PHRED 33-63），
31-61到第二个（5hmC，PHRED 64-94）
和62-92到第三个（m3C，PHRED 95-125）修饰的C。

现在想象一下，在某些读取中，碱基C（下面的位置2635）对于m5C，5hmC和m3C的概率分别为12，247，9。由于每个碱基只能存储一个概率，最有信息量的将是存储5hmC的概率，因为它具有最高的概率。我们将247的概率重新缩放到31单位范围，如下所示：

Q = 255 * max{p1, p2, p3} / 31 + 31 * Mi

其中Mi（修饰索引）是

• 0如果第一个（m5C），
• 1如果第二个（5hmC）
• 或2如果第三个（m3C）修饰的C
• 具有最高的修饰概率（p）。

这种信息存储方式保证了简单性和多功能性。由于修饰概率存储在FastQ中，不需要外部数据库或额外文件来计算修饰，因为所有读取的每个碱基修饰概率也将存储在从FastQ派生的BAM文件中。此外，我们的编码系统是灵活的 - 它可以轻松调整以编码更多修饰。例如，每个碱基9个修饰，总共最多36个修饰，可以在每个修饰的10值范围内编码（而不是默认的3个修饰的31值）。

理论上，我们的系统允许在FastQ格式中直接编码多达368个修饰（每个碱基92个）的信息，前提是我们将修饰状态的信息限制为二进制形式，即碱基携带368个修饰中的一个或不携带。

modPhred: 文档

什么是modPhred？

modPhred是一个用于检测、注释和可视化DNA/RNA修饰的管道。管道由四个步骤/模块组成：

modEncode：在FastQ中编码修饰概率（mod_encode.py）
modAlign：构建比对，保留BAM中的修饰信息（mod_align.py）
modReport：提取RNA修饰信息（bedGraph）和QC报告（mod_report.py）
modAnalysis：

a. 绘制QC统计、配对图和维恩图（mod_plot.py），
b. 修饰的共现（mod_correlation.py）
c. 基于修饰配置的每读取聚类（mod_cluster.py）

所有这些脚本可以单独运行或通过执行modPhred/run作为管道运行。你可以在方法部分找到更多细节。