【4.1】minimap2序列比对

ngs_topic 阅读量：次

Minimap2是李恒大牛在2018年开发的针对于三代测序数据进行比对的工具，minimap2的优势是速度快，而且听说比对的结果也比较不错，不知道甩samtools几条街；缺点呢，就是耗费内存。

由于二代测序序列长度短（~300 bp），错误率低（<0.02%），其测序错误主要为替换错误，可通过哈希查找或数据压缩等方法直接查询每条序列在基因组中的准确比对区域，得到比对结果。而三代测序序列一方面错误率较高（~15%），且错误类型主要为插入或删除错误，难以通过文本查找直接找到序列在基因组中的比对区域，需要采用针对性的查询策略找到比对区域。另一方面，在比对阶段，由于三代序列读段长（平均长度超过10 kbp），直接采用全局比对（Needleman-Wunsch，NW）或局部比对（Smith-Waterman，SW）比对算法需构建L×L大小的得分矩阵（其中L表示序列长度），计算空间复杂度较高。基于以上两方面分析，目前绝大多数针对二代测序的比对方法不适合处理三代测序数据。因此， 需要开发针对三代测序数据的序列比对算法。

minimap2的主要思想是：首先将基因组序列的minimizer存储在哈希表中(minimizer指一段序列内最小哈希值的种子)；然后对于每一条待比对序列， 找到待比对序列所有的minimizer，通过哈希表找出其在基因组中的位置， 并利用chaining算法寻找待比对区域；最后将非种子区域用动态规划算法进行比对，得到比对结果。minimap2方法只对最小哈希值的种子进行存储，可有效降低时间复杂度。

一、安装

cd /data/software
git clone https://github.com/lh3/minimap2
cd minimap2 && make

二、用法

minimap2有多种比对功能，处理支持nanopre数据之外，还支持pacbio数据。比对模式可以是reads与reads之前比对，reads与基因组比对，基因组与基因组比对，以及短序列与基因组的比对，不同的比对有不同的作用，千万不能设置错误了。下面我们拿具体案例演示一下。 minimap2最常用的功能就是将测序数据比对到基因组上，这个过程与bwa比对类似，需要先建立索引，然后比对，最终得到sam格式的比对结果，如果熟悉bwa比对，那么这个操作就非常容易。 minimap2的比对输入文件为测序的reads，一般为fastq或者fasta格式，参考基因组，一般为fasta格式。minimap2可以输出paf格式以及sam格式，默认为paf格式。

第一步：建立索引

minimap2 -d co92.min co92.fna

第二步：比对

minimap2 -ax map-ont co92.min ../4.filter/clean.filtlong.fq.gz >s1037.sam

主要分成五大类，索引（Indexing），回帖（Mapping），比对(Alignment),输入/输出（Input/Output）,预设值（Preset）。

-x ：非常中要的一个选项，软件预测的一些值，针对不同的数据选择不同的值
map-pb/map-ont: pb或者ont数据与参考序列比对；
ava-pb/ava-ont: 寻找pd数据或者ont数据之间的overlap关系；
asm5/asm10/asm20: 拼接结果与参考序列进行比对，适合~0.1/1/5% 序列分歧度；
splice: 长reads的切割比对
sr: 短reads比对

-d :创建索引文件名
-a   ：指定输出格式为sa格式，默认为PAF
-Q   ：sam文件中不输出碱基质量    
-R ：reads Group信息，与bwa比对中的-R一致
-t：线程数，默认为3

其他用法：

# long sequences against a reference genome
./minimap2 -a test/MT-human.fa test/MT-orang.fa > test.sam

# create an index first and then map
./minimap2 -x map-ont -d MT-human-ont.mmi test/MT-human.fa
./minimap2 -a MT-human-ont.mmi test/MT-orang.fa > test.sam

# use presets (no test data)
./minimap2 -ax map-pb ref.fa pacbio.fq.gz > aln.sam       # PacBio CLR genomic reads
./minimap2 -ax map-ont ref.fa ont.fq.gz > aln.sam         # Oxford Nanopore genomic reads
./minimap2 -ax map-hifi ref.fa pacbio-ccs.fq.gz > aln.sam # PacBio HiFi/CCS genomic reads (v2.19 or later)
./minimap2 -ax asm20 ref.fa pacbio-ccs.fq.gz > aln.sam    # PacBio HiFi/CCS genomic reads (v2.18 or earlier)
./minimap2 -ax sr ref.fa read1.fa read2.fa > aln.sam      # short genomic paired-end reads
./minimap2 -ax splice ref.fa rna-reads.fa > aln.sam       # spliced long reads (strand unknown)
./minimap2 -ax splice -uf -k14 ref.fa reads.fa > aln.sam  # noisy Nanopore Direct RNA-seq
./minimap2 -ax splice:hq -uf ref.fa query.fa > aln.sam    # Final PacBio Iso-seq or traditional cDNA
./minimap2 -ax splice --junc-bed anno.bed12 ref.fa query.fa > aln.sam  # prioritize on annotated junctions
./minimap2 -cx asm5 asm1.fa asm2.fa > aln.paf             # intra-species asm-to-asm alignment
./minimap2 -x ava-pb reads.fa reads.fa > overlaps.paf     # PacBio read overlap
./minimap2 -x ava-ont reads.fa reads.fa > overlaps.paf    # Nanopore read overlap
# man page for detailed command line options
man ./minimap2.1

三、案例

a. 长序列比对

minimap2 -ax map-pb  ref.fa pacbio-reads.fq > aln.sam   # for PacBio CLR reads

minimap2 -ax map-ont ref.fa ont-reads.fq > aln.sam      # for Oxford Nanopore reads

b. 给参考序列建立索引

minimap2 -d ref.mmi ref.fa                     # indexing
minimap2 -a ref.mmi reads.fq > alignment.sam   # alignment

minimap2比对文件处理

利用minimap2将nanopore测序数据与参考序列比对得到的sam之后，需要对sam文件行处理，主要包括排序与建立索引。排序主要是依据参考序列位点进行排序。

#samtools处理
samtools sort -@ 4 -O bam -o s0137.sorted.bam s1037.sam
samtools index s0137.sorted.bam
samtools faidx co92.fna

参考资料

药企，独角兽，苏州。团队长期招人，感兴趣的都可以发邮件聊聊：tiehan@sina.cn
个人公众号，比较懒，很少更新，可以在上面提问题，如果回复不及时，可发邮件给我： tiehan@sina.cn