samtoolsflagstat：从多个方面进行详细阐述

一、原理解析

samtoolsflagstat命令是SAMtools软件包中用于检查SAM/BAM格式文件的工具。SAM（Sequence Alignment/Mapping）格式是常用于存储测序数据的文件格式之一，而BAM格式是SAM格式的压缩版本。samtoolsflagstat命令能够快速统计测序数据的关键指标，通常用于质量控制、评估测序数据的成果以及后续分析的数据处理。

具体来说，samtoolsflagstat可以针对SAM/BAM文件中的一系列flag标志进行统计，其中包括：

• Total：总的序列数
• Secondary：次要的序列数，即没有通过重名比对或没有通过质量阈值筛选的序列数
• Supplementary：补充的序列数，即通过重名比对但不是最优比对或没有通过质量阈值筛选的序列数
• Duplicates：复制的序列数，即PCR扩增引入的序列重复次数
• Mapped：比对上的序列数
• Paired in sequencing：成对测序的序列数
• Read1/Read2：成对测序中Read1或Read2比对上的序列数
• Properly paired：正常成对测序的序列数，即成对测序中两个序列一前一后比对上了相同染色体的不同位置，并且方向正确
• With itself and mate mapped：两个序列都比对上了相同染色体的不同位置，并且方向一致的序列数
• Singletons：只有一个序列比对上了相同染色体的不同位置的序列数
• With mate mapped to a different chr：成对测序中两个序列比对上了不同染色体的序列数
• With mate mapped to a different chr (mapQ>=5)：比对质量大于等于5的成对测序中两个序列比对上了不同染色体的序列数

通过以上标志位的统计，我们可以翻译得到很多关键的测序参数，如比对率、重复率、成对比对率以及单端比对率等重要参数，这些参数能够帮助我们更好地判断测序数据的质量，进行下一步进一步的数据分析并会更精细。

二、使用方法

为了使读者更好地了解如何使用samtoolsflagstat，下面将使用三个例子进行阐述。

1. 统计BAM文件的总序列数和比对上的序列数

$ samtools flagstat example.bam

运行以上命令后，samtoolsflagstat会统计BAM文件example.bam中各个标志位的信息，并输出统计结果，如下所示：

1000 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
1000 + 0 mapped (100.00% : N/A)
1000 + 0 paired in sequencing
500 + 0 read1
500 + 0 read2
1000 + 0 properly paired (100.00% : N/A)
1000 + 0 with itself and mate mapped
0 + 0 singletons (0.00% : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

可以看到，在example.bam文件中，共有1000个序列，全部通过了质量控制，没有次要序列或补充序列，也没有PCR复制序列。其中1000个序列比对上了基因组序列，且比对成功率达到了100%。成对测序的数据中，read1和read2数量都是500，所有序列都正常成对比对。

2. 比较两个BAM文件的比对率和重复率

$ samtools flagstat -@ 4 -r example.bam control.bam > compare.txt

通过以上命令可以比较example.bam和control.bam两个BAM文件的比对率和重复率，-@选项指定使用4个线程处理任务，-r选项表示只考虑比对上的序列。结果会输出到compare.txt文件中。比对结果如下所示：

1000 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
1000 + 0 mapped (100.00% : N/A)
1000 + 0 paired in sequencing
500 + 0 read1
500 + 0 read2
1000 + 0 properly paired (100.00% : N/A)
1000 + 0 with itself and mate mapped
0 + 0 singletons (0.00% : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

900 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
900 + 0 mapped (100.00% : N/A)
900 + 0 paired in sequencing
450 + 0 read1
450 + 0 read2
900 + 0 properly paired (100.00% : N/A)
900 + 0 with itself and mate mapped
0 + 0 singletons (0.00% : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

可以看出，两个文件中序列数都是900和1000，全部通过质量控制，并且没有次要序列或补充序列。example.bam文件中有1000个比对上了基因组序列的序列，比对成功率达到了100%。control.bam文件中有900个序列比对成功了。这样，我们可以发现两个文件中比对率的差异，以及bame文件的重复率。

3. 计算比对质量分数大于20的序列比对率

$ samtools view -q 20 example.bam | samtools flagstat -

以上命令首先使用samtoolsview命令筛选出比对质量分数大于20的序列，并将结果传递给samtoolsflagstat命令进行统计。结果如下所示：

610 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
610 + 0 mapped (100.00% : N/A)
610 + 0 paired in sequencing
305 + 0 read1
305 + 0 read2
610 + 0 properly paired (100.00% : N/A)
610 + 0 with itself and mate mapped
0 + 0 singletons (0.00% : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

可以看到，example.bam文件中有1000个序列，其中610个比对质量分数大于20。在这610个序列中，全部通过了质量控制，没有次要序列或补充序列，也没有PCR复制序列。其中610个序列比对成功，比对成功率达到了100%。成对测序的数据中，read1和read2数量都是305，所有序列都正常成对比对。

三、常见问题

1. samtoolsflagstat能统计哪些参数？

samtoolsflagstat统计能够统计的参数包括：总序列数、次要序列数、补充序列数、PCR复制序列数、比对上的序列数、成对测序序列数、正常成对测序序列数、两个序列比对上了相同染色体的不同位置的序列数以及成对测序中两个序列比对上了不同染色体的序列数。

2. samtoolsflagstat的输出结果是什么？

samtoolsflagstat的输出结果包括各项标志，如Total、Secondary、Supplementary、Duplicates、Mapped、Paired in sequencing等，以及每项标志的统计数量，并添加统计比例。

3. samtoolsflagstat是否可用于统计单端测序数据？

可以，但是值得一提的是，实验室往往在关注单端测序数据比对情况时，不会统计properly paired、with itself and mate mapped、with mate mapped to a different chr连接器，因为单端数据自带信息有限。

结语

samtoolsflagstat是SAMtools软件包中非常常用的工具之一，能够方便地统计不同标志位的测序数据，获取测序数据的重要质量指标，并进行数据分析和后续的数据处理。运用相应的命令并正确解读输出信息，可以在测序分析过程中发挥重要作用。