4. Using MiToHiFi to Assemble Mitochondrial Genome from HiFi Data

现在的 HiFi reads 的长度普遍有10~15kb，几乎是正常哺乳动物的线粒体基因组的长度了。所以利用测得的 hifi reads 组装出来线粒体基因组是比较容易的，之前我的办法就是直接从 hifiasm 的结果中找到装出来的环 contig （也就是序列 ID 末尾是c的contig，c指的就是cyclic，末尾是l指的是线性序列，line）。

但现在也有一些专门用来组装 mitogenome 的工具，其中一个就是MiToHiFi：‣，本文就介绍如何使用 MiToHiFi 组装线粒体基因组并注释。

软件文章：

MitoHiFi: a python pipeline for mitochondrial genome assembly from PacBio high fidelity reads - BMC Bioinformatics

安装软件

提前安装好一个依赖的软件 MitoFinder：https://github.com/RemiAllio/MitoFinder，会用来注释组装好的基因组，将MitoFinder 添加到环境变量中。

采用官方给出的 conda 安装步骤安装，新建一个conda 的虚拟环境。
```
#Clone MitoHiFi git repo
git clone <https://github.com/marcelauliano/MitoHiFi.git>

#create a conda environment with our yml file that is inside MitoHiFi/environment
conda env create -n mitohifi_env -f MitoHiFi/environment/mitohifi_env.yml
```
在建好的虚拟环境中执行 MitoHiFi 路径下的mitohifi.py即可。亲测 v2.2 版本比较稳定，我试了 v3.0 和 v3.2 都会报错，测试文件都跑不出来。
准备数据

（1）从 NCBI 上找一个mitogenome 的 Refseq，将其 fasta 文件和 Genbenk 注释文件都下载下来，作为模板。如果同物种没有 Ref mitogenome，也可以找近缘物种的 Ref mitogenome 做参考。

（2）准备好 HiFi fasta 格式的 reads 文件，或者是上一个流程获得的 contig fasta 文件 pri.hic.asm.hic.p_ctg.fa。软件将从原始 reads 或者 contigs 中组装mitogenome。

3.HiFi + HiC 数据通过 hifiasm 组装嵌合体基因组
运行软件

我是直接输入 contigs.fa 进行组装，命令如下：
```
conda activate mitohifi_env
python /work_path/software/MitoHiFi-2.2/mitohifi.py -c pri.hic.asm.hic.p_ctg.fa -f NC_023457.1.fa -g NC_023457.1.gb -t 32 -o 2
```
-t 线程数， -o为NCBI中不同生物的线粒体遗传密码，2指的就是脊椎动物的密码表。
环状基因组 restart

软件运行成功结束后，同目录下的 final_mitogenome.fasta，即为组装的结果，但是组装好的基因组与 NCBI refseq 的起始位置不太一样，NCBI 的线粒体基因组，基本上都是以 tRNA-Phe 基因组为起始。所以我们要把组装结果从新划定起始位置。

可以将 final_mitogenome.fasta 在NCBI 上做一个在线的 blastn，网址如下，选择文件提交稍等就好。

NCBI - WWW Error Blocked Diagnostic

在比对好的结果里，查看 Alignment 模块，如果看到组装结果和 ref 并不是起止相同的，而是对上了两截，那么就需要调整，具体根据 Alignment 对应到 Ref 起始的 final_mitogenome.fasta 位置做 restart。假如 final_mitogenome.fasta 的432 位碱基需要作为起始，那可以用seqkit restart完成：
```
seqkit restart -i 432 final_mitogenome.fasta > final_mitogenome.restart.fasta
```
restart 完成后，再进行 blastn，应该是与 Ref 起止完美对应的，如下图所
mitofinder 注释

经过 restart 的 mitogenome：final_mitogenome.restart.fasta，再跑一次 mitofinder 进行注释，命令如下：
```
mitofinder --max-contig-size 82220 -j final_mitogenome.restart.annotation -a final_mitogenome.restart.fasta -r NC_023457.1.gb -o 2 -p 20 --circular-size 8000
```
-r 就是 ref 的注释文件，-o 与 mitohifi 的-o参数含义一样， -p 为线程数，--max-contig-size 指最大的 contigs 的长度，默认为 2.5 Kb，--circular-size 检查的循环次数，最后两个参数为 mitoHiFi 调用 mitofinder 时设置的，单独跑不用设置也可以。

最终的结果为：
```
final_mitogenome.restart.annotation/final_mitogenome.restart.annotation_MitoFinder_mitfi_Final_Results/final_mitogenome.restart.annotation_mtDNA_contig.gb
```
此文件上传 NCBI 时会用到。也可到在线平台 OGDRAW 可视化。

MPI-MP CHLOROBOX - OGDRAW
后续还可以把线粒体基因组 final_mitogenome.restart.fasta 与 pri.hic.asm.hic.p_ctg.fa 做一个 blastn，找出时线粒体序列的比较小的 contigs，并去重。