序列的多比对情况大家都懂，因为NGS时代，序列都很短，也就是50-250bp范围，而且参考基因组本来就是会有很多低复杂度区域，那么我们的reads比对到参考基因组的多个区域，就很好理解了。
最近有粉丝咨询，因为有些比对工具为了保证输入多少reads就输出多少条比对记录，所以会随机挑选一个最好的比对，然后问我是不是hisat2也会对多比对的reads随机输出一条吗？我觉得有必要帮忙探索一下，分享这个过程。
其实是可以设置参数的，需要看说明书，如下：
代码语言：javascript复制 Reporting:
(default) look for multiple alignments, report best, with MAPQ
OR
-k <int> report up to <int> alns per read; MAPQ not meaningful
OR
-a/--all report all alignments; very slow, MAPQ not meaningful
不过考虑到大部分人对软件的使用很初级，应该是默认参数，我也这样演示：
代码语言：javascript复制hisat2=/trainee/jmzeng/tools//hisat2-2.0.0-beta/hisat2
fasta=/trainee/jmzeng/Probe_seqfasta/lncRNA/human/GPL15314_seq2fa.fasta
sample=GPL15314
index=/teach/database/index/hisat/hg38/genome
$hisat2 -f $fasta -x $index -S ${sample}.sam
就是芯片探针的序列，比对到参考基因组。
首先看我们的比对日志
输入的fasta序列是60699 个reads，有 54578 (89.92%)条reads都是精准匹配到参考基因组的唯一位置。
代码语言：javascript复制60699 reads; of these:
60699 (100.00%) were unpaired; of these:
519 (0.86%) aligned 0 times
54578 (89.92%) aligned exactly 1 time
5602 (9.23%) aligned >1 times
99.14% overall alignment rate
因为它是靠是否含有 NH:i:1 来判断是否是唯一比对。
hisat2认为是唯一比对的其实也有可能是多比对
下面的这个60bp长度的探针，因为标记了 NH:i:1，所以认为是唯一比对，其成功比对到了参考基因组的chr1的23527046坐标，而且整个比对的sam文件里面的确只能找到这一个记录。 这个时候，肯定是认为其是唯一比对啦！
代码语言：javascript复制GPL15314:ASHG19A3A000433 0 chr1 23527046 255 25M334N35M * 0 0 CAAGCCCAGGAGGATCAGCATCAAGGCAGCCCCTTTGGCTGGATCTTCAAGGCCTGGTCA IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII AS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:60 YT:Z:UU XS:A: NH:i:1
但是我把该序列拿到blat工具去看了看：https://genome.ucsc.edu/cgi-bin/hgBlat
结果是其在参考基因组的两个区域都是完美匹配，但是呢，这样的情况其实是参考基因组本身的问题，包含了那些不是染色体的片段的碱基序列。
image-20191201170907300
因为我们的hisat2工具使用的参考基因组里面并没有blat这个工具里面的那个染色体，所以出现冲突也不意外。
再看看另外一个比对到6个地方的例子
可以看到这个序列，被比对了6次，所以记录了 NH:i:6 。
代码语言：javascript复制4249:GPL15314:ASHG19A3A004249 0 chr6_GL000252v2_alt 3219967 1 25M85N35M * 0 0 GCCCGGCCTGGCGGAGGTGATGCTGGAGGGACGCCCGGGGAGACCGTACGTCACTGCTCT IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII AS:i:0 ZS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:60 YT:Z:UU XS:A: NH:i:6
4249:GPL15314:ASHG19A3A004249 256 chr6_GL000254v2_alt 3314228 1 25M85N35M * 0 0 GCCCGGCCTGGCGGAGGTGATGCTGGAGGGACGCCCGGGGAGACCGTACGTCACTGCTCT IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII AS:i:0 ZS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:60 YT:Z:UU XS:A: NH:i:6
4249:GPL15314:ASHG19A3A004249 256 chr6_GL000255v2_alt 3228162 1 25M85N35M * 0 0 GCCCGGCCTGGCGGAGGTGATGCTGGAGGGACGCCCGGGGAGACCGTACGTCACTGCTCT IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII AS:i:0 ZS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:60 YT:Z:UU XS:A: NH:i:6
4249:GPL15314:ASHG19A3A004249 256 chr6_GL000256v2_alt 3273381 1 25M85N35M * 0 0 GCCCGGCCTGGCGGAGGTGATGCTGGAGGGACGCCCGGGGAGACCGTACGTCACTGCTCT IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII AS:i:0 ZS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:60 YT:Z:UU XS:A: NH:i:6
4249:GPL15314:ASHG19A3A004249 256 chr6 31972192 1 25M85N35M * 0 0 GCCCGGCCTGGCGGAGGTGATGCTGGAGGGACGCCCGGGGAGACCGTACGTCACTGCTCT IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII AS:i:0 ZS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:60 YT:Z:UU XS:A: NH:i:6
4249:GPL15314:ASHG19A3A004249 256 chr6_GL000251v2_alt 3449619 1 25M85N35M * 0 0 GCCCGGCCTGGCGGAGGTGATGCTGGAGGGACGCCCGGGGAGACCGTACGTCACTGCTCT IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII AS:i:0 ZS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:60 YT:Z:UU XS:A: NH:i:6
这个时候可以看到， 其实这个序列也是在hisat2使用的参考基因组里面了比对到了多条非染色体片段。
同样的，我们也是在blat工具检查看看；
image-20191201171304037
实际上，我们并不想看到这样的事情发生，我们只需要染色体序列即可，并不需要那么多的非染色体片段。就是chr6_GL000251v2_alt这样的不需要出现在我们的参考基因组。
我们看看比对次数最多的
有很多序列都是可以比对10次， 我随便找了一个，如下：
代码语言：javascript复制TAACATTGGGAGAAATAGCCAGCTGAATCTGTAACTCAACAGAAACAAGTGATCCATATA
在blat结果是：
image-20191201220735165
hisat2默认是允许错配的
比如下面这个序列，虽然是60M，但是里面的MD:Z:5G54表明这个序列的第5个位置碱基错配。
代码语言：javascript复制chr10 78562138 0 60M * 0 0
CAGAATGAGTGGGAGGAGAGAAATGCATTGCTCCAAGTCCAAGAAAATGATCATTATCAA
AS:i:-6 ZS:i:-6 XN:i:0 XM:i:1 XO:i:0 XG:i:0 NM:i:1 MD:Z:5G54 YT:Z:UU NH:i:1
同样的，blat也可以看到结果：
image-20191201223222866
还有更多其它形式的错配，自己可以慢慢查看：
代码语言：javascript复制chr5 57973331 0 60M * 0 0
TTGGTCGTTTAACAGCAACCCCCTCCCCTATAATAATCAGAACCATTCACTTCAGCTAAT
AS:i:-12 ZS:i:-12 XN:i:0 XM:i:2 XO:i:0 XG:i:0 NM:i:2 MD:Z:5T32T21
MD:Z:5T32T21 代表的错配如下：
代码语言：javascript复制00000001 ttggtcgtttaacagcaaccccctcccctataataatcagaaccattcac 00000050
>>>>>>>> ||||| |||||||||||||||||||||||||||||||| ||||||||||| >>>>>>>>
57973331 ttggttgtttaacagcaaccccctcccctataataatctgaaccattcac 57973380
00000051 ttcagctaat 00000060
>>>>>>>> |||||||||| >>>>>>>>
57973381 ttcagctaat 57973390
可以看到，是两个碱基错配。
