–https://doi.org/10.1038/s41588-024-01682-1

Tissue-specific enhancer–gene maps from multimodal single-cell data identify causal disease alleles

研究團隊和單位

Alkes L. Price–Broad Institute of MIT and Harvard

Soumya Raychaudhuri–Harvard Medical School

研究簡介

SCENT（單細胞增強子靶基因圖譜）是一種用于分析10x Genomics單細胞多組學ATAC + RNA 多模態數據的模型，揭示增強子染色質可及性與基因表達之間的關聯。該模型原理是泊松回歸和非參數自舉法，繪制出增強子-基因對，并預測與表達相關的增強子可能具有重要功能。

增強子-基因圖譜的構建：

SCENT應用于9個多模態數據集，包含超過120,000個單細胞或細胞核，構建了23個細胞類型特異性的增強子-基因圖譜。

這些圖譜在表達定量基因座（eQTL）和1,143種疾病及性狀的全基因組關聯研究（GWAS）中高度富集了致病變異。

功能驗證：

1.SCENT模型認定的增強子在統計上精細定位的因果變異中表現出富集。

2.SCENT模型結果通過了CRISPR-Flow FISH和H3K27ac HiChIP實驗驗證。

3.SCENT預測的增強子-基因連接具有較高的功能性和進化保守性。

4.模型成功識別了常見疾病和罕見疾病的潛在致病基因，并將體細胞突變熱點與目標基因聯系起來。在類風濕性關節炎、1型糖尿病和炎癥性腸病等自身免疫性疾病中，SCENT揭示了多個潛在的因果變異和基因。

5.與其他方法的比較：SCENT在識別因果變異方面優于現有的單細胞多模態方法（如ArchR和Signac）和基于bulk組織的方法（如EpiMap和ABC模型）。

模型方法

1.數據處理

單細胞轉錄組數據：標準質控流程后，基因表達 count 矩陣使用NormalizeData函數進行歸一化，使用ScaleData函數進行縮放，并使用Harmony進行批次校正。

單細胞ATAC數據：標準質控流程后，ATAC峰矩陣被二值化，如果計數>0則為1，否則為0。

什么peak對應什么基因呢？

研究是在利用基因組位置，把和基因本體位于同一染色體上，且距離基因本體500kb以內的peak劃分為該基因的peak。

2.SCENT方法

泊松回歸

SCENT方法原理是泊松回歸，核心目標是找到染色質可及性（ATAC-seq數據）和基因表達（RNA-seq數據）之間的關系。具體來說，想探究某個染色質區域（peak峰）的可及性是否會影響某個基因的表達。公式的主要內容如下：

${\beta }_0$ 是偏置項；${\beta }_{peak}$ 、${\beta }_{mito}$ 、${\beta }_{nUMI}$ 、${\beta }_{batch}$ 都是對應的系數；

$E_i$ 是原本測序的基因表達值；

${\lambda }_i$ 是模型預測出來的基因表達值；

$X_{peak}$ 是該基因對應的多個peak值（0/1）；

$X_{mito}$ 是細胞種線粒體基因表達的比例；

$X_{nUMI}$ 是細胞種唯一分子識別符的數量；

$X_{batch}$ 是實驗批次效應。

模型訓練的目標是找到最佳的參數（${\beta }_0$、${\beta }_{peak}$ 、${\beta }_{mito}$ 、${\beta }_{nUMI}$ 、${\beta }_{batch}$），使得模型的預測值${\lambda }_i$盡可能接近實際觀察值$E_i$。

模型擬合：

使用最大似然估計來擬合泊松回歸模型。

評估顯著性：

為了評估 βpeak 的顯著性（即染色質峰對基因表達的影響是否真實），SCENT使用了自舉法（Bootstrapping）。自舉法的核心思想是通過重采樣來模擬數據的分布。

具體來說：

從原始數據中隨機抽取細胞（有放回），生成一個新的樣本；

在新樣本上重新擬合模型，得到一個新的 ${\beta }_{peak’}$

重復這個過程很多次（比如100到50,000次），得到 ${\beta }_{peak’}$ 的分布

通過比較 ${\beta }_{peak’}$ 的分布與零假設（${\beta }_{peak’}$=0）來判斷 ${\beta }_{peak}$ 是否顯著。

多重檢驗校正：

由于SCENT同時測試了大量的峰-基因對，可能會產生假陽性結果。因此使用了FDR校正來控制假陽性的比例。

研究結果

1.SCENT模型概覽

通過統計方法（如泊松回歸）來測試每個峰的可及性是否與目標基因的表達水平顯著相關。如果某個峰的可及性與基因表達顯著相關，那么這個峰可能是一個功能性的調控元件，能夠影響該基因的表達。

作者在文章中認定這些peak是增強子，所以文章中以SCENT增強子去稱呼相關性系數高的peak。

SCENT精確地識別了調控區域染色質可及性與單細胞中單個基因表達之間的顯著關聯。

然而，對于高表達和分散的基因，泊松回歸可能不是最優選擇：

通過置換細胞條形碼，可以打亂ATAC-seq和RNA-seq數據之間的對應關系，從而破壞染色質可及性與基因表達之間的真實關聯，等于制造出空的數據集。

在空數據集中，由于ATAC-seq和RNA-seq數據之間的關聯已經被破壞，理論上不應該檢測到任何顯著的關聯。然而泊松回歸模型在空數據集中仍然檢測到了大量顯著關聯，說明模型的I型錯誤率過高，即假陽性率失控。

為了解決I型錯誤失控的問題，作者在后續的分析中引入了雙尾非參數自舉法（nonparametric bootstrapping），通過重采樣細胞來估計統計量的分布，從而更準確地評估顯著性。

2.發現新的細胞類型特異性SCENT增強子-基因關聯

數據包含9個單細胞多組學數據集，涵蓋了13種細胞類型（免疫相關、造血、神經元和垂體）。其中作者從11名類風濕性關節炎患者和1名骨關節炎患者的滑膜組織中自測了一個炎癥組織數據集（30893個細胞），并獲取了8個公共數據集，包含129672個細胞。

質量控制（QC）后，分析了16,621個基因和1,193,842個順式開放染色質峰（共4,753,521個峰-基因對，每個基因中位數為28個峰）。

SCENT模型用于每個細胞類型（ncells>500），最后一共構建了23個細胞類型特異性增強子-基因圖譜，共包含87,648個峰-基因關聯（錯誤發現率（FDR）<10%）。基因的關聯峰數量各不相同（從0到97個，平均4.13個）。

ATAC-seq片段數量越多，可能意味著更高的染色質開放區域檢測靈敏度，從而有助于更準確地識別與基因表達相關的染色質區域；RNA分子數量越多，有助于更準確地識別與染色質開放區域相關的基因：

為了評估SCENT選出的peak是否具有功能，作者開展以下幾點分析：

• （1）它們是否與傳統順式調控注釋共定位；
• （2）它們對表達的影響是否在更接近的峰-基因對中更強；
• （3）它們是否具有較高的序列保守性；
• （4）峰-基因關聯是否有實驗支持。

2.1 測試SCENT選出的peak峰與ENCODE cCRE的重疊mapping情況

平均98.0%的SCENT峰與cCRE重疊，而隨機大小匹配的順式區域和非SCENT峰的重疊率分別為23.3%和89.0%

還使用18種狀態的chromHMM結果對免疫細胞類型中的SCENT峰進行了注釋；97.4%的SCENT峰在41個免疫相關樣本中與啟動子或增強子注釋重疊。

2.2 檢查增強子-基因關聯的強度

作者基于更強的關聯會更接近目標基因的轉錄起始位點（TSS）的假設。隨著SCENT峰接近TSS，回歸系數βpeak（峰可及性對基因表達的影響大小）變得更大且更正向。

2.3 評估了SCENT峰是否具有更高的phastCons分數

更高的phastCons分數，反映了跨物種的序列保守性；進化保守的調控區域在功能上活躍且富含復雜性狀的遺傳性。

結果如下圖，外顯子（Exonic）區域比所有非編碼順式區域進化保守性更高。

SCENT調控區域相對于非編碼順式區域也具有保守性。相比之下，所有順式ATAC峰與非編碼順式區域的ΔphastCons分數較低。

2.4 測試SCENT峰的目標基因是否富集于實驗驗證的增強子-基因關聯

作者使用了CRISPR-Flow FISH結果，其中包括278個陽性和5,470個陰性增強子-基因關聯。與實驗相關的組織中的SCENT峰顯著富集于陽性關聯。

其次，使用了naive T細胞、Th17 T細胞和調節性T細胞中分辨率高達1 kb的H3K27ac HiChIP數據。在我們的6個T細胞數據集中，SCENT峰-基因關聯在H3K27ac HiChIP增強子-基因環中富集1.6倍。

功能喪失不耐受概率（pLI）和功能喪失觀察/預期上限分數（LOEUF）分析

作者假設，具有最多SCENT峰的基因可能是功能上最受限且對功能喪失突變最不耐受的基因。

研究基于人類群體中缺乏有害變異的概率評估了突變約束性，包括功能喪失不耐受概率（pLI）和功能喪失觀察/預期上限分數（LOEUF）。每個基因的SCENT峰數量與基因的平均約束分數顯著相關（pLI的β=0.37，P=4.9×10?90，分數越高表示約束性越強；LOEUF的β=?0.35， P=?0.35×10?106，分數越低表示約束性越強）。

研究結果與之前的報告一致，即具有大量調控區域的基因在批量表觀基因組數據中富集于功能喪失不耐受基因。

3.SCENT峰中eQTL推定因果變異的富集

為了檢測SCENT峰是否在統計學上精細定位了表達數量性狀位點（eQTL）的推定因果變異，研究使用了GTEx中49種組織的eQTL數據，并將后驗包含概率（PIP）>0.2的變異定義為推定因果變異。

結果顯示順式區域內由ATAC-seq定義的所有可及區域在精細定位變異中均顯示出2.7倍的適度富集。SCENT峰在所有數據集中平均顯示出9.6倍的精細定位變異富集。

由于許多SCENT峰靠近TSS區域，作者考慮了這種富集是否由TSS鄰近性驅動。結果顯示TSS距離是因果eQTL變異的最重要特征之一：

此外，SCENT方法能夠以細胞類型特異性的方式檢測順式調控區域，并通過整合多個數據集構建了四個主要細胞類型的增強子-基因圖譜。而且，細胞類型特異性的SCENT增強子在GTEx相關樣本中的eQTL變異中富集程度最高：

這些結果展示了SCENT以細胞類型特異性的方式優先考慮因果eQTL變異的能力。

4.SCENT增強子中GWAS因果變異的富集

SCENT可用于通過將其應用于疾病相關組織的多模態數據來構建疾病特異性增強子-基因圖譜。

研究檢測了SCENT峰是否可用于優先考慮疾病因果變異。從兩個大規模生物庫（包括FinnGen中的1,046種疾病性狀和UK Biobank中的35種二元性狀及59種數量性狀）中獲取了PIP >0.2的候選因果變異。

結果可見SCENT增強子在FinnGen和UK Biobank中的GWAS因果變異中顯著富集：

研究從9個數據集和23種細胞類型中構建了SCENT，僅使用了28個樣本，遠少于構建EpiMap的833個樣本和組織以及構建ABC模型的131個樣本和細胞系。

盡管數據集較小，SCENT峰在給定數量的峰-基因鏈接中顯示出比ABC模型和EpiMap更高的GWAS變異富集，而且更嚴格的PIP閾值增加了富集：

研究假設SCENT識別的推定因果變異能夠調控染色質可及性。如果是這樣，SCENT增強子與caQTL的交集可能會進一步富集GWAS因果變異。

為了驗證這一點，研究使用了具有基因型數據的單細胞ATAC-seq樣本，并通過利用等位基因特異性染色質可及性或結合等位基因特異性和個體間差異進行caQTL定位。

分析結果觀察到SCENT與caQTL交集區域中的富集高于僅使用SCENT的區域。隨著caQTL峰閾值的嚴格性增加，富集度最高可達333倍。

例如，15q22.33處的一個哮喘GWAS位點包含一個位于SMAD3基因內含子內的SCENT增強子，該增強子攜帶一個推定因果變異rs172936320。該SCENT增強子具有顯著的caQTL效應，來自等位基因特異性效應和染色質可及性的個體間差異。

替代等位基因T降低了染色質可及性，并被報告破壞了保守的AP-1結合位點。等位基因T還降低了SMAD3的表達。SCENT識別的目標基因SMAD3參與TGF-β信號傳導，該信號在哮喘中重塑氣道。

5. 通過SCENT定義GWAS位點的機制

5.1 研究利用SCENT方法定義疾病的因果機制

研究分析了來自FinnGen、UK Biobank以及類風濕性關節炎、炎癥性腸病和1型糖尿病的GWAS隊列的精細定位變異。

SCENT將4,124個推定因果變異（PIP >0.1）與1,143個性狀中的潛在目標基因聯系起來。大多數目標基因靠近因果變異，其中20%是距離因果變異最近的基因：

但30.6%的SCENT關聯基因距離因果變異超過300 kb。

5.2 自身免疫位點的分析

由于SCENT軌跡主要來源于免疫細胞類型，研究首先關注自身免疫位點。

6q15位點：

T細胞特異性SCENT增強子內的單個精細定位變異rs72928038（PIP >0.3）。

相關疾病和目標基因：類風濕性關節炎、1型糖尿病、特應性皮炎和甲狀腺功能減退，距離該變異最近的基因是BACH2。

實驗驗證：在T細胞中將保護性等位基因編輯為風險等位基因的實驗證實了該變異對BACH2表達的影響。刪除rs72928038的小鼠初始CD8 T細胞更容易分化為效應T細胞。

4p15.2位點：

類風濕性關節炎和1型糖尿病種，候選變異21個，每個變異的PIP值較低（<0.14）。SCENT模型優先考慮的變異是來自炎癥滑膜組織的關節炎組織數據集中T細胞和成纖維細胞中的單個變異rs35944082，目標基因是 RBPJ。

RBPJ轉錄因子對NOTCH信號傳導至關重要，NOTCH信號傳導與類風濕性關節炎組織炎癥相關。Rbpj敲低導致T細胞分化異常并破壞調節性T細胞表型。

5.3 垂體數據的分析

11p14.1位點，與多種婦科性狀（子宮內膜異位癥、月經過多、卵巢囊腫和絕經年齡）相關，其目標基因是FSHB，在垂體中特異性表達，促進卵巢卵泡生成。

罕見的FSHB變異會導致常染色體隱性遺傳的低促性腺激素性性腺功能減退癥。來自其他組織的多模態數據和bulk組織方法未能優先考慮這一變異，因為它們遺漏了最相關的疾病組織——垂體。

6. SCENT增強子中的罕見疾病變異和體細胞突變

6.1 罕見非編碼變異的分析

目前，僅在約30-40%的孟德爾疾病患者中能夠識別出因果突變，許多病例中的變異被注釋為意義未明的變異（VUS），尤其是非編碼變異。

研究檢查了ClinVar報告的臨床非良性非編碼變異（總計400,300個變異）與SCENT增強子的重疊情況，ClinVar變異富集中，SCENT增強子中平均包含的ClinVar變異數量是相同基因組長度ATAC區域的2.0倍：

此外，ClinVar變異的密度分別是ENCODE cCREs和所有非編碼區域的3.2倍和12倍。SCENT為33,618個非編碼ClinVar變異定義了3,724個目標基因。

研究發現了40個與LDLR基因相關的非編碼變異，這些變異導致家族性高膽固醇血癥1；3個與IL10RA相關的非編碼變異，導致常染色體隱性早發性炎癥性腸病28（如下圖所示）；以及一個與ATM基因相關的內含子變異rs1591491477，導致遺傳性癌癥易感綜合征。

6.2 體細胞突變的分析

在白血病的17個非編碼突變熱點中，來自血液相關細胞類型的SCENT增強子包含了12個熱點。SCENT增強子-基因鏈接將這些熱點與已知的驅動基因（例如，BACH2、BCL6、BCR、CXCR4（下圖所示）和IRF8）聯系起來。

對比ABC模型：在某些情況下，SCENT為這些突變熱點提名的目標基因與原研究中使用的ABC模型不同。SCENT將白血病中chr14:105568663-106851785的一個體細胞突變熱點與免疫球蛋白重鏈相關基因（如IGHA1）聯系起來，這可能比ABC模型提名的ADAM6更具生物學相關性。

總結

1.SCENT方法優勢

1. I型錯誤控制：SCENT模型展示了良好的I型錯誤控制能力，優于常用的統計模型。

2. 因果變異富集：SCENT映射的增強子在eQTL和GWAS中顯示出高富集的推定因果變異，并優于之前的方法。

3. 樣本效率：盡管使用的樣本數量顯著較少，SCENT定義的增強子在因果變異富集方面與基于大塊組織的方法相當甚至更高。

4. 單細胞水平建模：SCENT在單細胞水平上進行建模，避免了通過將細胞聚合成單個樣本來模糊關聯。

2.SCENT的局限性

1. 增強子數量較少：與其他資源相比，SCENT增強子-基因圖譜中的增強子數量相對較少（圖2a）。

2. 擴展潛力：細胞數量與顯著SCENT峰-基因鏈接數量之間的線性關系（補充圖8）表明，將SCENT應用于更大的數據集將擴展當前的增強子-基因圖譜。

3. bulk組織方法的限制：基于大塊組織的增強子-基因圖譜可能更難以擴展，因為需要更多的樣本。

4. 因果機制的局限性：SCENT專注于基因順式調控機制以精細定位疾病因果等位基因。可能存在其他因果機制，例如通過反式調控效應、剪接效應或轉錄后效應起作用的等位基因。為了證明SCENT優先考慮的等位基因的因果關系，需要使用基因編輯技術進行實驗驗證。

5. 計算強度：由于泊松回歸和自舉法的使用，SCENT的計算強度高于之前的方法（例如，SCENT需要1.5×10^7 CPU秒，Signac需要2.5×10^5 CPU秒，ArchR需要2.2×10^2 CPU秒）。

6. 優化：在SCENT中實現了多線程和并行化選項，可以線性加快計算速度，但需要額外的計算資源。對于非常大的數據集，算法改進（如下采樣或聚合細胞）可能是有用的。

3. SCENT的實用性

1. 疾病組織相關圖譜：SCENT的真正實用性在于它能夠構建與疾病組織相關的增強子-基因圖譜。

2. 數據獲取：多模態單細胞數據可以從廣泛的人類原代組織中獲取，包括那些難以解離的組織，因為這些數據可以從核材料中獲取而無需組織解離。

3. 應用示例：FSHB位點：理解婦科性狀中的FSHB位點特別需要垂體圖譜；RBPJ位點：類風濕性關節炎中的RBPJ位點特別需要滑膜組織圖譜。