测序数据处理 —— 比对数据

测序数据处理 —— 比对数据

介绍

一般测序的原始数据需要先对其进行质控过滤，比如去接头，删除低质量的数据等。获得了比较干净的数据之后，就需要将序列比对到参考基因组，获取每条序列的位置信息，比对结果一般会存储为 SAM 或 BAM 格式的文件。

SAM 和 BAM 是两种相同格式的文件，只是一种是文本形式一种是二进制形式，二进制形式所占的存储空间更小，个人更喜欢用这种格式。

SAM（Sequence Alignment Map） 文件是以 TAB 作为分隔符的文本文件，分为可选的头部信息以及比对部分，如果包含头部信息，则必须放在比对结果之前。

注意 SAM 文件与 BAM 文件的坐标系是不一样的，也就是说 SAM 中比对到的参考基因组位置是从 1 开始计数的，通常使用闭区间，即同时包括起始位置和结束位置。

例如，区间 [1, 3] 表示的是序列的前三个元素。而 BAM 文件中的比对位置是从 0 开始计算的，使用左闭右开区间，即包括起始位置但不包括结束位置。例如，区间 [0, 3)表示的是序列的前三个元素。

常见的几种文件的坐标系统：

1. 1-base 的坐标系统：SAM、VCF、GFF 和 Wiggle
2. 0-base 的坐标系统：BAM、BCFv2、BED 和 PSL

BED 文件可能是最容易搞错的，不同人可能会使用不同的坐标系统，需要注意。

头部信息

每行头部信息都是以 @ 符号起始并跟上两个标识分类信息的字符，包括 HD|SQ|RG|PG|CO，除了 @CO 行之外，每种分类信息中都以 TAG:VALUE 的形式定义字段信息，不同分类中含有不同的字段，每个字段之间用 TAB 符分隔。其正则匹配形式为

/^@(HD|SQ|RG|PG)(\t[A-Za-z][A-Za-z0-9]:[ -~]+)+$/ or /^@CO\t.*/

下面列举一些常见的字段

HD

表示文件信息，如果有必须放在首行且只能有一行，字段包括

• VN：格式的版本，必须字段。

可选值的形式为 /^[0-9]+\.[0-9]+$/

• SO：比对结果的排序方式。

可选值为 unknown（默认值）、unsorted、queryname、 coordinate，如果值为 coordinate，排序的主键是 RNAME，其顺序由标题中的 @SQ 行顺序定义，次键是 POS 字段，如果主键和次键相同，则顺序是任意的。其中 RNAME 和 POS 为比对信息中的字段。

• GO：比对数据的分组方式。

可选值为 none (默认值)、 query ( 根据 QNAME) 和 reference（根据 RNAME/POS）

• SS：添加更多的排序字段。

形式为 (coordinate|queryname|unsorted)(:[A-Za-z0-9_-]+)+

SQ

表示参考序列信息，允许多行，其顺序代表了比对文件的排序顺序

• SN：参考序列名称（染色体或转录本名称），每行的名称必须唯一，必须字段
• LN：参考序列长度，必须字段
• AS：组装的基因组标识符
• M5：序列的 MD5 校验码
• SP：所属物种
• TP：分子拓扑结构，可为 linear（默认） 和 circular
• UR：参考序列的链接，以 http: 或 ftp: 开头

RG

表示reads 分组信息，允许多行，不需要排序

• ID：分组标识符，每行的标识符必须唯一，必须字段
• BC：样本的 barcode 序列
• DT：测序时间日期
• LB：文库名称
• PI：预测的中位插入片段长度
• PL：测序平台，ILLUMINA、ONT 和 DNBSEQ 等
• SM：样本名称

PG

表示使用到的软件信息，允许多行，比如 bowtie2、samtools 等

• ID：软件标识符，必须字段
• PN：软件名称
• CL：运行的命令行
• VN：软件版本

CO

表示注释信息，允许多行

简单示例

例如，下面的就是一个 bam 文件的头部信息

@HD     VN:1.6  SO:coordinate@SQ     
@SQ     SN:chr1 LN:249250621
@SQ     SN:chr2 LN:243199373
@SQ     SN:chr3 LN:198022430
@SQ     SN:chr4 LN:191154276
@SQ     SN:chr5 LN:180915260 
...
@SQ     SN:chrY LN:59373566
@RG ID:1    BC:GGACTCCT-TCTACTCT    PL:Illumina SM:test LB:Lib
@PG     ID:bowtie2      PN:bowtie2      VN:2.4.5 CL:"bowtie2-align-s --wrapper basic-0 --threads 20 -x hg19 --end-to-end --very-sensitive --no-mixed --no-discordant --phred33 -I 10 -X 700 -1 GSM4734764.clean.R1.fq.gz -2 GSM4734764.clean.R2.fq.gz"
@PG     ID:samtools     PN:samtools     PP:bowtie2      VN:1.14 CL:samtools sort -o GSM4734764.sorted.bam
@PG     ID:MarkDuplicates       VN:Version:2.27.5       CL:MarkDuplicates --INPUT GSM4734764.sorted.bam --OUTPUT GSM4734764.rmdup.bam --METRICS_FILE GSM4734764.metrics.txt --REMOVE_DUPLICATES true --TMP_DIR /tmp/tmpp68rehfh --MAX_SEQUENCES_FOR_DISK_READ_ENDS_MAP 50000 --MAX_FILE_HANDLES_FOR_READ_ENDS_MAP 8000 --SORTING_COLLECTION_SIZE_RATIO 0.25 --TAG_DUPLICATE_SET_MEMBERS false --REMOVE_SEQUENCING_DUPLICATES false --TAGGING_POLICY DontTag --CLEAR_DT true --DUPLEX_UMI false --FLOW_MODE false --FLOW_QUALITY_SUM_STRATEGY false --USE_END_IN_UNPAIRED_READS false --USE_UNPAIRED_CLIPPED_END false --UNPAIRED_END_UNCERTAINTY 0 --FLOW_SKIP_FIRST_N_FLOWS 0 --FLOW_Q_IS_KNOWN_END false --FLOW_EFFECTIVE_QUALITY_THRESHOLD 15 --ADD_PG_TAG_TO_READS true --ASSUME_SORTED false --DUPLICATE_SCORING_STRATEGY SUM_OF_BASE_QUALITIES --PROGRAM_RECORD_ID MarkDuplicates --PROGRAM_GROUP_NAME MarkDuplicates --READ_NAME_REGEX <optimized capture of last three ':' separated fields as numeric values> --OPTICAL_DUPLICATE_PIXEL_DISTANCE 100 --MAX_OPTICAL_DUPLICATE_SET_SIZE 300000 --VERBOSITY INFO --QUIET false --VALIDATION_STRINGENCY STRICT --COMPRESSION_LEVEL 5 --MAX_RECORDS_IN_RAM 500000 --CREATE_INDEX false --CREATE_MD5_FILE false --GA4GH_CLIENT_SECRETS client_secrets.json --help false --version false --showHidden false --USE_JDK_DEFLATER false --USE_JDK_INFLATER false     PN:MarkDuplicates

头部信息中比较常用的是序列名称（SN）和长度（LN）

比对信息

该部分的格式不同于头部，每一个比对结果一般为一行，每行包括 11 个必须字段以及可变数量的可选字段，使用 TAB 作为分隔符。如果必须字段中有不可用或不存在的信息，则该字段的值将使用占位符，由字段类型决定是 0 或者 *。

必须字段

11 个必须字段为

更具体来说

QNAME

表示查询模板名称：即片段化后的 DNA 或 RNA 片段，如果是双端测序，read1 和 read2 来自同一个模板，QNAME 一样

FLAG

序列标识符，由一个十二位的二进制数表示，每个位标识不同的含义

例如，一条 read 的 FLAG 值为 96 = 32 + 64，即该 read 为 read1，且 read2 比对到负链。

注意：

对于双端测序的 fastq 文件，由 R1 和 R2 进行标识，两端的测序结果分别保存在两个文件中，其 read 分别称为 read1 和 read2。可以理解为一个 DNA 或 RNA 片段，分别从两端开始向中间逐碱基测序，读取序列信息，会产生两条 read，是一个配对（pair）

MAPQ

比对质量，比对位置错误的概率取负对数（−10log10(P)），然后将其四舍五入到最近的整数，255 表示比对质量不存在

CIGAR

由操作符和数字组成的字符串，用于标识该 read 比对到参考基因组上的情况，如 100M75S 表示该 read 的前 100 个碱基匹配到当前位置，后面 75 个碱基匹配到基因组的其他位置，表示该处可能存在基因组结构变异

注意：

H 和 S 只能出现在两端；N 只在 RNA 比对到 DNA 上用到，表示内含子，其他情况并不会使用；M/I/S/=/X 的长度和就是该条序列的长度（query 列），M/D/N/=/X 的长度表示比对到的参考基因组的长度（reference 列）。

RNEXT

模版中另一条 read 比对上的参考序列的名字/编号

如果在头部信息中设置了 @SQ，则 RNEXT 的值必须在 SQ-SN 字段中出现（除非是 * 或者 =），其中 = 表示两条序列比对到同一条染色体

PNEXT

模版中另一条 read 比对上的参考序列上的起始位置。在 SAM 文件中，其为比对上的最左端从 1 开始的坐标

一般与 RNEXT 一起使用，用于描述配对的测序片段在基因组上的相对位置。在单端测序中，这两个字段通常都是星号（*）或等号（=），表示没有下一个比对。

TLEN

推测出的模版（上机测序的序列）片段的长度。

该值可正可负，正值表示第二条read 在参考基因组上的起始位置在第一条 read 的后面，需要根据两条配对 reads 的相对比对位置来推测

SEQ

测到片段的序列。如果不存储序列，该字段可以是 *，否则序列的长度必须等于 CIGAR 中 M/I/S/=/X 操作符对应长度的总和

QUAL

序列中每个碱基的测序质量，与序列一一对应

可选字段

SAM/BAM 文件预定义了非常多的可选字段，同时也支持自定义字段。这些字段的基本形式为：TAG:TYPE:VALUE。

其中 TAG 为字段名称（两个字符），TYPE 为字段值的类型，VALUE 则是具体的值。类型值可以是：字符（A）、数组（B）、实数（f）、十六进制数组（H）、整数（i）或字符串（Z）。
常用的几个字段有

SA 字段值一般在做 SV 分析中会用到，主要用于识别嵌合比对的 reads，也就是说序列中存在一个断点，断点两端的序列能够有效比对到基因组中不同的位置（可以是同染色体或不同染色体），这种序列指示了样本中可能存在结构变异。

其数据格式为：(rname ,pos ,strand ,CIGAR ,mapQ ,NM;)+ 可以有多条比对质量相同的记录，其中 pos 为 1 起始的位置。

自定义字段通常为 X?、Y?、Z? 或任何小写字符组成的字段（? 可以为任何 26 个字母），可以自由定义字段名和字段类型。

示例

下面列出了三条比对记录

SRR12491667.18894604    163     chr1    10074   30      65M     =       10074   65      AACCCTAACCCTAACCCTAACCCTAACCCTAACCCAACCCTAACCCTAACCCTAACCCTAACCCT       AAFFFJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJFJJJJJJ       MD:Z:65 PG:Z:MarkDuplicates     XG:i:0  NM:i:0  XM:i:0  XN:i:0  XO:i:0  AS:i:0  XS:i:-8 YS:i:0  YT:Z:CP
SRR12491667.18894604    83      chr1    10074   30      65M     =       10074   -65     AACCCTAACCCTAACCCTAACCCTAACCCTAACCCAACCCTAACCCTAACCCTAACCCTAACCCT       JAJJJJAFJJJJFJJJJJJJJJJFJAJJJJJJJJJJJJJJJJFJJJJJJJJJJAFJJJJJFFFAA       MD:Z:65 PG:Z:MarkDuplicates     XG:i:0  NM:i:0  XM:i:0  XN:i:0  XO:i:0  AS:i:0  YS:i:0  YT:Z:CP
SRR12491667.27435019    99      chr1    10114   30      73M     =       10114   73      TAACCCTAACCCTAACCCTAACCCTAACCCTAACCCCTAACCCTAACCCTAACCCTAACCCTAACCTAACCCT       AAFFFJJJJJJJJJJJJJJJJFJJFJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJ       MD:Z:73 PG:Z:MarkDuplicates     XG:i:0  NM:i:0  XM:i:0  XN:i:0  XO:i:0  AS:i:0  YS:i:0  YT:Z:CP

参考链接：

https://samtools.github.io/hts-specs/SAMv1.pdf

https://samtools.github.io/hts-specs/SAMtags.pdf