hisat2会对多比对的reads随机输出一条吗？

序列的多比对情况大家都懂，因为NGS时代，序列都很短，也就是50-250bp范围，而且参考基因组本来就是会有很多低复杂度区域，那么我们的reads比对到参考基因组的多个区域，就很好理解了。

最近有粉丝咨询，因为有些比对工具为了保证输入多少reads就输出多少条比对记录，所以会随机挑选一个最好的比对，然后问我是不是hisat2也会对多比对的reads随机输出一条吗？我觉得有必要帮忙探索一下，分享这个过程。

其实是可以设置参数的，需要看说明书，如下：

 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


不过考虑到大部分人对软件的使用很初级，应该是默认参数，我也这样演示：

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都是精准匹配到参考基因组的唯一位置。

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文件里面的确只能找到这一个记录。 这个时候，肯定是认为其是唯一比对啦！

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 。

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次， 我随便找了一个，如下：

TAACATTGGGAGAAATAGCCAGCTGAATCTGTAACTCAACAGAAACAAGTGATCCATATA


在blat结果是：

image-20191201220735165

hisat2默认是允许错配的

比如下面这个序列，虽然是60M，但是里面的MD:Z:5G54表明这个序列的第5个位置碱基错配。

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

还有更多其它形式的错配，自己可以慢慢查看：

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 代表的错配如下：

00000001 ttggtcgtttaacagcaaccccctcccctataataatcagaaccattcac 00000050
>>>>>>>> ||||| |||||||||||||||||||||||||||||||| ||||||||||| >>>>>>>>
57973331 ttggttgtttaacagcaaccccctcccctataataatctgaaccattcac 57973380

00000051 ttcagctaat 00000060
>>>>>>>> |||||||||| >>>>>>>>
57973381 ttcagctaat 57973390


可以看到，是两个碱基错配。