測序數(shù)據(jù)量對于NGS數(shù)據(jù)分析是非常重要的,，測序數(shù)據(jù)量過低,，不能有效覆蓋基因組完整信息，測序數(shù)據(jù)量過高,，則會造成冗余,，不夠經(jīng)濟(jì)。為了驗證當(dāng)前測序量能否滿足需求,，或者說加大測序量是否能夠進(jìn)一步挖掘的更大量的信息,，通常需要進(jìn)行飽和度分析。

在轉(zhuǎn)錄組和擴(kuò)增子測序中,，飽和度分析應(yīng)用的已經(jīng)非常成熟,，對于WES/WGS等基因組測序而言,，對測序深度也已經(jīng)基本形成了統(tǒng)一共識。對于ATAC而言,，測序飽和度分析要怎么做呢,？本片文章就來解密一下。

ATACseqQC作為一個專門針對ATAC文庫進(jìn)行QC的R包,，在其文章中,，提供了測序飽和度分析的思路和方法，文章鏈接如下

https://bmcgenomics./articles/10.1186/s12864-018-4559-3

如上圖所示,，ATACseqQC的飽和度分析思路如下,，對effectice fragment進(jìn)行隨機(jī)抽樣，比如10%,，20%到100%,， 對于每個梯度分別進(jìn)行peak calling, 統(tǒng)計其peak個數(shù)，然后以effectice fragment數(shù)目為橫坐標(biāo),，peak個數(shù)為縱坐標(biāo)繪制散點圖,，并進(jìn)行擬合,，如果擬合曲線的末端斜率上升平緩,，表明繼續(xù)加大測序數(shù)據(jù)量，新識別到的peak個數(shù)也不會新增太多,，就可以認(rèn)為當(dāng)前的測序數(shù)據(jù)量是比較飽和的,，不比再加測數(shù)據(jù)量了。如果曲線末端仍處于一個斜率急劇上升的階段,，說明測序數(shù)據(jù)量不夠,，還需加大測序數(shù)據(jù)量。

理解了分析思路,，再來看一下具體的操作過程,。第一步是對effectice fragment進(jìn)行抽樣。所謂effectice fragment,，即過濾掉PCR重復(fù),，去除了線粒體之后的有效序列。在實際分析中,，質(zhì)控后的clean reads首先比對到參考基因組上,，生成原始的bam文件，也稱為raw bam,。然后對原始的bam文件進(jìn)行過濾,，去除MAPQ較低的alignment, 去除PCR重復(fù)，去除比對到線粒體的序列,，經(jīng)過重重篩選,，得到了一個filter之后的bam文件,，這個filter bam文件中保存的就是effectice fragment。

對effectice fragment進(jìn)行抽樣,，本質(zhì)就是操作這個bam文件,。通過linux下的shuf命令對sam文件隨機(jī)抽取對應(yīng)的行，從而實現(xiàn)對effectice fragment抽樣的目的,，基本思路如下

1. 提取sam文件的header部分作為一個單獨的文件,；

2. 使用linux內(nèi)置的shuf命令對sam文件的正文部分隨機(jī)抽取一定比例的行數(shù)，追加到header文件中,，生成一個新的sam文件

在文章的附件中提供了一個bash腳本,，用于對bam文件進(jìn)行抽樣，用法如下

第一個參數(shù)為filter之后的bam文件,，第二個參數(shù)為輸出文件的前綴,，其他參數(shù)為抽取的百分比，0到1之間,，運行成功后,，會生成一系列文件，以0.1為了,，生成的文件如下

利用排序好的bam文件可以直接進(jìn)行peak calling, 然后統(tǒng)計每個抽樣百分比對應(yīng)的effectice fragment數(shù)目和peak個數(shù),，就可以繪制測序飽和度圖了。在ATACseqQC中,，提供了saturationPlot函數(shù),，可以直接讀取一系列peak文件，然后繪制測序飽和度圖,，具體的用法可以查看該函數(shù)的幫助文檔,。

ATACseqQC給我們提供了一個很好的思路，來進(jìn)行ATAC文庫的測序飽和度分析,，類似的,，這個做法也可以推廣到chip_seq, m6A_seq等IP類型的實驗中去。