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軟體包中非常常用的工具之一，能夠方便地統計不同標誌位的測序數據，獲取測序數據的重要質量指標，並進行數據分析和後續的數據處理。運用相應的命令並正確解讀輸出信息，可以在測序分析過程中發揮重要作用。