上一期敛助，跟大家簡單介紹了下單細胞ATAC的背景知識點及其10x ATAC基礎數(shù)據(jù)的獲取方式萌业。接下來就帶大家從fragment.csv核偿、singlecell.csv晋修、peaks matrix等數(shù)據(jù)出發(fā)恒界，做單細胞ATAC的亞群分析喂柒。
與單細胞轉(zhuǎn)錄組類似踩蔚，單細胞ATAC的分析流程也主要包括細胞質(zhì)控棚放、peaks標準化及其降維分群、marker基因的鑒定等幾個步驟馅闽。常用的單細胞ATAC分析流程軟件包含 cell-ranger-atac飘蚯、Signac和ArchR等馍迄。
一、細胞質(zhì)控
單細胞ATAC的質(zhì)控點一般包含以下幾個方面：樣本重復（biological replicates）局骤，bulkATAC vs scATAC的相關性攀圈、fragment length distribution、per nucleus read-depth峦甩、transcription start site (TSS) enrichment赘来、雙細胞比例等。
1.1 低質(zhì)量細胞的過濾
前面提到的樣本相關性和fragments的長度分布主要是從整體水平上檢查我們的單個樣本數(shù)據(jù)的可靠性凯傲。
而要去掉不符合質(zhì)控的細胞犬辰，我們主要從fragments 數(shù)目和TSS enrichment score這兩點出發(fā)。
?? fragments 數(shù)目：一般指單個細胞（barcode）所屬的total fragments數(shù)目冰单。這個不同的軟件具體的定義不同忧风，比如cell-ranger-atac和Signac指peaks所屬區(qū)域的fragments 數(shù)目，其中singlecell.csv文件中peak_region_fragments列便是指fragments 數(shù)目球凰，而ArchR是指全基因組所有的fragments 數(shù)（這個跟該軟件的分析策略有關，后面會提到）腿宰。
?? TSS enrichment score：相當于計算每個細胞的信噪比（signal-to-background ratio）呕诉，ENCODE項目已經(jīng)定義了一個ATAC-seq目標評分，該評分基于TSS中心的片段與TSS側(cè)翼區(qū)域的片段的比例(見https://www.encodeproject.org/data-standards/terms/)吃度。較差的ATAC-seq實驗通常會有較低的TSS濃縮分數(shù)甩挫。Signac軟件可以用TSSEnrichment()函數(shù)為每個細胞計算TSS enrichment score，而ArchR包也是利用類似的原理createArrowFiles()函數(shù)在讀取基礎數(shù)據(jù)時就為每個細胞計算了該指標椿每。
備注：fragments 數(shù)目&TSS enrichment score的閾值不僅與所用軟件具體的計算公式有關（不同的軟件具體的參數(shù)可能不同）伊者，也與自己數(shù)據(jù)的實際情況有關。比如哺乳動物和植物的單細胞ATAC數(shù)據(jù)TSS enrichment score就不能用相同的指標cutoff來衡量间护，一般來說哺乳動物的TSS enrichment score值要整體偏高些亦渗。
1.2 雙細胞的去除
雙細胞預測幾乎是所有單細胞測序技術都得考慮的一個問題，從原理來說汁尺，我們每個barcode就是一個細胞法精，但是因為所有的實驗技術都不是100%完美的，因此往往會有一個barcode所包裹的油滴進來2個細胞痴突。
對于10x數(shù)據(jù)來說搂蜓，即使在使用標準試劑盒時，也可能有超過5%的細胞屬于雙細胞辽装，這對聚類產(chǎn)生了重大影響帮碰。特別是在發(fā)育/軌跡分析中十分受影響，因為doublets看起來像是兩種細胞類型的混合物拾积，這可能與中間細胞類型或細胞狀態(tài)混淆殉挽。
為了預測哪些“細胞”實際上是雙細胞的丰涉，ArchR會從我們真實的數(shù)據(jù)中隨機模擬產(chǎn)生混合的“雙細胞”數(shù)據(jù)，這些“雙細胞”數(shù)據(jù)與我們所有細胞一起做降維并UMAP可視化（"雙細胞"會投影到UMAP中此再，并識別它們鄰近的細胞）昔搂，在這個過程中，ArchR會計算每個細胞的Doublet Enrichment输拇，值越大摘符，表示該細胞是雙細胞的可能性越大。
二策吠、降維分群
與單細胞RNA(scRNA-seq)相比逛裤，scATAC-seq數(shù)據(jù)由于其高維度和稀疏性而更具計算分析挑戰(zhàn)性。主要體現(xiàn)在標準化和降維猴抹，這兩大步驟跟單細胞轉(zhuǎn)錄組分析所用的統(tǒng)計學原理完全不同带族，以下為歸納總結(jié)的具體內(nèi)容，如下表所示：
備注：TF-IDF & LSI都是自然語言常用的統(tǒng)計學方法蟀给。
2.1 peaks標準化
獲得peak matrix后蝙砌，跟基因類似，我們必須對其標準化跋理。因為單細胞ATAC測的是DNA序列择克，對于二倍體物種來說，同一個位置最多有2套DNA序列前普，這便是單細胞ATAC peak matrix稀疏性的最大根源（單細胞轉(zhuǎn)錄組因測的是RNA肚邢，高表達的基因往往有多個轉(zhuǎn)錄分子）。因此拭卿，從數(shù)據(jù)實際情況出發(fā)骡湖，單細胞ATAC采取的是log(TF-IDF)( Term frequency-inverse document frequency) 標準化，簡稱文檔頻率法峻厚。
- TF-IDF：是一種用于信息檢索與數(shù)據(jù)挖掘的常用加權技術响蕴。TF是詞頻(Term Frequency)，IDF是逆文本頻率指數(shù)(Inverse Document Frequency)惠桃。單個詞匯在一篇文章中出現(xiàn)的次數(shù)越多换途，越重要。但是在語料庫多次出現(xiàn)刽射，重要性越來越低军拟。IDF : 計算A term 出現(xiàn)稀少度。越稀少誓禁，越重要懈息。
2.2 peaks降維
所有高維數(shù)據(jù)的分析都是采取降維的方式從多維到低緯的策略，之后還可以再次降維成2個維度并可視化（比如TSNE和UMAP）摹恰。我們對peaks是采取LSI降維的方式辫继。
- LSI：潛在語義索引(Latent Semantic Indexing,以下簡稱LSI)怒见，有的文章也叫Latent Semantic Analysis（LSA）。其實是一個東西姑宽，后面我們統(tǒng)稱LSI遣耍，它是一種簡單實用的主題模型。LSI是基于奇異值分解（SVD）的方法來得到文本的主題的炮车。
2.3 細胞分群
與單細胞轉(zhuǎn)錄組類似舵变，降維后的單細胞ATAC數(shù)據(jù)也同樣可以采取graph-based clustering的分群方法。Graph-based圖聚類算法包括兩步：首先用降維（PCA或者LSI）的數(shù)據(jù)構建一個細胞間的k近鄰稀疏矩陣瘦穆，即將一個細胞與其歐式距離上最近的k個細胞聚為一類纪隙，然后在此基礎上用Louvain算法進行模塊優(yōu)化(Blondel, Guillaume, Lambiotte, & Lefebvre, 2008)，旨在找到圖中高度連接的模塊扛或。最后通過層次聚類將位于同一區(qū)域內(nèi)沒有差異表達基因(B-H adjusted p-value 低于0.05)的cluster進一步融合绵咱，重復該過程直到?jīng)]有clusters可以合并。
備注：Signac和ArchR都是直接調(diào)用Seurat包的FindClusters()函數(shù)用不同分辨率來分群的熙兔。
三悲伶、marker基因的鑒定
細胞分群后，我們需要知道每個cluster屬于什么細胞類型住涉，也就是細胞命名麸锉。我們知道，單細胞轉(zhuǎn)錄組主要是依據(jù)每個cluster的marker基因來判斷細胞類型的秆吵。那么對于單細胞ATAC，是不是也可以定義出每個cluster的特異高表達的基因集呢五慈？
答案是肯定的纳寂，一般來說，我們是通過基因body區(qū)域加上一定范圍內(nèi)的上下游區(qū)域的整體ATAC信號來計算每個細胞每個基因的genescore值泻拦。
3.1 Genescore的計算
1）Signac是通過GeneActivity()函數(shù)https://satijalab.org/signac/reference/geneactivit來實現(xiàn)的毙芜，默認參數(shù)是基因上游2kb到TES區(qū)域。
2）而ArchR是通過addGeneScoreMatrix()函數(shù)https://www.archrproject.com/reference/addGeneScoreMatrix.html來實現(xiàn)的（createArrowFiles函數(shù)也會用默認參數(shù)得到genescore matrix矩陣）争拐，注意其計算原理稍微復雜腋粥，ArchR考慮到遠端調(diào)控元件對基因活性的影響，因此默認的upstream和downstream范圍更廣架曹。
在ArchR作者的發(fā)表文章中隘冲，他們測試了50多個不同的基因評分模型，并確定了一類在各種測試條件下表現(xiàn)始終優(yōu)于其他模型的模型绑雄。這個模型類展辞，在ArchR中作為默認實現(xiàn)，有三個主要組件:
- a.整個基因體內(nèi)的可及性有助于基因得分万牺。
- b.一種指數(shù)加權函數(shù)罗珍，以一種距離依賴的方式來解釋假定的遠端調(diào)節(jié)元件的活動洽腺。
-
c.施加基因邊界，使不相關的調(diào)控元素對基因得分的貢獻最小化覆旱。
圖3：ArchR包Gene Scores計算原理
3.2 marker基因的可視化
marker 基因的ATAC信號（genescore值）同樣可以在umap展示蘸朋，也可以用小提琴圖（VlnPlot），點狀圖（DotPlot）展示扣唱。與單細胞轉(zhuǎn)錄組相比藕坯，單細胞ATAC還多了基因區(qū)域的track的可視化展示。
-
1）以下為Signac包里單細胞ATAC marker 基因的ATAC信號（genescore值）結(jié)果展示圖画舌，Signac包與Seurat包一樣堕担，都是satijalab實驗室團隊開發(fā)的，因此該包繼承了很多我們熟悉的Seurat包的方法函數(shù)曲聂。以下為例圖舉例霹购。
圖4：Signac包單細胞ATAC可視化展示例圖 -
2）ArchR包同樣也可以做很多可視化的圖。左邊是CD14基因genescore值umap可視化展示朋腋，右邊是track可視化圖齐疙。
圖5：ArchR包單細胞ATAC可視化展示例圖
單細胞ATAC的亞群分析介紹就到這里，下一篇會給大家介紹單細胞ATAC的高級分析內(nèi)容旭咽，比如motifdeviation贞奋、 擬時間分析、 單細胞RNA與單細胞ATAC的整合分析等穷绵。
本分享更多是從知識點和分析原理來講解和歸納總結(jié)轿塔，具體實現(xiàn)方法和流程腳本可以查看下面參考資料軟件的官方文檔，里面都寫得都很詳細清楚仲墨。
四勾缭、參考資料
- 1.https://support.10xgenomics.com/single-cell-atac/software/pipelines/latest/what-is-cell-ranger-atac # cell-ranger-atac
- 2.https://satijalab.org/signac/articles/pbmc_vignette.html #Signac官方教程
- 3.https://github.com/GreenleafLab/ArchR/ #ArchR github網(wǎng)站
- 4.https://www.archrproject.com/bookdown/creating-arrow-files.html #ArchR官方教程
- 5.Granja, J. M., et al. (2021). "ArchR is a scalable software package for integrative single-cell chromatin accessibility analysis." Nat Genet 53(3): 403-411.https://www.nature.com/articles/s41588-021-00790-6
- 6.Single-cell multiomic analysis identifies regulatory programs in mixed-phenotype acute leukemia. Nature Biotechnology (Granja JM, Klemm SK, McGinnis LM*, et al. 2019)https://www.nature.com/articles/s41587-019-0332-7 #人scATAC文章
- 7.A cis-regulatory atlas in maize at single-cell resolution. https://www.cell.com/cell/fulltext/S0092-8674(21)00493-1#articleInformation #玉米scATAC文章。5月7日發(fā)表在Cell上目养。
