之前有过用二代测序的数据组装植物叶绿体基因组昆虫线粒体的经历，用的是单位的超算（Linux系统）。

这里的二代测序数据是全基因组的浅层测序数据，因为叶绿体和线粒体是多拷贝的，一般浅层测序数据就可以组装出完整的叶绿体和线粒体基因组。我的单个样本（昆虫）测序数据大小是4G， 仅供参考。

用到的软件为Getorganelle和Mitofinder，这里介绍Mitofinder。

Mitofinder

官网：https://github.com/RemiAllio/MitoFinder

下面的教程基本来自于对官网教程的翻译，如有需要可以去看官网原文。

Mitofinder is a pipeline to assemble mitochondrail genomes and annotate mitochondrial genes from trimmed read sequencing data.

Mitofinder is also designed to find and annotate mitochpondrail sequences in existing genomic assemblies.

Mitofinder用来组装和注释线粒体基因组。

Mitofinder在整个流程中会调用的其他程序包括：

用于BLAST的：

用于组装的：MEGAHIT、MetaSPAdes、IDBA-UD

用于tRNA 注释的：MiTFi、tRNA-scan、ARWEN

Mitofinder的安装

Mitofiinder最简单的安装方法就是用conda安装。

Mitofinder是在python2.7下编写的，所以安装的时候建议用conda新建一个python2.7的环境，方法参考之前的教程：

https://www.jianshu.com/p/3ef9e6041dee

当然，用coonda创建环境和安装Mitofinder的前提是你已经安装好了miniiconda

miniconda的安装可以参考之前的教程：

https://www.jianshu.com/p/9dc419c33f42

首先创建一个名为Mitofinder，安装2.7版本的python的环境

conda create -n Mitofinder python=2.7

激活环境

conda activate Mitofinder

conda安装mitofinder

conda install -c bioconda mitofinder

组装并注释线粒体基因组

用Mitofinder组装和注释线粒体基因组很简单，只需要一行命令。

如果是双向测序的原始数据，一个样本两个测序原始数据文件(left_reads.fastq.gz和right_reads.fastq.gz)，用如下命令：

mitofinder -j [seqid] -1 [left_reads.fastq.gz] -2 [right_reads.fastq.gz] -r [genbank_reference.gb] -o [genetic_code] -p [threads] -m [memory]

如果是单向测序的原始数据，一个样本一个测序原始数据文件(SE_reads.fastq.gz)，用如下命令：

mitofinder -j [seqid] -s [SE_reads.fastq.gz] -r [genbank_reference.gb] -o [genetic_code] -p [threads] -m [memory]

如果已经组装好了线粒体基因组，只需要测序，那么组装好的线粒体基因组序列作为输入文件，要求为fasta文件（assembly.fasta），用如下命令：

mitofinder -j [seqid] -a [assembly.fasta] -r [genbank_reference.gb] -o [genetic_code] -p [threads] -m [memory]

下面是我自己的例子：

~/Programfile/MitoFinder/mitofinder -j  XXXX -1 XXXX.fastq.gz -2 XXXX.fastq.gz -r reference.gb -o 5 -p 5 -m 10

各项参数的意义如下：

~/Programfile/MitoFinder/mitofinder：给出mitofinder所在的路径

-j  任务的id号，输出的结果文件也用id号命名

-r 近缘物种的已经注释的线粒体基因组文件，要求为.gb文件，可以自行去NCBI网站下载。

下载的方法参考之前的教程：https://www.jianshu.com/p/910c92a0d03b

-o 指遗传密码类型，5代表无脊椎动物线粒体的遗传密码

-p 允许Mitofinder在运行时使用的最大线程数

-m 在组装过程中（运行MEGAHIT 或 MetaSPAdes)允许使用的计算机最大存储量

如果我们有很多测序数据需要组装，那么就可以写一个循环，写在任务脚本run.sh文件里，然后直接运行run.sh文件就可以了。我用的循环如下（惭愧，不是我写的，组上老师写的，这里挪用一下）

解释下上半部分：

-cwd #指定当前路径为工作目录，sge的日志会输出到当前路径

-S #指定远程计算节点的shell路径

-l #指定资源请求，多个请求用逗号（，）隔开

这一部分只针对我们单位系统的超算，不适用于大家，可以不看。

用到的循环如下：

for file1 in $(ls *.R1.fastq.gz)

  do

  for file2 in $(ls *.R2.fastq.gz)

  do

  if [ ${file1:0:4} == ${file2:0:4} ]; then

~/Programfile/MitoFinder/mitofinder -j ${file2:0:4} -1 $file1 -2 $file2 -r reference.gb -o 5 -p 5 -m 10

  fi

  done

  done

简单解释一下这个循环：

我这里用到的每个样本的两个测序原始数据文件命名为：XXXX.R1.fastq.gz和XXXX.R2.fastq.gz，

XXXX是4个字符，把所有的输入文件都放在一个路径下，如果XXXX.R1.fastq.gz和XXXX.R2.fastq.gz具有相同的XXXX（在循环中为${file1:0:4}

== ${file2:0:4}），就对它们执行组装命令，输出文件命名为XXXXout（在循环中为 ${file2:0:4}"out"）。

然后运行任务脚本就可以了：

qsub run.sh

qsub只是我们单位超算运行任务脚本的命令，大家根据自己的linux系统输入运行脚本的命令。

如果已经用别的软件组装好了线粒体基因组序列，只是用mitofinder进行注释，在序列很多的情况下也可以用下面的循环：

for filename in *.fasta;

do

~/Programfile/MitoFinder/mitofinder -j ${filename#*.} -a ${filename} -r reference.gb -o 5；

done

j后面代表的是任务的名称，也是后面result文件的名称，${filename#*.} 意思是取输入文件最后一个“.”前面的字符，如果输入文件是F410.fasta. 这个任务的名称就是F410