mAFiA

mAFiA - (Another) m6A Finding Algorithm

这里我们提供了一个运行mAFiA的简要指南，以X染色体为例。

0. 预备步骤

以下步骤已在Python版本3.9上测试。获取代码并激活虚拟环境，例如：

git clone https://github.com/dieterich-lab/mAFiA.git
cd mAFiA
python3 -m venv mafia-venv
source mafia-venv/bin/activate

如果您的pip版本小于21.3，则将其升级到较新版本：

python3 -m pip install --upgrade pip

安装包：

pip install -e .

从此处下载模型和数据：

"models"文件夹包含：

• backbone.torch: 基于RODAN的神经网络，用于基础调用和特征提取
• backbone.config: backbone的训练配置
• classifiers: 打包的逻辑回归模型

"data"文件夹包含X染色体的输入数据子集：

• fast5_chrX: HEK293 WT mRNA的dRNA-Seq原始数据
• GRCh38_96.X: 基因组参考
• GLORI_chrX.bed: bed格式的查询修饰位点。此文件特别对应于GLORI中列出的位点。

假设数据和模型分别解压到${data}和${model}。您的输出目录是${output}：

backbone="${models}/backbone.torch"
classifiers="${models}/classifiers"
fast5dir="${data}/fast5_chrX"
ref="${data}/GRCh38_96.X.fa"
mod="${data}/GLORI_chrX.bed"

mkdir -p "${output}"
basecall="${output}/rodan.fasta"
bam="${output}/minimap.q50.bam"

1. 基础调用

基础调用脚本改编自RODAN仓库。假设${mafia}是您的代码目录。

python3 ${mafia}/RODAN/basecall.py \
--fast5dir ${fast5dir} \
--model ${backbone} \
--batchsize 4096 \
--outdir ${output}

在一台相当现代的GPU机器上，这一步应该在30分钟内完成。

2. 对齐

将基础调用结果与参考基因组对齐。过滤、排序并索引BAM文件。

minimap2 --secondary=no -ax splice -uf -k14 -t 36 --cs ${ref} ${basecall} \
| samtools view -bST ${ref} -q50 - \
| samtools sort - > ${bam}

samtools index ${bam}

3. mAFiA

经过标准程序后，我们现在可以测量${mod}中指定位点的m6A化学计量学。

test_mAFiA \
--bam_file ${bam} \
--fast5_dir ${fast5dir} \
--ref_file ${ref} \
--mod_file ${mod} \
--min_coverage 50 \
--max_num_reads 1000 \
--backbone_model_path ${backbone} \
--classifier_model_dir ${classifiers} \
--mod_prob_thresh 0.5 \
--out_dir ${output}

最后一步应该在1小时内完成。我们目前正在努力将特征提取步骤直接集成到基础调用中。运行时间应该会显著减少。

在${output}目录中，您现在应该看到两个文件：

• “mAFiA.sites.bed”: 覆盖度高于最低阈值（默认50）的位点列表。"modRatio"列列出了位点的化学计量学。
• “mAFiA.reads.bam”: 与输入BAM文件${bam}中相同的对齐读取，但增加了MM和ML标签，分别标记了每个单独读取中的位置和修饰概率。结果可以使用IGV等工具进行可视化。

所有染色体的完整HEK293 WT数据集可以从zenodo下载。