BWA关于生成MEM问题

BWA中在查找MEM集合概述

MEM（Maximal Exact Match）：
是一个 在给定文本中无法进一步向左或向右扩展 的匹配子串。

SMEM（Supermaximal Exact Match）：
是一种特殊的 MEM，它在 查询串中不被其他 MEM 覆盖。

• 第一步会根据QUERY序列起始开始，在当前查找匹配会进行两端扩展匹配。扩展匹配过程，会先正向匹配，每有MEM就会保存作为候选MEM，然后在这些候选MEM基础上反向扩展，每次扩展后必须。
• 第二步会在第一步得到的候选MEM上查找较长的，从中点进行左右扩展，然后加入到候选MEM中。
• 第三步从当前位置向右查找一个不可再扩展的MEM加入候选，更新当前位置再查找。

我在BWA源码中进行更改：

• 只保留BWA查找MEM的第一步，且扩展时只进行正向扩展，向右扩展时只有当前位置不可再向右扩展时才加入候选MEM。
  这是更改后与原始BWA比对输出：第一张是视作单端read比对结果，第二张是双端结果（read原本是双端）主要查看第5列的CIGAR输出（他是SW评分计算后生成，和SW评分效果一致，M代表完全匹配，S代表软截取，H代表硬截取，前面跟数字代表该结果碱基数量。如：72M代表72个完全匹配）。
  1.明显更改后的匹配结果评分更高，完全匹配的数量更多。
  2.原始BWA对于一个READ可能会输出多个结果，会通过第二个参数作为FLAG进行标识，关于flag可以在网站查阅，但评分一样的结果不会多次输出，即在BWA比对过程中的种子链不会都输出。https://broadinstitute.github.io/picard/explain-flags.html

注意

• 该测试结果可能并没有实际价值，因为SW评分只能判断某条读段比对的好坏，不能用于比较不同比对策略或者比对工具的优劣基准，应该选择其他的测试基准。

• 关于正向扩展时，如果选择只在read上对每个碱基位置进行向右扩展，以此得到的精准匹配可能不是最大精准匹配（MEM），向右扩展只能保证在右端是精准匹配的，不能确定在TARGET链上向左是不可扩展的，这不符合MEM的定义。
  比如TARGET ：AAGACCTCCT，READ：ACTCCT中
  READ3号位匹配到的TCCT就不是MEM，因为在TARGET链上比对到的TCCT位置7往前扩展的C碱基仍然是可以和READ碱基比对上的，TCCT比对被覆盖于TARGET 位置6的CTCCT。
  MEM的是否可扩展是TARGET，SEMS同时针对于TARGET和READ。

参考碱基序列 https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/000/395/GCF_000000395.1_GRCz11/GCF_000000395.1_GRCz11_genomic.fna.gz
read碱基序列 SRR390728（prefetch 下载截取前10条40行）

测试结果

GPU代码存在BUG，未能输出计算结果：