scSLAM-seq reveals core features of transcription dynamics in single cells
摘要
單細胞轉錄組測序（single-cell RNA-seq, scRNA-seq）強調了細胞間表達異質性在健康和疾病表型變異中的重要作用。然而，目前的scRNA-seq方法僅僅提供了基因表達的一個快照，很少傳達關於轉錄的真實時間動態和隨機特性的信息。scRNA-seq分析的另一個關鍵限制是每個細胞的RNA圖譜只能分析一次。原文作者開發了單細胞巰基連接的RNA烷基化代謝標記測序技術scSLAM-seq，它整合了代謝RNA標記、生化核苷轉化和scRNA-seq，通過區分每個單細胞數千個基因的新老RNA來直接記錄轉錄活性。通過使用scSLAM-seq來研究裂解性巨細胞病毒在單個小鼠成纖維細胞中的感染情況。從舊RNA推斷出的細胞周期狀態和感染劑量使得能夠基於新RNA進行劑量反應分析。因此，scSLAM-seq既可視化又解釋了單細胞水平轉錄活性的差異。此外，scSLAM-seq描述了宿主基因表達中的“開-關”控制和轉錄動力學，以及與啟動子內在特征(BP–TATA-box相互作用和DNA甲基化)相關的廣泛基因特異性差異。因此，基因特異性而非細胞特異性的特征解釋了單個細胞之間的轉錄組異質性和轉錄對擾動的反應。
<更多精彩，可關注微信公眾號：AIPuFuBio，和大型免費綜合生物信息學資源和工具平台AIPuFu：www.aipufu.com>
前言
SLAM-seq技術可將細胞短暫暴露於核苷類似物4-硫尿苷(4sU)。4sU在轉錄過程中被整合到新的RNA中，並在RNA測序前用碘乙酰胺轉化為胞嘧啶類似物。源自新RNA的測序讀段（reads）可以在總RNA讀段庫中根據特征性的鈾轉化為碳的轉化進行鑒定。通過應用SLAM-seq技術在單細胞水平上解決裂解性小鼠巨細胞病毒(CMV)感染的發作。在優化單細胞測序(scSLAM-seq)(圖1)后，原文作者在107個小鼠成纖維細胞上構建scSLAM-seq，並使用(大量)SLAM-seq並行分析匹配較大(1
× 10)細胞群體(n = 2)的全局轉錄變化。在對具有超過2，500個可檢測基因的細胞進行質量過濾后，剩余樣品(49個巨細胞病毒感染的，45個未感染的細胞)顯示了高質量scSLAM-seq文庫的所有特征，包括4%至6%之間的U-C轉化率。因此，在單細胞水平上，4sU的結合既有效又均勻。
結果和討論
由於4sU摻入率約為在50-200個核苷酸中有1個，因此所有來源於新RNA的SLAM序列讀段中，有高達50%可能不含鈾轉化成碳。為了克服這個問題，原文作者開發了“使用SLAM-seq技術對新轉錄的RNA和衰變率進行全局精確分析的方法（GRAND-SLAM）”，其為一種貝葉斯方法，以包括可信區間在內的完全量化的方式計算新RNA與總RNA的比率（new to total RNA，NTR）(圖1)。原文作者還開發了GRAND-SLAM 2.0用於並行分析數百個來源於單個細胞的SLAM-seq文庫。通過在成對端模式下分析長讀段(150個核苷酸)，進一步提高了定量的准確性，這使得4sU轉化能夠有效地與重疊序列中的測序誤差區分開。從而獲得了每個細胞數千個基因的精確測量值(90%可信區間< 0.2)，從而接單細胞轉錄組測序scRNA-seq的總靈敏度，並與多細胞SLAM-seq 的數據高相關(相關系數R > 0.73)。
高度可變細胞基因的無偏主成分分析（PCA）結果顯示，不能區分感染巨細胞病毒的細胞和未感染的細胞，無論是總RNA還是舊RNA，對新RNA只能有輕微的區分(圖2a)。因此，細胞間異質性超過了病毒誘導的變化，這種變化在感染后兩小時內很難在總RNA中檢測到，這是因為哺乳動物RNA的周轉緩慢。相比之下，NTR上的主成分分析可高精度從感染細胞中分離出未感染的細胞(圖2a)，並顯示出與病毒基因表達程度的明顯正相關(皮爾遜相關系數R = 0.59，P-value= 7.3 × 10)。
最近的研究表明，從單細胞轉錄組測序scRNA-seq數據中獲得的內含子讀段可以用來估計個體細胞中基因表達的時間導數，稱為“RNA速度”。這些表明了單個細胞的未來軌跡可投射到基因表達的低維空間中。然而，受感染的細胞不能通過無偏主成分分析從未受感染的細胞中分離出來，該主成分分析是根據相應的RNA速度計算的，或根據使用速度預測未來的表達譜計算的，或直接根據內含子/外顯子比率計算的。為了直接將scSLAM-seq與針對更大細胞群計算的RNA速度進行比較，原文作者在數百個未感染(n = 793)和巨細胞病毒感染(n = 353)的細胞上使用相同的實驗條件進行了10x Genomics Chromium droplet-based 的單細胞轉錄組測序。盡管成熟轉錄物(僅外顯子讀取)上的主成分分析不能分離未感染和感染的細胞，但使用內含子/外顯子比率可以進行區分。然而，在scSLAM-seq和10x數據中沒有觀察到有意義的RNA速率(圖2b)。原文作者使用scSLAM-seq獲得的新的和總的RNA水平來代替內含子和外顯子的read水平，並確定“NTR速度”。值得注意的是，這些進一步區分了感染細胞和未感染細胞(圖2c)。
為了直接比較NTRs和RNA的速度，原文作者問進一步分析了哪一個能最好地預測一個基因在大細胞群中是上調還是下調。雖然這在某種程度上是可能的，分別使用從scSLAM-seq或10x數據中的數十或數百個細胞計算的RNA速度(AUC值分別為0.68和0.74)，但是通過使用NTRs (AUC > 0.94)可獲得更好的結果(圖2d)。此外，NTRs相對於velocities (分別來自scSLAM-seq和10x數據，n =131和n = 295)，可以確定更多受調控的基因(n =583)。因此，NTR和RNA速率傳遞了不同但互補的信息，並可以合並到NTR速度中，以更可靠地預測細胞的未來狀態。
病毒學中的一個基本問題是為什么一個細胞的感染會導致廣泛的裂解性病毒復制而感染的第二個細胞沒有。在感染2小時后，大多數細胞的巨細胞病毒感染已經從“立即早期”(僅限於少數基因)發展到“早期”感染階段，其中大多數病毒基因已經轉錄。盡管所有細胞中的大多數病毒轉錄物都是新的，但原文作者也觀察到一些病毒基因的大量舊RNA。這代表病毒粒子相關的RNA，由傳入的病毒粒子傳遞到受感染的細胞。它與從病毒源中分離的病毒顆粒相關RNA有很好的相關性(相關系數R =0.48)，因此，提供了每個細胞接受的感染劑量的替代標記。相比之下，新的病毒RNA反映了感染的效率。感染劑量解釋了52%的感染效果差異。所以，幾乎不表達任何新病毒RNA的細胞中的舊病毒RNA的量比其它細胞中的低得多，這表明這些細胞對巨細胞病毒感染的許可性並沒有降低，而是接受了低得多的病毒劑量。基於細胞周期特征基因，可以分別從舊的和總的RNA推斷4sU代謝標記開始和結束時的細胞周期狀態，從而提供細胞周期軌跡(圖2e)。雖然裂解性感染始於所有細胞周期階段，但在G1期感染的細胞導致明顯更強的病毒基因表達和細胞周期破壞(P-value <0.05)。這使得新病毒基因表達的解釋差異增加到59%(圖2f)。因此，在成纖維細胞中單個細胞水平啟動裂解性病毒基因表達的效率可以通過感染劑量和細胞周期的相互作用得到很好的解釋。
為了評估巨細胞病毒感染對細胞基因表達的影響，原文作者使用單細胞兩相檢測程序(SC2P)從總的、新的和舊的RNA中鑒定了巨細胞病毒感染和未感染的單細胞之間差異表達的基因。大多數下調(87個基因中的61個)和上調(309個基因中的188個)基因(超過60%)只有通過專門考慮新的RNA才能被發現(圖3a)。基因表達的雙峰性是單細胞中一個很好描述的特征。雙峰性表達的基因在一個細胞亞群中檢測不到表達，但可在其他細胞群中檢測到表達。scSLAM-seq直接觀察細胞中給定基因的啟動子在研究的時間范圍內是否“開啟”。值得注意的是，原文作者發現大多數巨細胞病毒誘導的新RNA的變化更多地與開關動力學相一致，而不是與上調或下調相一致(圖3b，c)。
巨細胞病毒感染在感染的前兩個小時誘導強烈的I型干擾素和NF-κB反應。根據預測的轉錄因子靶標和基因本體條目（Gene Ontology terms）進行的基因富集分析表明，干擾素和NF-κB的激活高度依賴於病毒劑量。然而，雖然干擾素的激活只限於大約一半的感染細胞，但在巨細胞病毒感染的大多數細胞都發生了NF-κB激活 (圖3d)。病毒劑量依賴性激活所有細胞中的NF-κB與M45被膜蛋白介導的IKK激酶復合體處或上游的NF-κB激活一致。相比之下，干擾素反應需要檢測病原體相關的分子模式，受自分泌和旁分泌信號的影響，因此可能在單個細胞之間表現出更大的可變性。為了進行更集中的分析，原文作者根據以前發表的關於NF-κB誘導的數據和干擾素治療的bulk SLAM-seq數據，定義了巨細胞病毒感染特異的NF-κB和干擾素應答基因組。干擾素(圖3e)和NF-κB(圖3f)反應的大小在單個細胞之間明顯不同，但兩者都與病毒基因表達高度相關(斯皮爾曼相關系數ρ > 0.52，P-value < 3 × 10)。因此，大多數NF-κB-和干擾素誘導型基因表達出現在最強感染的細胞中，在S期感染的細胞中誘導最弱的反應。
單細胞水平上的轉錄活性不是一個連續的過程，而是由間歇的轉錄爆發組成，中間相隔幾分鍾到幾小時的轉錄不活躍，表明啟動子暫時不允許。原文作者推斷，scSLAM-seq應該可以檢測到基因的這種爆發因子。為了全面評估轉錄爆發動力學，原文作者將基因爆發評分(B-評分)定義為所有未感染細胞的正常參考時間分布的標准偏差，其中相應基因的RNA可以可靠地量化(正常參考時間的90%可信區間< 0.2；n = 5540)。從兩個獨立的生物重復中獲得的B-分數高度相關(相關系數R =0.74)。在某些情況下，接近0.5的極端B-分數對應於在每個細胞中僅檢測到單個或非常少的mRNA分子(新的或舊的)的基因。
當前的單細胞轉錄組測序（scRNA-seq）技術要么提供了獨特分子的識別標簽(UMIs)來估計捕獲的mRNA分子的數量，並且是鏈特異性的，但是僅克隆轉錄了末端(例如，10x Genomics Chromium單細胞轉錄組測序)，或者可測全長的mRNAs，但是沒有UMIs並且失去鏈特異性(例如，Smart-seq2)。原文作者基於Smart-seq技術開發的scSLAM-seq方法至少部分包含了所有三個特性。4sU的隨機合並轉化為mRNA分子提供了基於核苷酸轉化的獨特分子標識(nUMIs)，這使得能夠估計每個細胞和基因取樣的新的mRNA分子數量的下限，並且——通過外推(EnUMIS)——也是舊的mRNA的下限。基於這些對樣本多克隆抗體數量的保守估計，發現所有可檢測基因中超過30%(5540個中的1718個；矯正后P-value< 0.01，χ檢驗)的方差(B-得分)大於采樣時的預期(圖4a)。B分數與表達水平的相關性可以忽略不計。此外，NTR中觀察到的異質性不是由細胞周期依賴性差異引起的，而是與mRNA半衰期相關的。
無偏基因本體（Gene Ontoloy）富集分析顯示，高B分數的基因與蛋白質磷酸化和泛素化等功能類別相關。啟動子分析確定了六個顯著富集的低 (TATA box motif)或高 (富含CG和嘌呤motif) B-分數。Correctly placed TATA boxes富集的最明顯(P-value< 10-8),與TATA box的啟動子一致，在幾分鍾的時間尺度上頻繁的轉錄爆發。沒有觀察到與其他核心啟動子motif的關聯(圖4b)。富含CG的motif可以對應於顯示振盪激活模式的特定轉錄因子的結合位點，或者反映相應啟動子內富含CpG的區域。超過50%的哺乳動物轉錄起始於靠近CpG島的啟動子(CGI啟動子)，這些啟動子代表幾十到幾百個核苷酸的富含CpG的區域。CpG島的高甲基化是基因沉默的表觀遺傳控制機制。來自小鼠成纖維細胞的亞硫酸氫鹽測序數據揭示了甲基化的降鈣素基因相關肽啟動子與低B-評分的顯著相關性，而甲基化的非降鈣素基因相關肽啟動子傾向於顯示高B-評分(圖4c)。對於含有TATA box的啟動子也觀察到同樣的情況(圖4d)。通過重復對1718個具有顯著B值的基因和前50%最強表達基因的分析，證實了這些結果。雖然不能完全排除這里所用細胞系的轉錄活性和核型復雜性的等位基因差異導致了一些觀察到的效應，但這並不能解釋B-評分與啟動子內在特征的強烈相關性。原文作者提出了一個模型，其中基因啟動子中的DNA甲基化通過在數小時內的時間尺度上暫時抑制基因的啟動子而參與短暫的啟動子活性調控。最后，B分數與觀察到的單個細胞轉錄組的異質性高度相關(圖4e)。因此，基因啟動子特異性特征是細胞間異質性的主要原因。
使用4sU的代謝標記適用於所有主要模式生物，包括脊椎動物、昆蟲、植物和酵母。嘌呤類似物6-硫代鳥嘌呤(6sG)現在也能通過氧化-親核-芳香取代(TimeLapse-seq chemistry)實現從G到A的轉換。短而連續的4sU和6sG脈沖，結合巰基-(SH)介導的核苷轉化，可以實現單細胞轉錄活性的兩個獨立記錄。最后，scSLAM-seq與基於CRISPR的擾動相結合，將可極大地提高各種方法破譯分子機制的靈敏度，這對發育生物學、感染和癌症具有重要意義。
<更多精彩，可關注微信公眾號：AIPuFuBio和大型免費綜合生物信息學資源和工具平台AIPuFu：www.aipufu.com>