文獻詳解欄目
每個人的時間精力有限，必須優(yōu)先閱讀相關(guān)文獻，開設(shè)這個欄目也是希望為大家推薦高質(zhì)量的單細胞相關(guān)文獻。如果大家對單細胞轉(zhuǎn)錄組感興趣可以關(guān)注一下，哪怕每天只學(xué)一點點，積土成山，積水成淵。
當(dāng)然一個人的力量終歸是小的，我也希望匯聚一群人，形成一個場，這里頭最重要的生產(chǎn)力不是單個人多聰明，多厲害，而是每個人相互作用，形成的那個氛圍。
希望大家能有所收獲！
文章信息
文章由北醫(yī)三院付衛(wèi)團隊、喬杰團隊和北大生命科學(xué)學(xué)院湯富酬團隊于2018年11月合作發(fā)表在Science上，文章題目是：Single-cell multiomics sequencing and analyses of human colorectal cancer
這一篇理解難度稍大，當(dāng)做背景知識的了解，至于其中數(shù)據(jù)的分析細節(jié)文中沒有提及
導(dǎo)語
之前一次只能研究單細胞層面的基因組、轉(zhuǎn)錄組或DNA甲基化組其中一種，不能在一個細胞中同時研究多個組學(xué)，2016年湯富酬研究組將三重組學(xué)研究方法scTrio-seq（single-cell triple omics sequencing technique）發(fā)表在Cell Research上。2018年11月北醫(yī)三院付衛(wèi)團隊、喬杰團隊和北大生命科學(xué)學(xué)院湯富酬團隊又在Science合作發(fā)表了Single-cell multiomics sequencing and analyses of human colorectal cancer
文章優(yōu)化了單細胞多組學(xué)測序技術(shù)，并對原發(fā)瘤、淋巴轉(zhuǎn)移和遠端轉(zhuǎn)移區(qū)域分別采樣，首次從單細胞分辨率解釋了結(jié)直腸癌的發(fā)生與轉(zhuǎn)移過程中中基因組拷貝數(shù)變異、DNA甲基化異常及基因表達改變的特點，證明了用單細胞多組測序重建遺傳譜系和追蹤其表觀基因組和基因表達動力學(xué)的可行性
背景知識
淋巴轉(zhuǎn)移（lymphatic metastases）和遠端轉(zhuǎn)移(distant metastases)
癌癥可以以2種方式出現(xiàn)在淋巴結(jié)中：一種是從淋巴結(jié)形成的腫瘤叫淋巴瘤，另一種是從其他部位擴散叫淋巴轉(zhuǎn)移（更為常見）。淋巴轉(zhuǎn)移一般會和乳腺癌、前列腺癌、肺癌、結(jié)直腸癌的不良預(yù)后相關(guān)，淋巴轉(zhuǎn)移雖然不是致死因素，但會導(dǎo)致癌細胞擴散到重要器官。遠端轉(zhuǎn)移也叫惡性轉(zhuǎn)移，腫瘤細胞從原始發(fā)生的部位借由侵入循環(huán)系統(tǒng)，轉(zhuǎn)移到身體其他部位繼續(xù)生長，幾乎不可能使用外科手術(shù)切除根治，多半只能用大范圍循環(huán)全身的放射治療或化療等手段來抑制已轉(zhuǎn)移的癌細胞。
2017年發(fā)表在Science的Origins of lymphatic and distant metastases in human colorectal cancer 中描述了結(jié)直腸癌之前的腫瘤擴散的TNM層級是：primary tumor(T)=》lymph node system(N)=》distant metastases(M)，但是有臨床證據(jù)表明移除淋巴結(jié)并不會提高病人存活率，因此N和M之間的關(guān)系可能并不是簡單的上下級。文章發(fā)現(xiàn)淋巴結(jié)和遠端器官存在獨立起源的證據(jù)，只是淋巴結(jié)轉(zhuǎn)移機制形成的更快，因此淋巴轉(zhuǎn)移形成更早，發(fā)生更頻繁。因此作者不推薦直接假設(shè)淋巴結(jié)會引發(fā)遠端轉(zhuǎn)移而直接切除(https://www.medscape.com/viewarticle/882502)
結(jié)直腸癌colorectal cancer(CRC) ：是結(jié)腸癌和直腸癌的統(tǒng)稱，是消化道惡性腫瘤之一。2018年Cancer Stats統(tǒng)計顯示：結(jié)直腸癌在男性中發(fā)病率第2死亡率第3，女性發(fā)病率第4，死亡率第3。
約95％的結(jié)直腸癌是由結(jié)腸和直腸內(nèi)壁的腺細胞發(fā)展而來，癌癥通常開始于內(nèi)壁最內(nèi)層，并緩慢生長到外層。(http://www.cancer.org/Cancer/ColonandRectumCancer/DetailedGuide/colorectal-cancer-what-is-colorectal-cancer)
結(jié)直腸癌發(fā)生率在40歲開始增加，60～75歲達高峰，結(jié)腸癌在女性患者較常見；直腸癌在男性患者常見；大約5%的結(jié)腸癌或直腸癌患者在結(jié)腸和直腸有兩個或更多病灶，并非簡單從一個病灶轉(zhuǎn)移至另一個所致(https://www.msdmanuals.com/)
結(jié)直腸癌典型的分子特征是：基因組不穩(wěn)定性、表觀遺傳學(xué)異常、基因表達紊亂
瘤內(nèi)異質(zhì)性 Intratumoral heterogeneity (ITH) ：惡性腫瘤的特征之一，腫瘤異質(zhì)性包括腫瘤間異質(zhì)性（不同腫瘤細胞之間的基因與表型不同）和腫瘤內(nèi)異質(zhì)性（相同腫瘤細胞以內(nèi)的基因與表型也不同），其中腫瘤內(nèi)異質(zhì)性又包括空間異質(zhì)性（相同腫瘤不同區(qū)域不同，如未擴增的細胞背景中有成簇擴增細胞；少量擴增背景中有未擴增的細胞；孤立的細胞擴增【利用多位點取樣方案或者tissue microarrays (TMAs) 調(diào)查】）與時間異質(zhì)性（原初腫瘤與次生腫瘤不同）。
異質(zhì)性的產(chǎn)生是因為同一腫瘤由多種不同基因組特征的細胞組成，每一種細胞構(gòu)成一個亞克隆(subclone)。腫瘤組織會存在對治療藥物有抗性的亞克隆，但比例不高。當(dāng)治療的藥物除去敏感的亞克隆時，抗性的亞克隆細胞會不受藥物抑制并且少了空間競爭，因此會加快生長速度，導(dǎo)致腫瘤復(fù)發(fā)或者發(fā)生轉(zhuǎn)移，而且轉(zhuǎn)移后的亞克隆對同種治療方案也會有抗性。因此研究腫瘤細胞的亞克隆以及不同的亞克隆的轉(zhuǎn)移是一個熱點，尤其是亞克隆是如何從原位癌轉(zhuǎn)移到其他臟器而形成轉(zhuǎn)移癌。2017的一篇文章專門研究了轉(zhuǎn)移癌亞克隆與原位癌亞克隆的進化關(guān)系 ，他們發(fā)現(xiàn)結(jié)直腸癌腫瘤轉(zhuǎn)移癌高深度測序就可以找到腫瘤的大部分基因組變異，另外驗證了"轉(zhuǎn)移癌多克隆起源說(轉(zhuǎn)移癌是由多個起源于原位癌的亞克隆發(fā)展而來，而非由單個細胞發(fā)育)"，發(fā)現(xiàn)了結(jié)直腸癌的淋巴與遠端并行轉(zhuǎn)移 。
方法
scTrio-seq2技術(shù)：用于somatic copy number alterations (SCNAs)拷貝數(shù)變異、DNA甲基化特征、細胞連續(xù)的轉(zhuǎn)錄信息；整合了單細胞重亞硫酸鹽測序（scBS-seq）用于全基因組甲基化分析；研究的細胞數(shù)量從之前的25個增至1900個
分析了12個CRC患者（III期或IV期）的約1900個單細胞，7.6Tb高質(zhì)量測序數(shù)據(jù)。DNA甲基化研究中平均每個細胞測序量為4.1Gb，平均覆蓋到全基因組內(nèi)870多萬CpG位點；轉(zhuǎn)錄組研究中每個細胞測序量為0.9Gb，平均覆蓋3700多個基因
多區(qū)域采集了10個患者的原發(fā)瘤、淋巴轉(zhuǎn)移瘤或遠端轉(zhuǎn)移瘤樣本(利用兩種不同來源的細胞，可以發(fā)現(xiàn)每個患者因突變而產(chǎn)生的遺傳譜系)
文章實驗流程
圖A是取樣：治療前后的腫瘤區(qū)域（包括原發(fā)瘤primary tumor，PT；淋巴結(jié)轉(zhuǎn)移位 lymph node metastasis, LN；肝轉(zhuǎn)移位 liver metastasis, ML；化療后肝轉(zhuǎn)移位 posttreatment liver metastasis, MP），然后測序分析了基因組、轉(zhuǎn)錄組、甲基化組；
圖B是化療6個周期后的患者CRC01取樣：一共取了ML（4個）、MP（5個）、LN（3個）、PT（4個）共16個腫瘤區(qū)域
單細胞甲基化數(shù)據(jù)處理
首先raw reads去接頭、引物、低質(zhì)量堿基，然后利用Bismark（V0.7.6）clean reads比對到hg19基因組，PCR重復(fù)利用samtools rmdup(V 0.1.18)去除，數(shù)據(jù)統(tǒng)計(比對數(shù)、比對率、CpG位點數(shù)、亞硫酸氫鹽轉(zhuǎn)化率等)[其中亞硫酸氫鹽轉(zhuǎn)化率是由非甲基化的lamda DNA的spike-in計算的]，CpG位點小于200萬個或亞硫酸氫鹽轉(zhuǎn)化率小于98.5%的細胞被排除。僅使用甲基化水平大于0.9或小于0.1的CpG位點進行總體甲基化水平計算。
啟動子區(qū)設(shè)定為轉(zhuǎn)錄起始位點的上游1 kb到下游0.5 kb。為了計算RefSeq基因各基因體的DNA甲基化水平，將每個基因體劃分為100個等分，并將其上下游側(cè)翼區(qū)域(15 kb)分別劃分為10個等分?；蚪M注釋信息從UCSC獲取。利用bedtools(V 2.17.0)和自定義腳本(沒放Git鏈接)計算平均甲基化水平，設(shè)置滑動窗口大于等于3個CpG位點
根據(jù)甲基化測序數(shù)據(jù)估計拷貝變異數(shù)
主要基于Garvin等人開發(fā)的Ginkgo算法 ，基因組被分成10856個不等長的bins，長度中位數(shù)為250kb，并根據(jù)算法的過濾器排除了一些異常的bins。
BED文件是利用bedtools從BAM文件得到。每個bin的read counts值利用所有bins的count平均值進行標(biāo)準(zhǔn)化，采用低水平均一化(Lowess normalization)來校正基因組GC含量的偏差。另外，以正常的二倍體細胞作為對照，減少scRS-seq的其他誤差。利用Circular binary segmentation (CBS) 對copy number文件進行分隔，參數(shù)為"alpha = 0.0001"和"undo.prune = 0.05" 。CBS分隔后，每一段的所有bins的計數(shù)重置為這一段的bin count的中位數(shù)。每個單細胞的基本拷貝數(shù)由smallest sum-of- squares (SoS) error和 scaled copy-number profile (SCNP)決定，其中SCNP又進一步四舍五入，最后得到了整數(shù)值copy-number profile (FCNP)。利用 GISTIC2.0 (https://software.broadinstitute.org/cancer/cga/gistic) 鑒定了重要的SCNAs和潛在的基因靶點。
TCGA數(shù)據(jù)分析
從https://tcga-data.nci.nih.gov/docs/publications/coadread_2012/獲得已發(fā)表的人類CRC的SCNA片段數(shù)據(jù)，用于SCNA頻率統(tǒng)計(不包括X染色體)。將本文的數(shù)據(jù)與TCGA的進行比較時，本文研究的CMS3類型的患者CRC02是找不到對應(yīng)的，只有Affymetrix SNP 6.0芯片得到的178個non-hypermutated樣本。SCNA譜進一步轉(zhuǎn)變成長度不等的bins用于單細胞SCNA統(tǒng)計，拷貝數(shù)為>2.5的bin表示擴增，小于1.5表示缺失。Circos圖是根據(jù)https://github.com/venyao/shinyCircos制作的。
WGS數(shù)據(jù)處理
raw reads =》trimmed =》BWA mem(V 0.7.12) 比對到hg19 =》 samtools sort =》Picard(V1.139) merge BAM文件 + 標(biāo)記重復(fù) =》BAM文件用GATK(V 3.4-46) 預(yù)處理=》muTect (V 1.1.4) call SNVs ，自己腳本過濾=》取每個患者的外周血或鄰近正常組織作為對照(somatic variants)，在線Venn圖做出SNV數(shù)量
單細胞RNA-seq處理
利用湯教授自己的方法得到的數(shù)據(jù)處理：預(yù)處理過程都一樣，然后用STAR(V 2.5.0) 2步比對到hg19，Cufflinks(V 2.2.1) 使用默認參數(shù)進行FPKM定量；
利用multiplexed scRNA-seq方法得到的數(shù)據(jù)，先利用read2的barcode信息將reads分配到每一個細胞，每個細胞中read2對應(yīng)的read1利用read ID分隔，read1中的TSO序列 利用自己的腳本過濾掉，然后利用Tophat(V2.1.1) 單端比對到hg19，利用UMI實現(xiàn)TPM標(biāo)準(zhǔn)化(文中說道：大部分的表達量都使用log2(FPKM + 1)或者log2(TPM/10 + 1) ) ，然后統(tǒng)計了mapped read numbers, mapping ratios, RefSeq gene numbers等信息，根據(jù) FPKM > 1 or TPM > 1去除了比對率 < 20%或有表達量的Refseq基因數(shù)量 < 1500
結(jié)果
Genomic alterations in tumors provide markers for lineage tracing. 克隆變異出現(xiàn)在腫瘤早期階段，亞克隆拷貝數(shù)變異標(biāo)志著亞型的出現(xiàn)
結(jié)直腸癌患者單個癌細胞的染色體拷貝數(shù)變異譜+高精度的染色體內(nèi)斷點信息=》譜系追蹤=》原發(fā)位腫瘤（PT）的亞克隆結(jié)構(gòu)通常比轉(zhuǎn)移位腫瘤更復(fù)雜
結(jié)果得到了5個患者的90個細胞以上的甲基化數(shù)據(jù)，細胞被分成了不同的基因亞型
其中，CRC01基于21個亞克隆的拷貝數(shù)斷點，鑒定了來自2個不同譜系（A、B）的12個亞型，其中每個亞型都有4-8個亞克隆的斷點(斷點的上下位置和拷貝數(shù)變異數(shù)增加、減少對應(yīng))，A5亞型同時出現(xiàn)在了肝轉(zhuǎn)移位和淋巴轉(zhuǎn)移位，表明這兩種移位有共同起源，這5個病人的原癌亞克隆結(jié)構(gòu)比其他轉(zhuǎn)移類型更復(fù)雜
CRC01患者癌細胞的單細胞染色體拷貝數(shù)變異譜：
CRC01患者亞克隆拷貝數(shù)斷點：
CRC01患者亞克隆結(jié)構(gòu)：
結(jié)腸直癌細胞的DNA甲基化水平要低于癌旁的正常上皮細胞，同一腫瘤組織中同一譜系的甲基化程度相近，不同譜系出現(xiàn)差別。低甲基化基因組區(qū)域顯著富集在LTR（long terminal repeats）、LINE-1 (long interspersed nuclear elemnt 1)和異染色質(zhì)區(qū)域(H3K9me3)，而高甲基化的基因組區(qū)域顯著富集在CpG島、H3K4me3和開放染色質(zhì)區(qū)域。
以CRC01患者為例：
甲基化的異質(zhì)性主要來自同一個患者腫瘤內(nèi)不同亞克隆之間的DNA甲基化差異，而不是同一個亞克隆內(nèi)部不同細胞間的DNA甲基化差異
啟動子區(qū)域的甲基化與相應(yīng)基因的表達呈顯著的負相關(guān)，而基因區(qū)的甲基化與相應(yīng)基因的表達呈正相關(guān)
【The gray lines represents individual cells. The blue line represents the mean value for each patient. TSS, transcription start site; TES, transcription end site.】
同一患者同一譜系的腫瘤細胞從原發(fā)灶到轉(zhuǎn)移灶全基因組DNA甲基化水平基本穩(wěn)定，組內(nèi)局部區(qū)域可能會有比較大的波動
每個亞型內(nèi)的去甲基化 程度是一致的，不同亞型程度不同 ；
正常上皮細胞的基因組區(qū)域甲基化越高，它就越容易發(fā)生去甲基化；
去甲基化的程度與基因組重復(fù)序列L1（long interspersed nuclear elemnt 1）以及癌旁正常組織中H3K9me3修飾的密度呈正相關(guān)，與H3K4me3標(biāo)記和正常組織的開放染色質(zhì)區(qū)域密度呈負相關(guān)；
有趣的是，L1比LINE-2更活躍，在所有病人的癌細胞中顯示了更強的去甲基化能力，這個與胚胎發(fā)育中情況相反(胚胎發(fā)育過程中L1一般比L2去甲基化能力弱) ，說明在腫瘤發(fā)生與發(fā)展過程中，L1和異染色質(zhì)區(qū)域產(chǎn)生了異常的去甲基化過程，打破了正常的發(fā)育規(guī)律
癌細胞相比于癌旁細胞的DNA甲基化水平：
結(jié)直腸癌細胞中6條染色體（4號、5號、8號、13號、18號、和X染色體）傾向于發(fā)生更強烈的DNA去甲基化，其中三條低甲基化染色體（8、13和18）在TGCA和研究的患者中都有較高的拷貝數(shù)變異。結(jié)合WGS結(jié)果發(fā)現(xiàn)，有5條第甲基化的染色體(4號、5號、8號、13號、和X染色體)的單核苷酸變異(SNVs)發(fā)生顯著富集
