掃碼聯系客服
單細胞ATAC和RNA測序揭示了與前列腺癌復發有關的原始和持續細胞
大家好，今天給大家分享的文章是今年9月發表在Nature Communications(IF:14.919)上面的一篇文章。前列腺癌具有異質性，那些對系統治療有反應的患者，如果存在某些方法對這些患者進行分層，那將使患者受益。這篇文章就是采用ATAC-seq和RNA-seq的單細胞檢測方法，對恩雜魯胺(一種抗癌藥物)的早期治療反應和耐藥性模型進行研究。
方法
RNA測序及預處理
LNCaP和VCaP細胞系是從ATCC機構獲得。VCaP細胞的RNA測序是用Illumina HiSeq3000進行的，實驗中對3個重復樣品進行了測序，平均每個樣品獲得1.11億個雙端reads。在預處理過程中，基于Ensembl參考基因組GRCh38使用STAR對reads進行了比對。使用featureCounts和Gencode注釋對read counts進行了定量。
單細胞樣品制備和測序
對于scRNA-seq，使用CellRanger將測序reads處理為FASTQ格式和單細胞特征counts。同樣地，使用CellRanger ATAC將scATAC-seq的測序reads處理為FASTQ格式，并計算峰-條形碼counts。LNCaP-ENZ168、VCaP和VCaP-ENZ48的Drop-seq(液滴測序)樣本經過預處理、比對，并使用Drop-seq工具處理成細胞計數矩陣，預計每個樣本包含1000個細胞。該流程使用STAR和Picard工具，并使用了人類參考基因組GRCh38和Gencode注釋。
甲醛輔助分離調控元件(FAIRE)測序和分析
FAIRE分離的DNA片段用羅氏KAPA試劑盒進行文庫制備，在Illumina HiSeq2500上測序，產生50bp的單端reads，并使用bwa與hg19進行比對。使用Picard對重復進行了標記和重比對。使用MACS2對比對文件進行peak calling，DiffBind被用來探索峰的重疊情況，并推導出共識峰集。通過計算每個樣本條件下重復的reads平均數，對共有峰位點、MYC(髓細胞增生原癌基因)結合位點和AR(雄激素受體)結合位點進行讀長分布分析，并使用t檢驗比較樣本之間的分布中心值，以評估這些位點的染色質開放性差異。
單細胞RNA預處理和質量控制
CellRanger的輸出被用作Seurat的輸入，用于進一步分析scRNA-seq樣品。對于每個樣品，根據檢測到的基因數量、檢測到的分子總數和來自線粒體基因組的reads百分比篩選出質量差的細胞。為了解決數據中細胞周期異質性的影響，使用Seurat CellCycleScoring函數對每個細胞的S期或G2/M期相關基因的表達進行評分。使用sctransform對G2/M期和S期之間的評分差異進行評估。在單細胞RNA聚類中，研究者使用相互最近鄰方法fastMNN對4個LNCaP樣本進行整合，使用了2000個整合特征并考慮了批次效應。使用Seurat進行聚類和UMAP非線性降維。每個簇的特征基因和樣本間的差異表達基因是用Seurat的廣義線性模型MAST框架來識別的。如果Bonferroni校正的P值<0.01，簇中至少有10%的細胞表達該基因，且平均倍數變化至少為0.25，則該基因被認為是差異表達的。此外，研究者利用MSigDB中的hallmark基因集，根據其差異表達的基因來表征聚類和樣本。具體的步驟是使用GSVA包進行了基因集變異分析(GSVA)，以表征每個聚類的平均表達量。然后使用fgsea包對MSigDB hallmark基因集進行基因集富集分析。使用cytoTRACE預測了每個樣本中每個細胞的分化潛能。使用scVelo評估了scRNA-seq樣本中單細胞的RNA速度。
單細胞ATAC預處理和質量控制
CellRanger ATAC流程的輸出被用作Signac軟件包的輸入，用于進一步分析scATAC-seq樣品。對于每個樣品，篩選出劣質細胞。使用Signac和潛在語義索引(LSI)進行數據歸一化和降維。使用LSI嵌入的harmony方法對scATAC-seq樣本進行整合聚類。得到的結果被用作Signac軟件包的輸入，用于UMAP非線性降維和聚類，使用默認參數和SLM算法進行模塊優化。通過匯集每個樣品中所有優質細胞的reads，對scATAC-seq樣品中的染色質可及性變化進行分析。使用Signac中的TSSEnrichment函數生成轉錄起始位點的富集，CoveragePlot函數生成基于片段覆蓋率的染色質可及性軌跡。使用ReactomePA評估Reactome通路的過表達。用Signac的邏輯回歸模型來識別聚類中的差異可及區域，該模型根據每個基因來預測群組成員，并使用似然比檢驗來比較結果和零模型，以峰的總數作為潛在變量。如果Bonferroni校正的p值<0.05，至少有10%的細胞在該區域顯示出可及性，并且平均倍數變化至少為0.25，則該區域被認為是可及的。使用SignacClosestFeature函數用最接近的基因對差異可及區域進行注釋。使用R包ggradar對樣品間簇的差異可及區域進行可視化。
scATAC-seq的轉錄因子基序富集和轉錄輸出
使用R包TFBSTools、BSgenome.Hsapiens.UCSC.hg38和從RJASPAR2018數據包中檢索的JASPAR數據庫的位置-頻率矩陣，在樣品條件之間和每個樣品簇之間的差異可及染色質區域中使用Signac進行轉錄因子基序富集。考慮到染色質區域的序列特征(如GC頻率)，因此使用了超幾何檢驗來檢驗顯著性基序富集。通過評估富集的轉錄因子的靶基因和scRNA-seq簇中差異表達基因之間的重疊，將scATAC-seq中的染色質狀態與scRNA-seq中的轉錄輸出相關聯。
基于標簽轉移整合scRNA-seq和scATAC-seq數據集
LNCaP，LNCaP-ENZ48，RES-A和RES-B的每個樣本都有scRNA-seq和scATAC-seq數據。這些數據類型使用Signac和Seurat中實現的聚類-標簽轉移程序進行整合。每個scRNA-seq樣本被單獨聚類，其聚類標簽被投射到匹配的，單獨聚類的scATAC-seq樣本上，每個樣本的聚類分辨率使用clustree進行評估和確定。另外，RNA-seq的表達水平從scATAC-seq數據中推算出來。通過測試不同的預測分數閾值，發現在所有樣本中，使用0.3的閾值，大約有50%以上的細胞在數據類型之間進行了標簽轉移。
RNA-seq和臨床數據分析
使用GSVA包對每個基因特征或基因集進行富集評估并在樣本中評分。在單細胞水平評估基因集表達的情況下，Seurat中的AddModuleScore函數被用來生成每個細胞的平均表達分數。生存分析使用survival包進行，Kaplan-Meier曲線使用survminer包繪制。對于單一特征的生存分析，使用中位GSVA評分將患者分為低表達和高表達的特征組。對于多個特征的生存分析，利用每個特征的GSVA富集分數，使用歐氏距離和層次聚類法對樣本進行聚類。聚類結果隨后被用來確定生存分析中的兩組樣本。
原發性前列腺癌組織的空間轉錄組學分析
從一名患者身上獲得了兩片冷凍的PC組織，在Illumina NovaSeqPE150測序儀上完成測序。測序數據首先用10xGenomics的SpaceRanger進行處理，以獲得兩個部分的表達矩陣。然后用Seurat進行下游處理和聚類，用sctransform對數據進行歸一化處理，以考慮測序深度差異。Seurat的AddModuleScore函數被用來對基因特征的spot進行評分。研究者比較了管家基因集和scRNA-seq特征集的基因表達分數分布，以確定第90個百分點作為cutoff，在這個分界點上把spot視為具有高表達的scRNA-seq特征。
結果
染色質重編程是恩雜魯胺耐藥的基礎
為了研究PC中AR信號抑制和耐藥動力學的分子后果，利用LNCaP親本細胞系、LNCaP衍生的ENZ耐藥細胞系RES-A和RES-B，這些細胞系通過長期暴露于AR靶向藥物而產生，以及其他獨立生成的LNCaP和VCaP衍生模型(圖1a)。為了確定染色質結構對ENZ抗性的影響，對4個樣本進行了scATAC測序(圖1a)。與親本細胞相比，ENZ抗性細胞中轉錄起始位點(TSS)的ATAC-seq信號下降(圖1b)。這種趨勢也被發現在管家基因、雄性激素反應通路的基因和參與MYC信號傳導的基因中，表明這種模式不限于一個特定的基因集。此外，與RES-B和LNCaP相比，RES-A細胞的獨特開放位點比例更高(圖1c)。
研究者通過對親代LNCaP和RES-A細胞進行FAIRE測序，證實了ENZ抗性細胞中染色質開放和重編程的程度。雖然開放染色質位點的總數沒有明顯差異，但與親代相比，在雄性激素存在和雄性激素匱乏的條件下，ENZ抗性樣本有更高比例的獨特開放位點。接下來，研究者基于scATAC-seq樣本來生成具有不同染色質可及性特征的細胞亞群的聚類可視化(圖1d)，由此確定了各樣品中獨特或共享的簇(圖1e)。獨特的簇是專門針對RES-A和RES-B，共有的簇以相似的比例存在于樣本中，被命名為持續簇(圖1e)。將每個簇與其他所有的簇進行比較，根據差異可及的染色質區域(DAR)來確定其獨特的染色質特征。就細胞數量而言，最普遍的基于染色質的scATAC-seq簇是持續的(圖1e)，這表明在ENZ抗性的發展過程中，74%的細胞共享一個整體類似的染色質可及性譜。然后，評估了親代LNCaP、LNCaP-ENZ48以及RES-A和RES-B之間簇染色質DAR的變化。在對恩雜魯胺的短效反應期間，在幾個簇的MYC和TP53周圍觀察到DAR。
已有研究表明，在無雄性激素的情況下長時間培養的PC細胞系往往顯示出類似神經內分泌的表型。結果發現，RES-A和RES-B細胞中NEPC衍生的特征表達升高，以及原始簇中NEPC下調的基因表達升高。此外，在相同細胞系的RNA測序中，NEPC特征的基因集變異分析顯示，只有RES-A細胞中NEPC上調的基因表達量較高。這些數據顯示，在出現對AR靶向藥物的耐藥性過程中，染色質發生了廣泛的重編程。
恩雜魯胺的耐藥性重構了染色質中TF結合DNA基序的可用性
染色質的可及性通過暴露TF DNA結合基序的軌跡來確定細胞的轉錄輸出。根據reads分布分析，觀察到ENZ抗性細胞中MYC結合位點的開放染色質明顯增加(圖2a)，去勢條件下AR結合位點的開放染色質減少(圖2b)。這些發現表明，ENZ抗性的染色質失調與AR和MYC染色質結合的重構有關，這與之前報道的一致。
為了解決染色質重編程如何在單細胞水平上影響TF DNA基序的暴露，研究者對每個樣品中scATAC-seq細胞簇的標記DAR進行了TF基序富集分析(圖2c)。在RES-A和RES-B中，簇3和簇5富集了相同TF基序的子集，簇5在所有樣品中顯示出一致的富集趨勢(圖2c)。簇3以FOXA1和JUND為特征，而簇5以CTCF、ETS和MYC為特征，盡管它們擁有不同的DAR，但ENZ誘導的簇6和簇7在RES-A或RES-B中并沒有顯示出TF基序的富集。值得一提的是，MYC和ESR1分別是RES-A和RES-B中所有簇中最常見的(圖2d)。這些分析表明，ENZ耐藥性與TF DNA基序軌跡的重構有關。
恩雜魯胺耐藥性的轉錄模式由不同的染色質重編程引起
為了研究與單細胞水平的染色質結構重組有關的轉錄模式，對LNCaP親本、RES-A和-B模型進行了scRNA-seq。4個LNCaP樣本的綜合聚類(圖3a)顯示了7個持續、3個ENZ誘導和3個初始細胞簇，由標記的差異表達基因集構成(圖3b)。為了證實這些細胞亞群與ENZ抗性的其他獨立模型有關，使用標簽轉移方法探索獨立scRNA-seq數據集中的匹配細胞群。轉移scRNA-seq簇標簽證實了LNCaP親代和RES-C中原始簇的存在，在RES-C和LNCaP-ENZ168中證實了ENZ誘導簇的存在。由于大多數scATAC-seq簇的細胞比例與scRNA-seq的簇3細胞相對應(圖3e)，這表明該簇的細胞可能代表了ENZ抗性的基因組。對VCaP細胞的分析證實了VCaP親代細胞以及VCaP-ENZ48的原始和ENZ誘導細胞的普遍性(圖3c)。接下來，將scRNA-seq簇與其匹配的scATAC-seq簇連接起來，再次利用標簽轉移的方法，在相同的樣本條件下確定匹配的scRNA-seq和scATAC-seq細胞狀態。結果發現，一個染色質狀態可以對應多個轉錄狀態。通過探索scATAC-seq數據中的scRNA-seq簇，結果可以在scRNA中找到所有scATAC簇的匹配細胞狀態，其中scATAC簇中的細胞通常對應于多個scRNA簇(圖3d-e)。總之，這些數據表明，ENZ抗性細胞的轉錄構型，特別是在治療期間的持續細胞，是由染色質結構和TF介導的轉錄重編程驅動過程產生的，這些過程影響許多的細胞命運調節因子。
具有干性特征的前列腺癌細胞亞群先于恩雜魯胺耐藥
細胞周期階段可能是scRNA-seq數據整合聚類的一個重要決定因素。在Seurat中使用細胞周期評分法，持續簇8、9和11的S期和G2/M期相關基因得分很高(圖4a)，表明這些簇中的細胞循環和增殖更為活躍。然而，簇11中的細胞不僅具有細胞周期基因的特征，還參與染色質重塑等過程(圖4b)。簇5和11顯示了一個基因集的高表達，該基因集表征了具有干細胞樣、雄激素不敏感和細胞周期驅動特征的細胞亞群(圖4c)，將這一基因特征命名為Persist。另外，在親代LNCaP的ENZ治療前，原始簇10對PC相關基因特征的評分很高，并將其命名為PROSGenesis(圖4d)。最后，將PROSGenesis和Persist基因在scRNA-seq樣品中的表達進行了可視化，以證實這些亞群細胞在其他模型中的存在(圖4e)。
使用CytoTRACE根據表達基因的數量估計分化狀態，在大多數樣品中，簇11中的細胞顯示出高發展潛力(圖4f)，表明其他細胞亞群可以從該簇的細胞中衍生出來。RNA速度分析預測簇10是恩雜魯胺誘導的簇前體(圖4g)。RES-A和RES-B的簇差速分析顯示，許多PC相關基因如ATAD2下調，以及UBE2T等基因的上調。另外，ATAD2和UBE2T在持續簇8、9和11中上調，這表明在ENZ誘導的簇中發生了額外的轉錄重編程。這些分析得出了ENZ耐藥之前兩個不同的PC細胞亞群：一個與Persist匹配的持續細胞簇(簇11)和一個與PROSGenesis匹配的原始簇(簇10)。
前列腺癌RNA測序中基因特征的表征
根據基因組變異分析，大多數持續簇和簇10顯示E2F靶點、G2M檢查點和MYC靶基因的富集，表明在RNA-seq數據中可以檢索到細胞亞群的信號。簇內的差異表達和基因集富集分析進一步顯示，氧化磷酸化在LNCaP-ENZ48中被上調，這一過程在RES-A中被選擇性的維持，但在RES-B細胞中未被維持。此外，在ENZ耐藥性的產生過程中，受活化mTORC1信號調節的基因在大多數簇中一致上調，這與之前的報道一致。大部分情況下，ENZ誘導的DEG選擇性地出現在RES-B細胞中。同樣，當用DHT誘導時，持續簇只在RES-A和RES-B中與持續特征相關。另一方面，PROSGenesis特征僅在RES-B中升高(圖5)。總的來說，持續簇，原始簇，PROSGenesis和Persist基因特征具有從RNA-seq來識別前列腺癌侵襲性、再生特征的潛力。
轉錄信號富集分析識別了前列腺癌患者的治療持續細胞和預后基因特征
研究者獲得了已有研究ENZ治療CRPC患者的臨床數據，對單個基因特征的PFS分析顯示，Persist特征基因得分較高的患者與PFS較短有關(圖5c，d)，而PFS較長的患者在PROSGenesis和原始簇特征上得分高。這些數據表明，對ENZ耐藥的單細胞分析中，與持續細胞(簇11)相關的Persist特征是一個有效的分類器，可能對患者進行分層，以確定對AR靶向治療的反應(圖5f)。此外，研究者基于最近發表的PC相關的scRNA-seq數據集，分析顯示與成纖維細胞相比，LNCaP模型衍生細胞簇在管腔和基底/中間細胞中得分更高。此外，與基底/中間細胞和原始細胞相比，管腔細胞具有更高的與原始scRNA-seq簇相關的基因表達。然后，對來自Persist和PROSGenesis特征的基因表達細胞及其相關簇和一些對照特征進行評分(圖6b)。在Persist特征評分較高的細胞中，約48%為管腔細胞(圖6c)。在PROSGenesis特征中得分較高的細胞大多為基底細胞/中間細胞(圖6d)。每個腫瘤平均有8%的細胞在Persist和PROSGenesis特征上得分較高(圖6e)。
為了探索這些細胞的存在和它們的相對組織病理學位置，用空間轉錄組學評估了原發性PC切片內的基因表達(前列腺A和B)，利用聚類分析重建了基質和上皮成分的基因表達信號，并在基質組織(ST)、良性上皮(BE)和腺癌(PC-AC)的5個簇中注釋了組織結構(圖6f)。PROSGenesis和Persist特征，以及伴隨模型衍生的簇10特征，與來自管家基因特征的分數相比，在切片內顯示出高表達分數(圖6g)。在兩個切片中，高信號的spot分布在所有5個簇中(圖6h)。然而在前列腺A中，與ST相比，PC-AC簇中 Persist 特征得分高的spot更為普遍(圖6h)。這些數據表明，在PC患者未經治療的原發性前列腺組織中，存在著治療的持續細胞，具有很高的轉移潛力。
最后，研究者驗證了是否能利用從這些細胞中提取的特征基因預測原發性PC患者的復發，由此獲得了原發性腫瘤TCGA PRAD(圖7a)和早發性PC ICGC RNA-seq數據。使用所有特征進行聚類，TCGA PRAD隊列分離出54%的GS(Gleason分級)-7和15%的GS-8+患者，他們不會從額外的治療中受益，因為他們的預后相對較好(圖7b)，在ICGC隊列中也觀察到類似的趨勢。在TCGA隊列中，ENZ誘導的簇(圖7c)、PROSGenesis(圖7d)、Persist(圖7e)和持續簇(圖7f)基因特征是簇分離的最重要因素，而NEPC下調基因是ICGC隊列的主要決定因素。這些腫瘤中反映AR活性的特征在TCGA隊列中與較長的進展時間相關(圖7g，h)，表明AR驅動的腫瘤對AR信號的抑制有更好的反應。在早發性PC隊列中，與TCGA PRAD隊列相比，持續簇和PROSGenesis特征對GS-7患者進行了顯著的分層，表明這些特征能夠進一步確定基于GS的患者風險分層，避免過度治療。高PROSGenesis評分與良好的預后以及來自原始簇10的基因集相關(圖7i)。在TCGA隊列中，13個簇衍生特征中有8個與PFS相關(圖7i)，表明這些特征在PC患者風險分層中的重要作用。
參考文獻: Single-cell ATAC and RNA sequencing reveal preexisting and persistent cells associated with prostate cancer relapse
