LASTZ说明文档(v1.04.15, Aug27 2021 build)
本文参考自LASTZ发行版中的说明文档，并根据实际情况进行了一些修改。原文档较长，此处仅翻译了文档前面的一小部分。
介绍
本文档用于介绍LASTZ序列比对程序的安装和使用。LASTZ是BLASTZ的替代实现，并对BLASTZ的命令行语法保持前向兼容。这意味着LASTZ支持BLASTZ的所有选项，并有一些额外的功能，并且可能产生略有不同的比对结果。
- LASTZ：一个工具，用于（1）对两条DNA序列进行比对，（2）自动推断适当的评分参数
- 平台：本工具在macOS X上开发，但可以几乎不经过修改而运行在其他Unix或Linux平台上。LASTZ用C语言编写，用gcc编译，如果用户想使用其他编译器，仅需修改Makefile的内容。一些辅助工具以python编写，仅需要典型python配置中可用的模块既可正常运行（也就是说，只要安装了python就能使用这些工具，不需要额外安装一些软件包）。
- 作者：Bob Harris <rsharris@bx.psu.edu>
- 更新时间：2021年8月
- 邮件列表：http://lists.bx.psu.edu/listinfo/lastz-users
可用性
LASTZ已经在github上开源，项目地址：https://github.com/lastz/lastz/
宾夕法尼亚州立大学的Miller实验室提供了一个包含LASTZ源代码的打包存档。注：1.04.15是Miller实验室网站上提供的最新版本。
安装
如果您以打包存档的方式收到了发行版，请按照适合您的计算机的方式解开包装。这将会得到一个目录<somepath>/lastz‑distrib‑X.XX.XX，其中包含一个名为src的子目录，以及一些其他的东西。为了方便使用，你也可以对目录进行重命名，将版本号-X.XX.XX删去。
在构建或安装本程序之前，你可能需要告诉安装器在哪里放置可执行文件，这可以通过设置shell环境变量$LASTZ_INSTALL，或在make-include.mak文件中修改installDir实现。并且，记得将可执行文件的安装目录放进$PATH环境变量中。
之后构建LASTZ可执行文件。在shell终端（可以是bash，zsh或其他）中依次执行下列的指令（对于Solaris系统的用户，使用gmake指令代替make指令）。这将构建出两个可执行文件，分别是lastz和lastz_D，并将这两个可执行文件复制到installDir目录下面。
cd <somepath>/lastz-distrib-X.XX.XX/src
make
make install
这两个可执行文件基本上算是一个程序，唯一的差别在于lastz使用整数得分，而lastz_D使用浮点数得分。
构建过程不应报告任何警告或错误。因此，设置了makefile，以便将警告视为错误并停止构建。如果遇到这种情况，则可以使用下列的指令进行构建
make -f Makefile.warnings
这应该允许构建完成，同时仍报告警告。您需要确定警告是否表明确实有问题。通常情况下是不会有问题的，但仍然需要报告给作者。
在发行版的打包存档中应该还包括一个简单的自检程序，因此您可以测试构建是否成功。要运行它，请输入以下命令：
make test
如果自检程序通过，将不会有任何输出。否则您将看到预期的输出和构建的输出之间的差异，以及一行如下所示的信息：
make: *** [test] Error 1
构建选项
可以构建一个额外的可执行文件（lastz_32），以处理大于2GB碱基的基因组。有关详细信息，请参见比对到全基因组一节。
可以构建任何可执行文件以允许相邻的索引（默认情况下，不允许使用这些索引）。有关详细信息，请参见相邻索引一节。
处理阶段和术语的概述
LASTZ被设计用于处理一条或多条序列（我们将其称之为 目标（target） ），并将多条 查询（query） 序列与之进行比对。总体上看，本软件的工作流就像是一条流水线（pipeline）：一个阶段的输出会成为下一个阶段的输入。用户可以通过命令行选项跳过大部分的阶段；任何被跳过的阶段都会将其输入原封不动地传递到下一阶段。其中有两个阶段是特殊的，分别是分数推理阶段和插值阶段，这两个阶段中各自执行了一个微型版本的流水线。
注意，下面的讨论是一个概括，旨在描述LASTZ操作的基本思想。有许多例外情况取决于指定的特定选项。
这些阶段包括
- 从目标文件中读取序列
- 推测得分参数（如有必要）
- 为目标构建一个种子单词位置表
- 对于查询文件中的每条序列：
- 寻找种子
- 将种子进行无间隔的扩展以得到HSP
- 链接HSP
- 将HSP进行有间隔的扩展以得到局部比对
- 基于特定标准过滤序列比对结果
- 在序列比对之间进行插值
- 打印输出结果
- 对查询序列的反向互补序列重复上述过程
通常的流程如下所示（尽管这些步骤中的大多数是可选的，并且-‑anyornone等一些设置可能会影响处理顺序）。
- 我们首先将目标序列读入内存，并使用它来构建种子单词位置表（seed word position table），该表将允许我们快速将目标中的任何词映射到它出现的所有位置。（为了本讨论的目的，您可以将一个 单词（word） 看作一个12聚体的DNA）。
- 然后我们依次读取每个查询序列，或多或少地独立处理它们。我们检查查询序列中从每个碱基开始的单词，并使用位置表在目标中查找匹配项，后者被我们称为 种子（seed） 。
- 将种子进行扩展，从而得到更长的匹配，后者被我们称为 HSPs（高得分片段配对，high-scoring segment pairs） ，之后我们基于分数对HSP进行过滤。
- HSP被链接到共线性比对的最高得分集合中，然后减少到单个位置，后者被我们称为 锚点（anchors） 。之后，这些锚点将扩展到局部比对（可能包含间隙），并再次按分数进行过滤，然后进行后端过滤，以丢弃不符合某些特征指定标准的序列比对区块。
- 然后我们进行插值，在序列比对区块之间的洞（holes）中以更高的灵敏度重复整个过程。
- 最后，我们将比对信息写入文件。然后，我们使用查询序列的反向互补序列重复上述步骤，之后开启对下一个查询序列的处理过程。
通常不执行评分推断阶段。一般来说，它仅在获取新物种序列时使用，以便为这些物种的后续比对创建评分文件。
示例
对于那些渴望尝试的人，这里有一些示例来帮助您开始。详细参考资料从下一节开始。
一、比较人类染色体和鸡染色体
与LASTZ的默认值相比，使用较低的灵敏度水平通常是足够的。例如，为了比较两个完整的染色体，即使是人类和鸡这样遥远的物种，即使在非常低的灵敏度设置下，序列比对的格局也很明显。这可以大大加快比对的过程。
本例比较了人类4号染色体和鸡4号染色体。这些序列可以在UCSC基因组浏览器的下载页面中找到，长度分别为191Mbp和94Mbp。要运行这些序列的快速低灵敏度对齐，请使用以下命令：
lastz hg18.chr4.fa galGal3.chr4.fa \
--notransition --step=20 --nogapped \
--format=maf > hg18_4_vs_galGal3_4.maf
在一台2GHz的工作站上运行上述指令，只消耗了大约400Mb的RAM，耗时约为两分半钟。下图展示了结果，其中这张图使用参数--format=rdotplot选项，并使用R statistical软件包绘制。
- 注：上图由笔者电脑生成，在运行
lastz时额外添加了参数--rdotplot=hg18_4_vs_galGal3_4.plot，以获得一个R.plot文件。这个文件可以使用下列R语言代码进行可视化dots = read.table("hg18_4_vs_galGal3_4.plot",header=T) plot(dots,type="l")
使用--notransition降低了对种子的敏感性，从而降低了运行时间。--step=20同样降低了对种子的敏感性，不仅降低了运行时间，还降低了内存消耗。--nogapped消除了对有间隔的序列比对（gapped alignment）的计算。如果使用默认选项，将会消耗1.3Gb的RAM，并需要消耗大约4.5小时的运行时间。（由于笔者精力所限，没有对这一方法进行尝试）
二、将鸟枪法测序结果比对到人类基因组上
相近物种的短读长reads的映射需要一些与LASTZ的默认值非常不同的参数。本例将一组模拟的灵长类鸟枪法测序reads与人类21号染色体进行比较。人类的21号染色体可在UCSC基因组浏览器的下载部分找到（约47Mbp）。通过从黑猩猩chr21中提取长度为60bp的序列，对其进行轻度突变（包括短间隔，short gaps），然后将其截断为50bp，从而生成了10000个模拟reads（这些都包含LASTZ发行版的test_data/fake_chimp_reads.2bit中）。
要查看这些reads映射到人类染色体的位置，请使用以下命令：
lastz hg18.chr21.fa[unmask] fake_chimp_reads.2bit \
--step=10 --seed=match12 --notransition --exact=20 --noytrim \
--match=1,5 --ambiguous=n \
--filter=coverage:90 --filter=identity:95 \
--format=general:name1,start1,length1,name2,strand2 \
> hg18_21_vs_reads.dat
将[unmask]附加到染色体文件名，从而指示LASTZ忽略掩码信息，并将重复视为染色体的任何其他部分，以便准确评估reads映射的唯一性。因为我们知道这两种物种很接近，所以我们想降低灵敏度。使用-‑step=10，我们将只在每10bp中寻找种子。我们使用-‑seed=match12和-‑notransition来代替默认的种子模式，从而我们的种子将精确匹配12个碱基。我们使用-‑exact=20来代替默认的x-drop extension method，从而需要20个碱基的精确匹配才能作为HSP。由于我们正在比对短读长序列，因此我们指定-‑noytrim，以便在有间隔的扩展过程中，不会将比对结果的末端修剪回最高得分位置。
我们使用更加严格的参数--match=1,5替换默认的分数集，从而用于遥远物种的比对。在这个参数下，每个匹配到的碱基会得+1分，而错配碱基会-5分。我们还使用--ambiguous=n，以便对Ns（即序列中的缺失数据，均已字母N代替）进行适当的评分。我们只对几乎涉及整个read的比对感兴趣，并且由于该物种很近，因此我们不希望具有低相似度的比对。因此，我们使用‑‑filter=coverage:90和‑‑filter=identity:95进行过滤。
对于输出，我们只对reads比对的位置感兴趣，因此我们使用‑‑format=general选项，并指定我们想知道reads在染色体的什么位置（name1,start1,lenth1），reads的名称和方向（name2,strand2）。这将创建一个制表符（tab）分隔的输出文件，每个比对结果区块占一行，这是一种非常适合其他程序进行下游处理的格式。例如，要计算我们映射的不同reads的数量，我们可以运行此Unix shell命令：
cat hg18_21_vs_reads.dat | grep -v "#" | awk '{print $4}' | sort -u | wc
三、种子，HSP，有间隔的比对（gapped alignment），以及链接（chaining）
本例使用牛和人的α-珠蛋白区域演示了初级比对处理阶段。该数据包含在LASTZ发行版的test_data/aglobin.2bit目录中，由人DNA的70Kbp片段和牛DNA的66Kbp片段组成。我们将遵循此示例完成种子寻找（seeding）、HSP的获取（对种子的无间隔的扩展）、HSP的链接，以及局部比对的获取（对HSP的有间隔的扩展）。
下面的代码以及图片显示了在这些片段中间的一个小窗口（3K bp）上默认seeding的结果。种子很短，且接近匹配；在这种情况下，每个种子为19 bp，可能有多达8个错配。此窗口中有338个种子，但有许多种子的区域与线段无法区分。
lastz \
aglobin.2bit/human[34000..37000] \
aglobin.2bit/cow[35000..38000] \
--nogfextend --nochain --nogapped
下面的代码和图片显示了高得分片段对（high-scoring segment pairs），这是将种子进行无间隔的扩展的结果。共有11个HSP（图中只有10个是明显可见的，因为其中一个被1-bp移位到下一个对角线）。请注意，许多种子被丢弃，因为它们的扩展得分低或存在重叠。
lastz \
aglobin.2bit/human[34000..37000] \
aglobin.2bit/cow[35000..38000] \
--gfextend --nochain --nogapped
下面的代码和图片显示了将HSP进行有间隔的扩展，得到的局部比对区块。共有四个区块。
lastz \
aglobin.2bit/human[34000..37000] \
aglobin.2bit/cow[35000..38000] \
--gfextend --nochain --gapped
然后我们缩小并显示整个序列的结果；红色框表示前面图中所示的小区域。
下面的代码和图片依次显示了HSP、有间隔的序列比对区块（gapped alignment blocks）、将序列比对区块链接为一条同线性的直线（或两条线，如果两条线上都有匹配）。请注意，仅从HSP就可以看出序列是如何排列的。
# HSPs
lastz \
aglobin.2bit/human \
aglobin.2bit/cow \
--gfextend --nochain --nogapped
# gapped alignment blocks
lastz \
aglobin.2bit/human \
aglobin.2bit/cow \
--gfextend --nochain --gapped
# chaining
lastz \
aglobin.2bit/human \
aglobin.2bit/cow \
--gfextend --chain --gapped
四、将一条序列与自身进行比对
当序列与自身对齐时，完整结果将包含每个对齐块的镜像副本。处理两个副本在计算上是浪费的。LASTZ可以用四种不同的方式处理这种情况。
- 只需为目标和查询提供相同的序列。在这种情况下，LASTZ不知道它正在将序列与自身对齐，并对两个副本执行完整计算。
- 指定
-‑notrivial选项。这在两个副本上执行完整计算，但不报告沿主对角线的琐碎比对区块。 - 指定
-‑self选项代替查询序列。LASTZ将通过仅使用每个镜像对的一个区块进行计算以节省工作，尽管它仍然通过从第一个镜像对重构第二个副本来在输出中报告两个副本。它还自动调用-‑notrival，以省略主对角线上的琐碎比对区块。这给出了与前一种方法相同的输出，但运行速度更快。 - 指定
-‑self代替查询序列，并添加-‑nomirror选项。在这种情况下，LASTZ只报告每个镜像对的一个副本，并省略了一个区块。
