medaka

一、Getting Started

1. Sequence correction

安装软件后，medaka可以使用其默认设置运行medaka_consensus。该程序使用samtools和minimap2，如果medaka是使用from-source方法安装的，那么这些将出现在medaka环境中，否则需要用户提供。

source ${MEDAKA}  # i.e. medaka/venv/bin/activate
NPROC=$(nproc)
BASECALLS=basecalls.fa
DRAFT=draft_assm/assm_final.fa
OUTDIR=medaka_consensus
medaka_consensus -i ${BASECALLS} -d ${DRAFT} -o ${OUTDIR} -t ${NPROC} -m r941_min_fast_g330

当medaka_consensus运行完成后，consensus将保存到${OUTDIR}/consensus.fasta中.

Warning
为了获得最佳结果，根据使用的basecaller，指定正确的模型-m 是很重要的。允许的值可以通过运行medaka tools list_models找到。

Medaka模型的命名
i) the pore type,
ii) the sequencing device (MinION or PromethION),
iii) the basecaller variant,
iv) the basecaller version:

{pore}_{device}_{caller variant}_{caller version}

例如，模型r941_min_fast_g303
i) R9.4.1 flowcells
ii) 测序设备MinION
iii) 使用fast Guppy basecaller
iv) basecaller 版本3.0.3

模型r941_prom_hac_g303
i) R9.4.1 flowcells
ii) 测序设备PromethION
iii) 高精度basecaller （Guppy配置文件中称为“hac”）
iv) basecaller 版本3.0.3

如果使用的Guppy版本没有完全对应的medaka模型，则应选择最高版本等于或小于Guppy版本的medaka模型。

2. 提高并行性

medaka_consensus 程序适用于简单的数据集，但对于运行大型数据集可能不是最佳选择。可以通过独立运行medaka_consensus的组件步骤，实现更高级别的并行性。该程序执行三项tasks：

  1. reads比对到基因组（通过mini_align，它是minimap2上的一个包装）
  2. 在assembly 区域运行consensus算法（medaka consensus，注意不要下划线！）
  3. 对第二步产生的结果聚合。创建consensus sequences (medaka stitch)

这三个步骤是分离的，可以独立运行。在大多数情况下，步骤2是瓶颈，可以简单地并行化。medaka consensus可以提供一个--regions参数，该参数将限制其在步骤1中输出的.bam文件中的特定组装基因组序列。因此，可以同时为多个组装基因组序列运行单个作业。在最后一步中，medaka stitch可以将步骤2输出的一个或多个.hdf文件作为输入。

# align reads to assembly
mini_align  \
    -i basecalls.fasta  \
    -r assembly.fasta  \
    -P -m \
    -p calls_to_draft.bam  \
    -t <threads>


# run lots of jobs like this, change model as appropriate
mkdir results
medaka consensus calls_to_draft.bam results/contigs1-4.hdf \
    --model r941_min_fast_g303 --batch 200 --threads 8 \
    --region contig1 contig2 contig3 contig4
...


# wait for jobs, then collate results
medaka stitch results/*.hdf polished.assembly.fasta

不建议为medaka consensus指定一个大于8的--threads值，因为超过这个值计算效率很差。还请注意，medaka consensus可能会使用相当于<threads> + 4的资源，因为另外4个线程用于读取和准备输入数据。

二、Benchmarks

三、Walkthrough

下面展示了如何根据从R9.4 flowcells和ligation sequencing试剂盒获得的测序数据训练和使用一种简单形式的 consensus 模型。

下面是一个简单的演示过程。它并不代表 best-practices 或state-of-the-art workflow。为了能很好的训练出适用性更广的模型，需要准备更大的数据集。Medaka的标准模型是用katuali训练的.

获取数据和软件

下载ONT的basecalls的训练数据：

WALKTHROUGH=${PWD}/medaka_walkthrough
mkdir -p ${WALKTHROUGH} && cd ${WALKTHROUGH}
wget https://s3-eu-west-1.amazonaws.com/ont-research/medaka_walkthrough_no_reads.tar.gz
tar -xvf medaka_walkthrough_no_reads.tar.gz
DATA=${PWD}/data

提取的文件还包含在以下过程中创建的所有中间输出文件。只需将必需的子文件夹从${DATA}目录复制到${WALKTHROUGH}目录，就可以跳过任何步骤。

The necessary software can be sourced using:

cd ${WALKTHROUGH}
git clone https://github.com/nanoporetech/pomoxis --recursive
git clone https://github.com/nanoporetech/medaka

# While it is possible to install pomoxis and medaka into the same
#   virtual environment, we will install each package into its own
#   environment for simplicity. For more details see the readme for
#   each of the packages.

cd pomoxis && make install && cd ..
cd medaka && make install && cd ..

POMOXIS=${WALKTHROUGH}/pomoxis/venv/bin/activate
MEDAKA=${WALKTHROUGH}/medaka/venv/bin/activate

Creating a Draft Assembly

基于pomoxis，使用miniasm和racon 生成draft assembly，可以从提供的basecalls生成草案程序集。或者，在这一步可以使用canu 。

cd ${WALKTHROUGH}
NPROC=$(nproc)
BASECALLS=data/basecalls.fa
source ${POMOXIS}
mini_assemble -i ${BASECALLS} -o draft_assm -p assm -t ${NPROC}

这将生成draft_assm/assm_final.fa.

mini_assemble脚本有两个有用的选项，此处未使用：

  • 指定-c将对reads 运行porechop去除测序接头;
  • 指定-e 10将在组装之前对前10%长的reads执行error correction（类似于canu的策略）;

这两个步骤都可以在牺牲速度的前提下，提高装配质量。

可以检查assembled contigs 的数量和长度

DRAFT=draft_assm/assm_final.fa
awk '{if(/>/){n=$1}else{print n " " length($0)}}' ${DRAFT}

预期输出是contig长为4,703,280 bp（utg00001c）。

Polishing a Consensus

在执行完Creating a Draft Assembly的所有步骤之后，可以使用medaka的默认模型运行以下命令来生成consensus 。使用从大肠杆菌、酿酒酵母和智人样本中获得的数据对该模型进行训练。

cd ${WALKTHROUGH}
source ${MEDAKA}
CONSENSUS=consensus
DRAFT=draft_assm/assm_final.fa
medaka_consensus -i ${BASECALLS} -d ${DRAFT} -o ${CONSENSUS} -t ${NPROC}

要使用另一个模型polish an assembly，请使用-m 选项指定模型的文件路径。

可使用pomoxis的assess_assembly计算Alignment statistics ：

cd ${WALKTHROUGH}
source ${POMOXIS}
TRUTH=${DATA}/truth.fasta
DRAFT2TRUTH=draft_to_truth
CONSENSUS2TRUTH=${CONSENSUS}_to_truth
echo "Draft assembly"
assess_assembly -i ${DRAFT} -r ${TRUTH} -p ${DRAFT2TRUTH} -t ${NPROC}
echo "Medaka consensus"
assess_assembly -i ${CONSENSUS}/consensus.fasta -r ${TRUTH} -p ${CONSENSUS2TRUTH} -t ${NPROC}

应观察到错误率从0.367%下降到0.070%。

Training a Consensus Network

为了训练一个bespoke network，首先执行包括创建上面的Creating a Draft Assembly的所有步骤。consensus network的最终目的是预测从 basecalls到draft的truth sequence。这需要理解basecalls如何与draft 对齐，以及如何编辑draft 以获得truth。draft充当了basecalls和 truth之间的共同参照系。

basecalls和truth序列与draft align。对于后者，这是分块（chunks）执行的。

cd ${WALKTHROUGH}
DRAFT=draft_assm/assm_final.fa
TRUTH=${DATA}/truth.fasta
source ${POMOXIS}
CHUNKSIZE=100000
CALLS2DRAFT=calls2draft
TRUTH2DRAFT=truth2draft

mini_align -P -m -r ${DRAFT} -i ${BASECALLS} -t ${NPROC} -p ${CALLS2DRAFT}
mini_align -c ${CHUNKSIZE} -P -m -r ${DRAFT} -i ${TRUTH} -t ${NPROC} -p ${TRUTH2DRAFT}

这些raw alignments必须转换为features输入到神经网络。为了减少训练期间的任何IO瓶颈，可以使用 -\\-batch_size 选项将训练数据成批写入HDF5文件。选项-\\-read_fraction用于随机subsample reads，当使用模型进行预测时，其效果是使生成的模型more robust to variations in pileup depth 。

cd ${WALKTHROUGH}
source ${MEDAKA}
REFNAME=utg000001c
TRAINEND=3762624
TRAINFEATURES=train_features.hdf
BATCHSIZE=100
medaka features ${CALLS2DRAFT}.bam ${TRAINFEATURES} --truth ${TRUTH2DRAFT}.bam --threads ${NPROC} --region ${REFNAME}:-${TRAINEND} --batch_size ${BATCHSIZE} --chunk_len 1000 --chunk_ovlp 0

现在一切就绪，可以通过medaka train训练consensus网络：

cd ${WALKTHROUGH}
source ${MEDAKA}
TRAINNAME=training
medaka train ${TRAINFEATURES} --train_name ${TRAINNAME} --epochs 10

根据可用的计算资源，此步骤可能需要一些时间。在训练过程中，定期检查模型，以便在训练中断时可以很容易地恢复训练。在训练结束时，我们有一些输出模型，尤其包括：

  • model.best.hdf5: 在训练集中具有最佳精度的模型
  • model.best.val.hdf5: 在验证集上具有最佳精度的模型

要使用模型，请运行medaka_consensus，使用-m选项指定模型的完整绝对路径：

cd ${WALKTHROUGH}
source ${MEDAKA}
CONSENSUS=consensus_trained
MODEL=${PWD}/${TRAINNAME}/model.best.val.hdf5
medaka_consensus -i ${BASECALLS} -d ${DRAFT} -o ${CONSENSUS} -t ${NPROC} -m ${MODEL}

Training with multipe data types

Medaka支持创建推理网络，处理不同来源的数据，例如两个不同的纳米孔。通过在输入.bam 文件中使用alignment标记来启用。

为了训练一个可以处理多个数据类型的模型，在medaka feature步骤之前，对所有的数据类型都遵循相同的操作。在此之前，合并.bam 文件，同时向文件添加DT字符串标记。例如，要组合R9和R10孔隙的数据：

samtools view r9.bam | awk 'BEGIN{OFS="\t"}{print $0, "DT:Z:r9"}' > r9.sam
samtools view r10.bam | awk 'BEGIN{OFS="\t"}{print $0, "DT:Z:r10"}' > r10.sam
samtools view -H r9.bam | cat - r9.sam r10.sam | samtools view -bS | samtools sort > all_reads.bam
samtools index all_reads.bam

只需要对medaka features附加一个参数，就可以对组合数据类型训练模型：

medaka features all_reads.bam ... --feature_encoder_args dtypes=r9,r10

为了使用多个数据类型训练的模型来创建consensus序列，medaka consensus的输入.bam文件应该将`DT标记适当地添加到所有reads中。如果不将此标记添加到输入文件中，多个数据类型训练的模型将无法起作用。

Automated training pipeline

有了 katuali，现在可以用一个命令从fast5s文件夹开始训练medaka模型：

katuali medaka_train_replicates --keep-going

执行上面的操作将会

  • basecall all multiple runs of data,
  • align所有basecall与reference序列
  • 在reference序列和深度上创建basecalls的subsampled sets
  • 将这些basecalls组装成draft assemblies
  • 为所有assemblies创建medaka training features
  • 在多个重复中训练一个medaka
  • 根据测试数据评估模型。

四、Origin of the draft sequence

Much has been discussed in the popular press on the nature and origin of the expected draft sequence for medaka. Here we give some background of previous workflow recommendations by exploring how medaka works and what this may imply for how it should be used.

Note
建议在运行medaka之前至少运行一轮racon。racon的进一步应用可能是徒劳的。

1. Medaka’s Algorithm

到目前为止，所有的medaka 发行版都在很大程度上利用了一个对草图序列进行read堆积的思想，为神经网络创建输入数据。medaka 的训练方法是创建这样的堆积，并通过将参考序列与相同的草图序列对齐，用真正的base (或 gap)对每列进行注释。

在机器学习领域，一个相当普遍的概念是将模型泛化为用于训练的集合之外的数据，或者是概念上相似的数据，或者是具有特征的新形式。下面的讨论重点介绍了medaka的输入数据可以变得与模型训练时相似或明显不同的方式。

Aligners

在alignment 过程，创建草稿序列的read堆积。medaka性能的第一个关注点是（至少在概念上，如果不是实际的话）结果可能对用于align reads的aligner敏感。此外，对任何一个aligner使用的alignment 参数可能有一个敏感度。

medaka_consensus，使用的aligner 与medaka模型训练期间使用的aligner 相同，即minimap2。然而，medaka consensus将把所有创建的alignment文件.bam 作为输入。

这种对aligner结果的敏感性很大程度上源于gaps是如何计算和记录在alignment中的：具有不同参数化的不同aligner可能导致堆积中不同的gaps结构。这将导致神经网络的输入不同，可能是网络不期望的输入。举个例子，想象一个机器学习算法，训练它从照片中识别不同的苹果品种。当看到一张梨的照片时，这个算法要做什么？它肯定不会给使用这个系统的人一个有意义的答案。medaka的情况可能不像苹果和梨那样严重，尽管关于什么是有效输入的概念仍然存在。

Basecallers

这并不是简单地选择aligner。当然，选择basecaller是一个考虑因素：不同的basecaller将产生具有不同error profiles特征的call（相对于真正序列）。同样，这将影响神经网络输入的性质，并可能导致 basecallers之间缺乏泛化。正是因为这个原因，medaka保守地为各种basecaller和basecaller模型提供了模型。

The draft sequence

也许更微妙的是draft input sequence的影响。draft sequences errors的可变性，例如由使用不同组装方法或程序引起的误差，可能导致不同的神经网络输入。关键点在于，如果草图和reads之间的替换、插入和删除的相对速率与在训练过程中呈现给算法的速率不同，那么性能可能会降低。

Sequencing depth

read堆的结构也受测序深度的影响。这在alignment 查看器（如IGV）中并不十分明显，它们倾向于在视图中squeeze-out insertions。Consider the case where an insertion occurs in one read with respect to the draft sequence or to other reads.。目前在medaka中的实现使用单独的简基创建了一个额外的pileup列。在更高的测序深度，出现更多的来自随机插入事件的插入列：在较高的测序深度，堆积层相对于truth逐渐被延伸。这种影响在用于训练medaka模型的技术中得到了明显的缓解。

Discussion

正是出于这些和类似的考虑，以前关于如何产生草图序列的建议有些严格。这些建议是基于medaka 模型训练的具体方式，预期将产生看起来像苹果而不是梨的神经网络输入！这些建议本质上是悲观的，在向medaka展示一个轻微擦伤的苹果时，推断错误似乎很小。

在twitter conversation中，Devin sounde提供了一个有用的例子，他在使用medaka之前，测试了运行racon不同次数的效果。结果是得到的准确度没有变化。迭代racon过程是否实际产生了不同的结果？当然，如果draft在迭代过程中没有演变，那么medaka 的输出也不会有任何变化。

Using the scripts for assembly benchmark from @nanopore (https://t.co/ezulj7RFN2) I also calculated consensus Q-scores. For this example, I achieved > Q35 pic.twitter.com/XLPdvUFZYr

— Devin Drown (@devindrown) August 7, 2019
需要注意的一点是，对于比Devin使用过的旧的basecaller，assembler 输出的质量较低，随着racon的后续应用，草图质量将得到改善。随着最新的basecallers提高草图的准确度可以使上述问题变得不那么敏感；medaka现在可能对草图的精确性的变化不那么敏感了

How should I create my draft sequence

鉴于上述讨论，问题仍然存在，如何组装草图序列以实现运行medaka后的最佳结果？如前所述，先前的建议是根据medaka模型如何训练而提出的。越来越多的经验证据表明，medaka对草图输入的性质并不过分敏感。Devin的数据确实显示出对总体质量水平仍有一定的依赖性。需要对这个例子进行更多的调试，但总体效果可能是racon 的单个应用比没有racon 的情况更好的原因：草图质量的提高导致medaka 的有效覆盖率，因为更多的reads与改进的草图 align 。应当强调的是，总体质量与前面讨论的草案性质和error profiles之间的概念存在差异。两个草图可以有非常相似的总体质量，但有非常不同的错误概况（error profiles）。

有一种已知的情况是，有可能对草图序列进行过度润色（over-polish）：

Why 4x, have seen good evidence >1 racon corrupts repeats in assembly, especially of diploids. Examples for canu of correctly resolved BACs: 0=74%, 1 = 77%, 2=70%, 3=63%. For miniasm: 0= 21%, 1=60%, 2=52%. Identity does go up but very slightly. I'd vote for 1 round.

— Sergey Koren (@sergekoren) August 6, 2019

其中一个原因是，因为overlap错误的repeat拷贝。reads可以被计算；racon (或 medaka)的迭代应用可能导致这些重复的平均值，而不是保持它们的不同序列。

Recommendations

很难推荐一个single best-practice workflow来获得所有数据集的最佳结果。但是，仍然建议在向medaka提供数据集之前至少运行一次racon。这样做的主要效果是无论使用的那款assembler 进行的组装，都标准化草图的error profile。随着basecall质量和assembler 输出质量的提高，当质量差异和error profiles 之间的区别变得模糊时，这个规范化步骤可能变得多余。

在训练medaka期间，与racon一起使用非默认参数值。上述使用的参数，对reads和草图之间的indel 频率具有相关性。因此，建议在使用 medaka之前将其与racon一起使用：

racon -m 8 -x -6 -g -8 -w 500 ...

在最近的basecallers，来源于这一方法的影响还没有进行过探讨，很可能medaka 对这些数值的变化不太敏感。

Future developments

一个悬而未决的问题是如何完全消除medaka结果对draft sequence的依赖性，至少在某种程度上，所需的只是一个合理准确的 draft。目前正在探索从计算中完全删除草图以及按 face value检查read数据的方法。此外，更充分地利用assemblers输出的方法正在开发中，以避免一些更微妙的问题。作为一种单独的、更简单的可能性，用户可能希望根据他们选择的assembler 的输出来训练medaka模型。
https://github.com/nanoporetech/medaka
https://nanoporetech.github.io/medaka/installation.html

最后编辑于 ：
©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 214,875评论 6赞 496
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 91,569评论 3赞 389
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 160,475评论 0赞 350
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 57,459评论 1赞 288
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 66,537评论 6赞 386
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 50,563评论 1赞 293
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 39,580评论 3赞 414
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 38,326评论 0赞 270
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 44,773评论 1赞 307
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 37,086评论 2赞 330
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 39,252评论 1赞 343
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 34,921评论 5赞 338
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 40,566评论 3赞 322
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 31,190评论 0赞 21
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,435评论 1赞 268
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 47,129评论 2赞 366
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 44,125评论 2赞 352