這篇文章是 Landscape of Infiltrating T Cells in Liver Cancer Revealed by Single-Cell Sequencing
第一個單細胞尺度的肝癌浸潤性T細胞圖譜

年前一作的鄭春紅博士在VIPaper有個在線講座,，我前一晚啃paper時不幸直接睡到開講,，聽得云里霧里~

In Brief
對肝癌患者的T細胞群體作了單細胞測序分析，看到了浸潤性淋巴細胞的獨特亞型和克隆擴增,。

Highlights
· 單細胞RNA-seq揭示肝癌T細胞的獨特功能組份
· 浸潤調(diào)節(jié)T細胞和耗竭性CD8 T細胞在肝癌的克隆富集
· LAYN基因與腫瘤調(diào)節(jié)T細胞以及耗竭性CD8 T細胞的抑制功能有關(guān)
· 結(jié)合表達和基于TCR的分析顯示T細胞子集間的連通性

SUMMARY
在癌癥里面,，對腫瘤浸潤淋巴細胞（TILs）的系統(tǒng)性問診是發(fā)展免疫療法和預(yù)測臨床反應(yīng)的關(guān)鍵。文章作者從六個肝癌患者的外周血,、腫瘤,、癌旁組織分離出5063個T細胞，做了深度單細胞RNA測序,。這些單個細胞的轉(zhuǎn)錄表達譜結(jié)合有T細胞受體（TCR）序列,，根據(jù)分子和功能性質(zhì)鑒定出11個T細胞子集，描繪了它們的發(fā)展軌跡,。特定的子集比如耗竭性CD8^+ T細胞和調(diào)節(jié)T細胞在肝癌中優(yōu)先富集,，可能存在克隆擴增。他們鑒定了每個子集的信號基因,，其中一個基因,，layilin（LAYN），在激活的CD8^+ T細胞和調(diào)節(jié)T細胞中上調(diào),，抑制CD8^+ T細胞在體外的功能,。這些轉(zhuǎn)錄組數(shù)據(jù)的概況為理解癌癥的免疫環(huán)境提供了有價值的見解和豐富的資源。

研究背景
癌癥免疫療法在過去數(shù)十年里極大促進了腫瘤的治療,。但是這些療法對于不同的癌癥患者以及癌癥類型,，效率并不一致。目前最重要的一點是找到預(yù)測治療反應(yīng)的穩(wěn)定生物學標記,，這需要深入理解腫瘤T細胞,，特別是導致CD8^+ T細胞耗竭和調(diào)節(jié)T細胞積累的機制和通路，才能為協(xié)調(diào)免疫系統(tǒng)消滅癌癥提供更好的策略,。
浸潤淋巴細胞的單細胞分析使得在高度復雜的腫瘤微環(huán)境中了解到這些細胞的細節(jié)成為可能,。scRNA-seq是一個很強大的工具，可以明確定義出每個細胞的TCR序列,，從而識別腫瘤特異的新抗原,。盡管TCR庫相當龐大，鑒定TCR序列可以解釋T細胞的克隆擴張模式和相關(guān)譜系,。
肝癌是致死人數(shù)最高的幾種癌癥之一,，在中國,，感染HBV是引發(fā)肝癌的主要原因。目前缺乏免疫治療的成功案例,，治療選擇有限,。肝癌腫瘤導致TILs水平的顯著提高，經(jīng)推測這些TILs已經(jīng)喪失了殺死腫瘤細胞的能力,。因此,，肝癌腫瘤的研究可以提供一個描述TILs調(diào)節(jié)異常的絕佳模型。

實驗&分析

Study Overview

實驗方法

對于6個未經(jīng)治療的肝癌患者,，首先做了外顯子組DNAseq和大量細胞RNAseq，尋找CNV和體細胞變異,，發(fā)現(xiàn)基因模式跟已有的肝癌研究是吻合的,。并且進一步確認了存在不同的T細胞亞型。
然后使用FACS分選出感興趣的細胞群體,，CD8+和CD4+,，做單細胞RNA-seq。

單細胞收集

外周血單個核細胞使用HISTOPAQUE-1077解決方案,。
新鮮的腫瘤和鄰近正常組織樣本先切成1立方毫米的小塊,，經(jīng)研磨過濾懸浮去紅重懸浮得到。
我以前處理過很多外周血和組織,，去除紅細胞的步驟是很重要很重要滴

單細胞分類,、反轉(zhuǎn)錄和擴增

使用FACS將細胞分選到有裂解液的96孔板里面，一個板子上都是同一類的細胞,，裂解液里面包含了dNTP,，OligodT引物以及加了RNA酶抑制劑的TritonX-100。分選后的板子凍存在-80℃,，有兩個病人還加了外源RNA質(zhì)控,。單細胞轉(zhuǎn)錄組擴增是按照Smart-seq2的流程走的，cDNA純化用Beckman的XP磁珠,。在第一輪磁珠純化后,，每個細胞的cDNA做CD3D qPCR和片段分析。
鑒定出高質(zhì)量的單細胞,，進一步純化去除小片段,，Qubit做濃度質(zhì)檢，使用Vazyme的試劑構(gòu)建文庫,?；颊逷0205、P0508 和 P0322使用HiSeq2500 測序（100bp讀長）,，患者P0407 和P1116使用Hiseq 4000測序（150bp讀長）,。而對于患者P1202,，每個單細胞是人工用吸管挑到孔里并且轉(zhuǎn)錄組擴增用的是Tang2010方法，因為明顯的3'端偏倚,，這些細胞的TCR序列不會組裝到單細胞轉(zhuǎn)錄組數(shù)據(jù)上,，該患者沒有做大量細胞測序。
測序的實驗設(shè)計很周到,，這種算珍貴樣本,，提前考慮一些可能出現(xiàn)的異常數(shù)據(jù)

大量細胞的DNA/RNA分離和測序

使用QIAGEN的Mini Kit提取外周血和組織樣本的基因組DNA，Qubit測濃度,，瓊脂糖電泳質(zhì)檢,。使用Agilent的試劑構(gòu)建外顯子文庫，Hiseq4000測序,。
對于RNA,，腫瘤和鄰近正常組織經(jīng)手術(shù)切除放在冰上，先用 RNAlater RNA stabilization reagent 保存,。用QIAGEN的RNeasy Mini Kit釋放出RNA,，NanoDrop測濃度，AATI片段分析質(zhì)檢,。然后使用NEB的試劑構(gòu)建文庫,，Hiseq4000測序。

免疫組化

人體組織標本經(jīng)過相關(guān)審核批準由北京世紀壇醫(yī)院提供,。腫瘤切除后福爾馬林處理48小時,，之后30分鐘內(nèi)收集，脫水包埋使用常規(guī)方法,。

離體T細胞激活和FACS分析

培養(yǎng)外周血單個核細胞（PBMCs）,，在指定的時間點采集細胞并染色，使用LSR-II分析儀和Flowjo軟件分析每一個存活/單線態(tài)亞群的LAYN基因和T細胞激活/耗竭標記的表面表達,。所有數(shù)據(jù)來自至少三組獨立實驗和6個不同的PBMCs源,。

LAYN在原代CD8^+ T細胞的過表達

把人layilin基因克隆到逆轉(zhuǎn)錄病毒載體，pAmpho包裝,，感染從人PBMCs分離的T細胞,，免疫磁珠激活。使用FACS分離GFP+的人layilin過表達CD8^+ T細胞,。

定量和統(tǒng)計分析

測序深度
5063個單細胞,，平均每個細胞129w單一匹配讀長（uniquely mapped reads）。
Smart-seq和Smart-seq2的原始文章建議在300w以上,，不同的項目需求不同,，一般超過100w就比較穩(wěn)定了

為了說明在這種測序深度下的檢測結(jié)果是可信賴的，做了一個飽和分析,，主要是看重要的基因類別,，特別是細胞激素和轉(zhuǎn)錄因子,。

單細胞RNAseq數(shù)據(jù)處理

從Hiseq2500/4000下機的數(shù)據(jù)，首先過濾掉低質(zhì)量的reads：1）有10%讀長的堿基無法識別,；2）50%讀長的堿基質(zhì)量小于5,；3）含有接頭序列。過濾后通常能剩下90%的原始reads,，先使用GSNAP工具把這些reads比對到RFam下載的核糖體RNA序列,。
沒有匹配上的reads（clean reads）再拿去比對參考序列，包括UCSC下載的hg19和spike-in的質(zhì)控RNA,。另外,，從UCSC下載了一個基因模型文件knownGene.txt，用可調(diào)參數(shù)“–novelsplicing 1 -n 10 -i 1 -M 2”來定義已知的exon-exon junctions,，其他參數(shù)使用默認設(shè)置,。使用R的函數(shù)findOverlaps檢索跟已定義基因重合的reads，得到每個基因的讀長計數(shù),。
使用DESeq方法對讀長計數(shù)作進一步的標準化。同一個患者,，首先計算每一個基因在所有細胞的讀長計數(shù)的幾何平均數(shù),，得到讀長計數(shù)與其幾何平均數(shù)的比值（只計算非零項）。然后得到比值的中位數(shù)作為每個細胞的大小因數(shù)（size factor）,，接著把標準化的讀長計數(shù)（$RC_{normalized}$ ）定義為RC/size-factor,。$log_2 (RC_{normalized} +1)$的值用于下游的分析，比如PCA和聚類,。
還使用下面的步驟過濾掉低質(zhì)量的細胞：1）去掉單一匹配讀長少于25w的細胞,；2）去掉CD3D基因的RPKM小于1.0的細胞；3）去掉FACS分選的CD8^+ T細胞里CD8A基因的RPKM小于1.0或者CD4基因的RPKM大于10的細胞,；4）去掉FACS分選的CD4^+ T細胞里CD4基因的RPKM小于1.0或者CD8A基因的RPKM大于10的細胞；5）去掉離群點分析(outlier analysis)失敗的細胞,。
3)和4)是不信任FACS么，還是別的原因,？

在離群點分析中,，PCA分析使用了$log_2 (RC_{normalized} +1)$,，擬合分析使用了方差平方系數(shù)（CV^2 ）、每個FACS-分選的T細胞亞型的$RC_{normalized}$ 的平均值,，手動檢查結(jié)果以確定異常細胞,。典型的異常細胞定位在其他細胞很遠的地方,，通常對主成分有很高的貢獻。此外,，異常細胞也表現(xiàn)出很低的大小因數(shù),，表示細胞里的RNA含量很少，通過擬合曲線可以更好地移除這些細胞,。過濾掉異常細胞之后,，6個患者總共剩下4128個細胞，其中5個HBV陽性患者共有4070個細胞,。

大量細胞的轉(zhuǎn)錄組和外顯子組數(shù)據(jù)處理

大量細胞的RNA-seq數(shù)據(jù)處理跟單細胞的RNA-seq數(shù)據(jù)使用同一套GSNAP流程,。對于外顯子測序數(shù)據(jù),，首先過濾掉低質(zhì)量的reads,。高質(zhì)量的reads使用bwa-0.7.12回貼到人的參考序列（從Broad下載，加了誘餌序列的b37版本）,。根據(jù)建議的最佳做法使用GATK作比對后處理,，包括bam文件分類、duplicate reads標記,，候選INDEl附近的reads本地重排,，以及堿基質(zhì)量重新校正。使用Strelka的默認參數(shù)分析體細胞突變,，包括SNV和INDEL。那些所謂的突變特別是低頻突變,，使用IGV工具作進一步的手動檢查,。用ADTex得到體細胞拷貝數(shù)變化和雜合缺失（LOH）。根據(jù)與已有的驅(qū)動基因（driver gene）名單的比較,，進一步得到這些患者的候選驅(qū)動突變（driver mutations）。一些已知的突變熱點：TP53 G245C (P0205), PIK3CA H1047R(P0508), PTEN C136R (P1116), TP53 R249S P0322) ,，TP53 I195F (P0427),。后期患者 P0322有兩個不同的TP53突變。

TCR分析

使用TraCeR給單細胞RNA-seq數(shù)據(jù)的每一個T細胞組裝上TCR序列,，進而鑒定出CDR3序列,、重排的TCR基因以及它們的表達豐度（TPM）。首先,，去掉沒有明顯TCR組成的細胞。然后按照后面的步驟分別整理TCR的α和β鏈,，第一條TCR的α（β）鏈定義為：1）所有單個T細胞中,，保證只能存在一條有用的α（β）鏈,。2）如果一個T細胞中不只一條α（β）鏈，保留只有一個α（β）鏈的主要形態(tài)被檢測到的細胞,。通常來講，一條α（β）鏈是有用的而其他的都沒用,，或者說其中一條的表達水平遠遠高于替代等位基因（alternative allele）,，以此可以識別出來,。接下來，根據(jù)表達分布的直方圖分析

過濾掉那些第二條TCR的α鏈TPM小于10或者β鏈TPM小于15的細胞。從成功組裝TCR序列的4032個細胞里面,，找到3792個細胞有TCR的α/β對,。

非監(jiān)督聚類

每個患者的給定基因的相對表達量通過Y-mean（Y）計算,，Y是標準化的表達量$log_2 (RC_{normalized} +1)$,，mean(Y)是這個患者所有細胞的Y的平均值,。這個過程移除了由患者差異引起的不想要的變異,。這種集中數(shù)據(jù)隨后結(jié)合所有的患者。前n個高標準差的基因被認為是高變異基因,，試驗過不同的n值：1500,、2000、2500和3000。所有基因產(chǎn)生的集中表達數(shù)據(jù)用于改進的SC3聚類流程,。首先創(chuàng)建基于相關(guān)分析（Spearman correlation）的距離矩陣,，然后計算拉普拉斯算子的特征向量。接著在最初的d特征向量（d取不同的值,，從細胞總數(shù)的4%到7%）多次運用k-means算法,。平均不同次運行的結(jié)果得到一個一致性矩陣（consensus matrix），從中可以看出兩個細胞在不同次運行中聚到一起的頻率是多少,。最后在一致性矩陣執(zhí)行完整的層次聚類（hierarchical clustering with complete agglomeration），得到k群（k clusters）,。用于k-means和層次聚類的SC3的參數(shù)k,，從2迭代到10,。每一次運行SC3,，都要計算剪影（silhouette）,，繪制一致性矩陣，鑒定集群特異基因（cluster specific genes）,。所有這三個方面可以幫助研究者根據(jù)經(jīng)驗得到最優(yōu)的k和n。一旦確定出穩(wěn)定的集群,，使用的程序?qū)⒈坏矫恳粋€集群,，用以揭示每個集群中細胞間的高變異基因，然后就可以使用這些變異基因來識別亞群,。
中間數(shù)據(jù)要用眼睛去看,，要根據(jù)經(jīng)驗去選，跑一個分析流程遠不是那么簡單

得到11個穩(wěn)定的集群,，包括5個CD8^+ 集群和6個CD4^+ 集群,，每一個集群都有獨特的標記基因，CD8^+ 在不同患者間可以比較,。

得到聚類結(jié)果之后,，用ANOVA的R函數(shù)aov識別差異表達基因（differential expressed genes），再用R函數(shù)TukeyHSD的杜凱氏方法（Tukey's range test）實現(xiàn)每個集群對的差異檢驗,。每個集群的差異表達基因符合這些標準：1）FDR調(diào)整后的F檢驗P值 < 0.01,；2）研究的集群與其他集群間的log2倍數(shù)變化絕對值="">1；3）比較 研究的集群與其他集群 的杜凱氏檢驗的P值<>
此外,，細胞周期也被認為是細胞群體分析的重要混雜因素,。用ScLVM檢查了細胞周期的貢獻度，發(fā)現(xiàn)對這里的聚類結(jié)果影響非常小,，細胞周期因素影響的方差只占生物學因素的4.5-13.6%（患者P0508, P0322, P0407, P1116）和37.7%（患者P0205）,，后面這個也遠低于之前報道的85%。聚類結(jié)果確實包含了一些細胞周期相關(guān)的基因,，像TUBB4B,、TUBA1B,、MCM7、STMN1和MKI67,，但這些只有62個細胞,，不是主要的亞群，也沒有進行tSNE可視化以及差異基因表達分析,。

基于綜合 表達和TCR克隆 的CD8^+ 和CD4^+ T細胞的狀態(tài)轉(zhuǎn)換分析

CD8^+ /CD4^+ T細胞在二維狀態(tài)空間的擬時間（pseudotime）排列使用Monocle2定義,。細胞的順序由CD8^+ /CD4^+ T細胞群（無MAIT）的最分散基因的表達推斷而來。每個點代表一個細胞,，每種顏色代表一個集群,。

CD8^+ /CD4^+ T細胞集群的細胞狀態(tài)轉(zhuǎn)換由共有的TCR推斷而來。只有在不同細胞集群共享TCR的才會著色,，每種顏色代表一種獨特的TCR克隆類型,，線條連接基于不同的TCR共享程度，線條厚度代表共享的TCR數(shù)目,。插入圖片顯示兩個集群的共享TCR數(shù)目,，條形中的兩種顏色代表特別的兩個集群共享的TCR。

T細胞耗竭狀態(tài)分析

為了揭示耗竭性T細胞的特征,，首先移除了一些腫瘤浸潤T細胞,，它們與正常組織的T細胞帶有相同的TCR序列。MAIT細胞因為特殊的特性也被移除掉,。根據(jù)SC3的聚類分析,，C4_CD8-LAYN 集群中的腫瘤浸潤CD8^+ T細胞被定義為耗竭性T細胞，其他的則是非耗竭性,。耗竭性與非耗竭性T細胞間的差異表達基因檢測使用R包limma的線性模型和經(jīng)驗貝葉斯方法實現(xiàn),，使用嚴格的顯著性閾值，調(diào)整后的P值（BH多重檢驗校正）<>降低假陽性找到可信的耗竭性標記,。使用limma是因為它對高度平行的基因組數(shù)據(jù)的杠桿作用,，可以借用gene-wise模型間的信息。計算每個患者每個DEG的相對表達（log2 fold-change）,，計算每個患者所有DEG的相對表達均值,，從而得到每個患者的耗竭性分數(shù)（per-patient exhaustion score）。耗竭性分數(shù)隨著臨床分期靠后而增加,。腫瘤調(diào)節(jié)T細胞特異的基因使用同樣的方法檢測,，但是采用不同的閾值（adjusted p value<0.01，fold>

TCGA數(shù)據(jù)分析

腫瘤基因組圖譜（TCGA）是美國NCI和NHGRI的聯(lián)合項目,。
TCGA數(shù)據(jù)用來檢測選取的基因與患者存活率的關(guān)聯(lián),，存活分析使用R包survival實現(xiàn),。

文章的raw data 和 processed data分別是EGA: EGAS00001002072 和 GEO:GSE98638。
http://hcc. 有這些數(shù)據(jù)的分析結(jié)果。 

讀了這篇文章的收獲是知道了具體的單細胞項目可以是這樣的,。單細胞測序也不是完全推新的技術(shù),，是現(xiàn)有技術(shù)的改進和細化。很多分析的套路都是延用的,，另需要照顧單細胞的數(shù)據(jù)量和異質(zhì)性,。怎么設(shè)計質(zhì)控，丟掉哪些數(shù)據(jù),，使用哪些數(shù)據(jù),，如何整合出一個可用的模型呢？ 
菜:bird:表示還有n多基礎(chǔ)的東西要刷~

點擊可以加入單細胞數(shù)據(jù)處理學習交流小組