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

Posted on by

序列的多比对情况大家都懂，因为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

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

​

Comments are closed.