首页分享计算生物学与生物信息学漫谈

计算生物学与生物信息学漫谈

来源：花匠小妙招时间：2025-01-17 22:22

我们之前讲完了序列比对的算法：

计算生物学与生物信息学漫谈-5-mapping算法-CSDN博客

我们先学习比对得到的BAM文件

1.mapping的前言

2.SAM和BAM文件

QNAME字段

FLAG

POS字段

MAPQ字段

CIGAR字段

RNEXT字段

SEQ字段

QUAL字段

1.mapping的前言

对于读取映射，我们通常会得到由高通量仪器产生的数百万条读取序列，并希望确定这些读取序列在参考基因组序列中的来源。对于大多数测序应用来说，读取序列的比对或映射是最慢且计算成本最高的步骤。这是因为映射程序会尝试确定每条读取序列相对于参考基因组的最可能来源点。从真核生物RNA测序产生的读取序列映射需要比对工具付出额外的努力。在真核生物中，基因的编码区（外显子）被非编码区（内含子）分隔开。由于RNA测序仅针对转录组或基因转录本，因此用于映射RNA测序读取序列的比对工具必须考虑到外显子的不连续性以及剪接区域检测的挑战。

确实，在进行比对之前，我们需要下载所研究物种的参考基因组序列，然后对参考基因组进行索引，以便在搜索和比对过程中可以轻松找到读取序列映射的位置。

计算生物学与生物信息学漫谈-4-参考基因组与Mapping准备-CSDN博客

对于大多数比对工具而言，对参考基因组进行索引是在执行读取映射之前的第一个步骤。上面，我们讨论了最常用的索引方法，用于存储和组织参考基因组序列，以便可以轻松搜索以确定比对读取序列的位置（坐标）。

比对工具面临的另一个挑战是：测序仪产生的读取序列可能由于碱基识别错误而无法精确比对到参考基因组序列中的某个位置，或者由于所研究个体DNA序列中的突变（替换、缺失或插入）而自然无法精确比对。为了克服这一挑战，比对工具必须采用一种策略来执行带空位的比对，而不仅仅是精确比对。

基本步骤如下：

读取映射是大多数测序应用所必需的，包括基于参考的基因组组装、变异发现、基因表达、表观遗传学和宏基因组学。如图所示，读取映射/对齐的工作流程包括下载所研究物种的正确FASTA格式参考基因组文件，使用“samtools index”命令对参考基因组序列进行索引，使用比对器对参考基因组进行索引，最后使用该比对器对清理后的读取进行映射。

请记住，在映射之前，必须执行质量控制步骤，如之前所述。

计算生物学与生物信息学漫谈-3-FastQC及FASTX-toolkit详解-CSDN博客

因此，当我们进入读取映射步骤时，应该已经清理了读取并修复了QC报告中显示的大部分失败指标。即使我们无法修复所有错误，我们也应该意识到这在最终结果中的影响。我们还讨论了如何下载一个生物体的参考基因组以及如何对其进行索引。

读取映射（Read mapping）是指在参考基因组中找到包含在FASTQ文件中的读取序列（reads）所对应位置的过程。然后将读取映射信息存储在SAM/BAM文件格式中，这是一种用于存储序列比对信息的特殊文件格式。大多数比对工具都能够执行精确匹配和非精确匹配，这对于找到可能具有某些碱基调用错误或与参考基因组遗传变异的读取序列的位置至关重要。不同的比对工具实现了不同的算法，以在索引的参考基因组中执行这两种查找，这些参考基因组存储在如后缀树、后缀数组、哈希表和BWT等数据结构中。虽然精确查找相对直接，但非精确匹配使用序列相似性来找到读取序列最可能的来源位置。尽管有多种方法可以测量序列相似性，但大多数比对工具使用汉明距离或莱文斯坦距离来根据阈值评分读取序列与参考基因组部分之间的相似性。一些比对工具使用种子扩展策略（seed-and-extend），通过多个不匹配的碱基扩展一个种子（完全匹配的子字符串），以允许映射具有碱基调用错误或变异的读取序列。大多数比对工具在局部序列比对中使用SW算法实现种子扩展策略。种子是通过从参考基因组序列创建重叠的k-mer（长度为k的子字符串或单词）来生成的。一些比对工具，如Novoalign和SOAP，使用trie或哈希表数据结构对k-mer进行索引，以实现快速搜索。

2.SAM和BAM文件

几乎所有的reads mapping程序将读取映射到参考基因组的对齐信息存储在序列对齐和映射（SAM）文件或二进制对齐和映射（BAM）文件中，后者是SAM的二进制形式。

SAM文件是一种可读的纯文本文件，主要用于存储与参考序列对齐的生物序列，但也可能包含未映射的读取。它是一个由TAB分隔的文本文件，主要包括两大部分：

（i）标题部分和（ii）对齐部分。

SAM文件的标题部分是可选的，如果存在，则必须位于对齐部分之前。标题部分中的每一行都必须以“@”符号开头，后面跟着两个字母的标题记录类型代码。记录类型代码可能还有两个字母的子类型代码。下表列出了SAM标题部分的两个字母代码并进行了描述：

下图显示了一个示例标题部分：

请注意，SAM文件以“@HD VN:1.0 SO:coordinate”开始，这表明使用了SAM版本1.0的规范，并且文件中的对齐按坐标排序。我们还可以注意到有几行@SQ标题行，每一行都是用于对齐的参考序列。@SQ标题包括序列名称（SN），即染色体编号，和序列长度（LN）。标题部分的最后两行应该包括@PG，它描述了用于对齐的程序，以及@CO，它描述了命令行。

对齐部分在标题部分之后开始。每一行对齐都有11个必填字段来存储基本的对齐信息。对齐部分可能有数量不等的可选字段，这些字段用于提供额外的和对齐器特定的信息。

下图显示了一个SAM文件的部分对齐部分：

对齐部分的列被分开了，因为它们不适合页面。

下表列出了SAM对齐部分的11个必填字段并进行了描述：

这11个必填字段在SAM文件中始终存在。如果这些必填字段的任何信息不可用，则该字段的值将被替换为“0”（如果其数据类型为整数）或“*”（如果数据类型为字符串）。大多数字段名称都是自解释的。

QNAME字段

用于显示读取序列的名称，该名称从FASTQ文件中获得。

FLAG

是一个按位整数，描述比对情况（例如，成对、未比对和重复）。

这些描述以代码形式存储，如下表所示：

SAM文件中的FLAG字段可能具有表中列出的十进制值之一，或者是这些十进制值的任意和。例如，下图显示了读取“SRR062634.6862698”的FLAG为99，这意味着这个FLAG结合了4个条件：1+2+32+64=99。

这表示读取映射到参考基因组的位置如下描述：读取是成对的（1），比对器正确地映射了两个配对（2），下一个序列（SEQ）是反向链（32），以及是配对中的第一个读取（64）。而不是通过心算来计算这样的FLAG数字，我们可以使用“samtools flags”命令，如下所示：

samtools flags 99

samtools flags

POS字段

指定了对齐读段的第一个碱基在参考基因组中的最左侧映射位置。例如，序列“SRR062634.6862698”的对齐位置为10001。未映射读段的POS值为“0”。

MAPQ字段

指定了Phred评分单位中的映射质量。如果无法获得质量分数，则将其设置为255。

CIGAR字段

包含CIGAR字符串，这是一个由基本长度和相关操作组成的序列。它用于指示在对齐的读段中，发生了哪些操作，如匹配/不匹配、删除或插入。下表列出了可以与读段对齐关联的CIGAR操作。

例如，序列“SRR062634.6862698”具有CIGAR字符串“30M70S”，这意味着该序列的前30个碱基与参考序列（30M）中的碱基匹配，而70个碱基显示软剪切（70S），这是由于参考序列中的错误而导致的特殊不匹配。

RNEXT字段

包含模板中NEXT读段的主要对齐的参考序列名称。如果不存在，则将其设置为“=”或“”。PNEXT字段指定了模板中NEXT读段的主要对齐的位置。当信息不可用时，此字段将被设置为“0”。 TLEN字段指定了有符号观测模板长度。对于最左侧的读段，该字段将为正数；对于最右侧的读段，该字段将为负数；对于任何中间读段，该字段将未定义。

SEQ字段

包含读段序列。如果未存储序列，则将其设置为“”。否则，序列必须等于CIGAR字符串中操作的长度之和。“=”符号表示读段的碱基与参考碱基相同。已映射的读段对齐到正向基因组链（加号链）。因此，如果一个读段映射到参考序列的反向链（减号链），则SEQ字段将包含未映射加号链的互补链。

QUAL字段

包含ASCII字符（Phred+33）中的Phred质量分数。此字段中的字符串与FASTQ文件中的字符串相同。如果没有存储质量字符串，则可以将其设置为“”。否则，它必须等于SEQ中序列的长度。 SAM文件的对齐部分可能包含多个可选字段。每个可选字段由标准标签、数据类型和值定义，格式如下： TAG:TYPE:VALUE TAG是一个两字符字符串。SAM可选字段有几个预定义的标准标签。

完整的列表可以在https://samtools.github.io/hts-specs/SAMtags.pdf找到。用户可以添加新标签。

TYPE是一个定义字段数据类型的单个字符。它可以是“A”表示字符数据类型，“B”表示一般数组，“f”表示实数，“H”表示十六进制数组，“i”表示整数，以及“Z”表示字符串。 VALUE是由标签数据类型定义的字段的值。

在图2.16所示的SAM文件中，最后四列是通过四个预定义的标准标签识别的可选字段：“NH”、“HI”、“AS”和“NM”。其中，“NH”标签显示了当前记录中包含读段序列的报告对齐数量（命中次数）。“HI”标签显示查询命中的索引。“AS”标签显示由对齐器定义的对齐分数。“NM”标签显示编辑距离，即定义为将读段序列转换为参考序列对齐片段所需的最小单核苷酸编辑次数（替换、插入和删除）。有关SAM文件的更多详细信息，请阅读序列对齐/映射文件格式规范，该规范可在“https://samtools.github.io/hts-specs/”上获取。