【單細胞入門必讀】單細胞RNA-seq技術(shù)和數(shù)據(jù)分析簡介

Single-Cell RNA-Seq Technologies and Related Computational Data Analysis

單細胞RNA測序(scRNA-seq)技術(shù)允許在單細胞分辨率下解析基因表達辽故,這極大地改變了轉(zhuǎn)錄組學(xué)研究掷漱。目前已經(jīng)開發(fā)了大量的scRNA-seq protocol怕吴,這些protocol各有特點赞哗,各有優(yōu)缺點羊赵。由于技術(shù)限制和生物因素,scRNA-seq數(shù)據(jù)比 bulk RNA-seq數(shù)據(jù)更嘈雜裳凸、更復(fù)雜碌奉。scRNA-seq數(shù)據(jù)的高度可變性給數(shù)據(jù)分析帶來了計算方面的挑戰(zhàn)。雖然越來越多的生物信息學(xué)方法被提出用于分析和解釋scRNA-seq數(shù)據(jù)碍论,但需要新的算法來確保結(jié)果的準確性和再現(xiàn)性谅猾。在本文中,我們概述現(xiàn)有的一些protocol和scRNA-seq單細胞分離技術(shù),并討論不同方法scRNA-seq數(shù)據(jù)分析包括質(zhì)量控制、read mapping, 基因表達量化鳍悠、批次效應(yīng)校正,規(guī)范化,imputation,降維,特征選擇,細胞聚類税娜、軌跡推理,差異表達calling,可變剪接,等位基因的表達,基因調(diào)控網(wǎng)絡(luò)重構(gòu)。此外藏研,我們概述了scRNA-seq技術(shù)的發(fā)展和應(yīng)用前景敬矩。

Introduction:

在過去的十年中,Bulk RNA-seq技術(shù)被廣泛應(yīng)用于研究群體水平上的基因表達模式遥倦。單細胞RNA測序(scRNA-seq)的出現(xiàn)為探索單細胞水平的基因表達譜提供了前所未有的機會谤绳。目前占锯,scRNA-seq已經(jīng)成為研究細胞異質(zhì)性和早期胚胎發(fā)育等關(guān)鍵生物學(xué)問題的良好選擇,因為bulk RNA-seq主要反映了數(shù)千個細胞的平均基因表達缩筛。近年來消略,scRNA-seq已被應(yīng)用于不同的物種,特別是不同的人體組織(包括正常和癌癥)瞎抛,這些研究揭示了有意義的細胞間基因表達變異.近年來艺演,隨著測序技術(shù)的不斷創(chuàng)新,一些不同的scRNA-seq protocal被提出桐臊,極大地促進了對單細胞動態(tài)基因表達的理解胎撤。其中一種是高效的LCM-seq (Nichterwitz et al., 2016),它結(jié)合了激光捕獲顯微鏡(LCM)和Smart-seq2 (Picelli et al., 2013)断凶,用于單細胞轉(zhuǎn)錄組學(xué)伤提,不用分離組織 。目前可用的scRNA-seq protocol 根據(jù)捕獲的轉(zhuǎn)錄覆蓋率主要可分為兩類:(i)全長轉(zhuǎn)錄測序方法[如Smart-seq2 (Picelli et al., 2013)认烁、MATQ-seq (Sheng et al., 2017)和SUPeR-seq (Fan X. et al.肿男, 2015)];和(ii) 3 -端, Drop-seq (Macosko et al. ,2015)却嗡, Seq-Well (Gierahn et al., 2017)舶沛, Chromium (Zheng et al., 2017), and DroNC-seq (Habib et al., 2017)] or 5 -end[如STRT-seq (Islam et al.窗价, 2011, 2012)]轉(zhuǎn)錄測序技術(shù)如庭。每個scRNA-seq協(xié)議都有其優(yōu)缺點,導(dǎo)致不同的scRNA-seq方法具有不同的特點和不同的性能(Ziegenhain et al., 2017)撼港。在進行單細胞轉(zhuǎn)錄組研究時坪它,考慮到研究目標和測序成本之間的平衡,可能需要使用特定的scRNA-seq技術(shù)餐胀。

思考:這些protocol中為什么會區(qū)分3'-end 和5'-end?

由于起始量少哟楷,所以scRNA-seq具有捕獲效率低瘤载、漏檢率高的局限性否灾。與 Bulk RNA-seq相比,scRNA-seq產(chǎn)生的數(shù)據(jù)更嘈雜鸣奔、更多變墨技。技術(shù)噪聲和生物變異(如隨機轉(zhuǎn)錄)對scRNA-seq數(shù)據(jù)的計算分析提出了實質(zhì)性的挑戰(zhàn)。已經(jīng)設(shè)計了各種工具來進行不同的bulk RNA-seq數(shù)據(jù)分析挎狸,但其中許多方法不能直接應(yīng)用于scRNA-seq數(shù)據(jù)扣汪。除了short-read mapping,幾乎所有的數(shù)據(jù)分析(如差異表達锨匆、細胞聚類和基因調(diào)控網(wǎng)絡(luò)推斷)在方法上都存在一定的差異崭别。由于高技術(shù)噪聲冬筒,質(zhì)量控制(QC)是識別和去除低質(zhì)量scRNA-seq數(shù)據(jù)的關(guān)鍵,以獲得可靠的和可重復(fù)的結(jié)果茅主。此外,一些分析包括可變剪接(AS, Alternative Splicing)檢測舞痰、等位基因表達勘探和rna編輯識別不適合scRNA-seq 3'或5'端測序 protocol, 但這些分析方法可以適用于whole-transcript scRNA-seq protocol 產(chǎn)生的數(shù)據(jù)。另一方面,越來越多的工具是專門為scRNA-seq數(shù)據(jù)分析提出诀姚,每種方法都有自己的優(yōu)點和缺點响牛。因此,為了有效地處理高可變性的scRNA-seq數(shù)據(jù),要注意選擇適當?shù)姆治龇椒ā?/p>

本文旨在總結(jié)和討論當前可用的scRNA-seq技術(shù)和各種數(shù)據(jù)分析方法。我們首先介紹了獨特的單細胞分離方法和各種scRNA-seq技術(shù)赫段,這是近年來發(fā)展起來的呀打。然后,我們重點分析了scRNA-seq數(shù)據(jù)糯笙,并強調(diào)了Bulk RNA-seq和scRNA-seq數(shù)據(jù)之間的分析差異贬丛。考慮到scRNA-seq數(shù)據(jù)的高技術(shù)噪聲和復(fù)雜性给涕,我們還對如何選擇合適的工具來分析scRNA-seq數(shù)據(jù)并確保結(jié)果的再現(xiàn)性提出了建議瘫寝。

單細胞分離(Isolation of Single Cells)

盡管對scRNA-seq來說,捕獲效率是一個很大的挑戰(zhàn)稠炬,但是scRNA-seq的第一步是分離單個細胞(圖1)焕阿。目前,分離單個細胞有幾種不同的方法首启,包括有限稀釋法(Limiting dilution technique)暮屡、顯微操作法(micromanipulation)、流式細胞分選(FACS)毅桃、激光捕獲顯微分離(LCM)和微流控(microfluidics)褒纲。有限稀釋技術(shù)(Limiting dilution technique)是利用移液管稀釋分離細胞,這種方法的主要缺點是效率低下钥飞。顯微操作是一種經(jīng)典的方法莺掠,用于從少量細胞樣本中提取細胞,如早期胚胎或未培養(yǎng)的微生物读宙,而這種技術(shù)是耗時和低通量的彻秆。流式細胞儀廣泛用于分離單個細胞,在懸浮狀態(tài)下需要較大的起始體積(>10,000個細胞)结闸。LCM是一種先進的技術(shù)唇兑,利用計算機輔助的激光系統(tǒng)將單個細胞從固體組織中分離出來。微流控技術(shù)以其樣品消耗量低桦锄、流體控制精確扎附、分析成本低等特點,越來越受到人們的重視结耀。這些單細胞分離方案具有各自的優(yōu)點留夜,在捕獲效率和目標細胞純度方面表現(xiàn)出明顯的性能匙铡。

image.jpeg

現(xiàn)有scRNA-seq技術(shù)(Currently Available ScRNA-Seq Technologies)

至今,許多scRNA-seq技術(shù)被提出用于單細胞轉(zhuǎn)錄組研究(表1)碍粥。Tang等人(2009)發(fā)表了第一個scRNA-seq方法慰枕,隨后又開發(fā)了許多其他scRNA-seq方法。這些scRNA-seq技術(shù)至少在以下一個方面不同:(i) 細胞分離; (ii)細胞裂解; (iii) 反轉(zhuǎn)錄;(iv) 擴增;(v) 轉(zhuǎn)錄本覆蓋度;(vi) 鏈特異性; (vii) UMI(獨特的分子標識符即纲,可用于檢測和量化獨特轉(zhuǎn)錄本的分子標記)可用性具帮。這些scRNA-seq方法顯著區(qū)別之一是,他們中的一些人可以產(chǎn)生全長(或幾乎全長)轉(zhuǎn)錄測序數(shù)據(jù)(例如,Smart-seq2 SUPeR-seq, MATQ-seq),而其他人只捕獲和序列的3' 端(如Drop-seq Seq-Well DroNC-seq, SPLiT-seq(羅森博格et al ., 2018))或 5' 端(例如,STRT-seq)記錄(表1),不同的scRNA-seq protocol 可能擁有不同的優(yōu)勢和劣勢,一些發(fā)表的文章已經(jīng)詳細比較了其中的一部分軟件。之前的研究表明低斋,Smart-seq2可以檢測到比其他scRNA-seq技術(shù)更多的表達基因蜂厅,包括cell-seq2、mas-seq膊畴、Smart-seq 和 Drop-seq protocol掘猿。最近,Sheng等人(2017)表明唇跨,另一種全長轉(zhuǎn)錄測序方法MATQ-seq在檢測低豐度基因方面可以優(yōu)于Smart-seq2稠通。

image.jpeg

與3’端或5’端計數(shù)方法相比,全長scRNA-seq方法在亞型使用分析买猖、等位基因表達檢測改橘、RNA編輯鑒定等方面具有無可比擬的優(yōu)勢。此外玉控,對于某些低表達基因/轉(zhuǎn)錄本的檢測飞主,全長scRNA-seq方法可能優(yōu)于3’測序方法。值得注意的是高诺,基于droplet-based的技術(shù)[例如碌识, Drop-seq (Macosko et al., 2015)、InDrop (Klein et al., 2015)和Chromium (Zheng et al., 2017)]通呈可以提供更大的細胞通量筏餐,而且與全轉(zhuǎn)錄scRNA-seq相比,每個細胞的測序成本更低牡拇。因此魁瞪,基于droplet-based的方案適合生成大量細胞來識別復(fù)雜組織或腫瘤樣本的細胞亞群。

值得注意的是诅迷,有幾種scRNA-seq技術(shù)可以同時捕獲polyA陽性和polyA陰性的 RNAs佩番,如SUPeR-seq (Fan X. et al., 2015)和MATQ-seq (Sheng et al., 2017)众旗。這些protocol對于長鏈非編碼rna (lncRNAs)和環(huán)狀rna (circRNAs)的測序非常有用罢杉。大量研究表明,lncRNAs和circRNAs在細胞的多種生物學(xué)過程中發(fā)揮重要作用贡歧,可能是癌癥的重要生物標志物滩租。因此赋秀,這種scRNA-seq方法可以提供前所未有的機會,在單細胞水平上全面探索蛋白質(zhì)編碼和非編碼rna的表達動態(tài)律想。

與傳統(tǒng)的Bulk RNA-seq技術(shù)相比猎莲,scRNA-seq protocol 存在較大的技術(shù)差異。為了評估不同細胞間的技術(shù)差異技即,spikein[如 External RNA Control Consortium (ERCC) controls (External, 2005)]和UMIs被廣泛應(yīng)用于相應(yīng)的scRNA-seq方法中著洼。RNA spike-in 是RNA轉(zhuǎn)錄本(具有已知的序列和數(shù)量),用于校準RNA雜交測定的測量結(jié)果而叼,如RNA-seq身笤,而UMIs在理論上可以用于估算絕對分子數(shù)。值得注意的是葵陵,ERCC和UMIs并不適用于所有的scRNA-seq技術(shù)液荸,因為它們之間存在固有的protocol差異。spikein用于Smart-seq2和SUPeR-seq等方法脱篙,但與基于droplet-based 方法不兼容娇钱,而UMIs通常用于3' 端測序技術(shù)[如Drop-seq (Macosko et al., 2015)、InDrop (Klein et al., 2015)和MARS-seq (Jaitin et al., 2014)]绊困。因此文搂,用戶可以根據(jù)技術(shù)特性和優(yōu)點、需要測序的細胞數(shù)和成本考慮來選擇合適的scRNA-seq方法秤朗。

思考:spike-in和UMI的作用是什么细疚?

單細胞數(shù)據(jù)的比對和質(zhì)控(Read Alignment and Expression Quantification of ScRNA-Seq Data)

reads的比對率是scRNA-seq數(shù)據(jù)整體質(zhì)量的重要指標。由于scRNA-seq和bulk RNA-seq技術(shù)通常都是將轉(zhuǎn)錄本序列轉(zhuǎn)換為read來生成fastq格式的原始數(shù)據(jù)川梅,所以這兩種類型的RNA-seq數(shù)據(jù)在read 比對方面沒有區(qū)別疯兼。最初為bulk RNA-seq開發(fā)的比對工具也適用于scRNA-seq數(shù)據(jù)。為比對RNA-seq數(shù)據(jù)贫途,已經(jīng)設(shè)計了許多比對軟件吧彪,這在之前已經(jīng)廣泛討論過。通常丢早,比對算法主要分為兩類:基于空間種子索引的和基于Burrows-Wheeler transform (BWT)的(Li and Homer, 2010)姨裸。目前流行的比對軟件:TopHat2,STAR,和HISAT 在速度和準確性上都表現(xiàn)很好,并且他們可以有效地數(shù)十億read 比對到參考基因組或轉(zhuǎn)錄組(表2)。STAR是基于suffix-array方法怨酝,比TopHat2更快,但它比對時需要一個巨大的內(nèi)存大小(28 g用于人類基因組)傀缩。Engstrom等人系統(tǒng)地評估了26種read比對protocol(不包括HISAT),發(fā)現(xiàn)不同的比對工具有不同的優(yōu)缺點农猬,有些程序的比對速度更快赡艰,但拼接連接檢測的準確性更低。HISAT是基于BWT和Ferragina-Manzini (FM)方法開發(fā)的斤葱。Kim等人(2015)指出慷垮,HISAT是目前最快的工具揖闸,可以達到與其他可用對準器相同或更好的精度。

image.jpeg

數(shù)據(jù)質(zhì)控(Quality Control of ScRNA-Seq Data)

scRNA-seq的局限性包括轉(zhuǎn)錄本覆蓋的偏差料身、低捕獲效率和測序覆蓋度低汤纸,導(dǎo)致scRNA-seq數(shù)據(jù)的技術(shù)噪聲水平高于bulk RNA-seq數(shù)據(jù)。即使對一些很靈敏的scRNA-seq protocol芹血,

也是個很常見的現(xiàn)象贮泞,一些特定的轉(zhuǎn)錄本無法被檢測到,這個術(shù)語叫dropout events幔烛。通常隙畜,scRNA-seq實驗可以從破損、死亡或與多個細胞混合的細胞中生成部分低質(zhì)量的數(shù)據(jù)说贝。

這些低質(zhì)量的細胞將阻礙下游的分析议惰,并可能導(dǎo)致數(shù)據(jù)的誤讀。因此乡恕,對scRNA-seq數(shù)據(jù)的QC質(zhì)量控制是識別和去除低質(zhì)量細胞的關(guān)鍵言询。

為了排除scRNA-seq中的低質(zhì)量細胞,在細胞捕獲步驟中應(yīng)密切注意避免多細胞或死亡細胞傲宜。測序后运杭,需要進行一系列的QC分析來剔除來自低質(zhì)量細胞的數(shù)據(jù)。由于測序深度不足可能導(dǎo)致大量低表達和中度表達的基因丟失函卒,因此這些樣本只包含少量的reads辆憔,應(yīng)該首先丟棄。然后报嵌,可以使用最初為bulk RNA-seq數(shù)據(jù)QC開發(fā)的工具(如FastQC)來檢查scRNA-seq數(shù)據(jù)的測序質(zhì)量虱咧。此外,在讀取比對之后锚国,應(yīng)該刪除比對率非常低的樣本腕巡,因為它們包含大量未必對上的read,這可能是由于RNA降解造成的血筑。如果在scRNA-seq中使用外源性噪聲(如ERCC)绘沉,則可以對技術(shù)噪聲進行估計。具有極高比例的比對率比對到spike-in的細胞表明豺总,它們可能在單元格捕獲過程中被破壞车伞,應(yīng)該被刪除。細胞質(zhì)rna通常會丟失喻喳,但線粒體rna會保留在受損細胞中另玖,因此,繪制到線粒體基因組的reads比例也有助于識別低質(zhì)量細胞(Bacher和Kendziorski, 2016)。此外日矫,在每個細胞中可以檢測到的表達基因/轉(zhuǎn)錄的數(shù)量也具有啟發(fā)性赂弓。如果在一個細胞中只能檢測到少量的基因绑榴,那么這個細胞很可能已經(jīng)受損或死亡哪轿,或者是RNA降解了∠柙酰考慮到scRNA-seq數(shù)據(jù)的高噪聲窃诉,通常采用1 FPKM/RPKM的閾值來定義表達基因。提出了一些scRNA-seq的QC方法包括SinQC 和Scater 赤套,這些工具對scRNA-seq數(shù)據(jù)的質(zhì)量控制是有用的飘痛。

批次效應(yīng)校正(Batch Effect Correction)

批次效應(yīng)是高通量測序?qū)嶒炛谐R姷募夹g(shù)變異來源。scRNA-seq的創(chuàng)新和成本的降低使得許多研究能夠分析大量細胞的轉(zhuǎn)錄組容握。

大規(guī)模scRNA-seq數(shù)據(jù)集可以在不同的時間由不同的操作人員單獨生成宣脉,也可以在多個實驗室使用不同的細胞分離protocol、文庫制備方法和或測序平臺生成剔氏。

這些因素會引入系統(tǒng)錯誤塑猖,混淆技術(shù)和生物變異,導(dǎo)致一個批次的基因表達譜與另一個批次的基因表達譜存在系統(tǒng)差異(Leek et al., 2010;(Hicks et al., 2018)谈跛。

因此羊苟,批次效應(yīng)是scRNA-seq數(shù)據(jù)分析的一個主要挑戰(zhàn),它可能掩蓋了潛在的生物學(xué)特性感憾,并導(dǎo)致虛假的結(jié)果蜡励。

為了避免不正確的數(shù)據(jù)整合和解釋,批次效應(yīng)必須在下游分析之前糾正阻桅。

由于scRNA-seq與bulk RNA-seq在數(shù)據(jù)特征上的差異凉倚,本文提出了針對bulk RNA-seq的批次校正方法[例如, RUVseq (Risso et al., 2014)和svaseq (Leek, 2014)]可能不適合scRNA-seq嫂沉。

最近設(shè)計了幾種方法來緩解scRNA-seq數(shù)據(jù)中的批處理效應(yīng)占遥,如MNN(相互最近鄰居)(Haghverdi et al., 2018)和kBET (k-最近鄰居批處理效應(yīng)測試)(Buttner et al., 2019)。

MNN使用來自不同批次中最相似細胞的數(shù)據(jù)校正批次效果输瓜。KBET χ2-based方法為量化scRNA-seq數(shù)據(jù)批處理的影響瓦胎。

這些針對scRNA-seq數(shù)據(jù)的特定批處理校正方法可以比針對批次RNA-seq開發(fā)的方法表現(xiàn)得更好(Haghverdi et al., 2018;Buttner et al., 2019)。

單細胞數(shù)據(jù)標準化(Normalization of ScRNA-Seq Data)

為了正確地解釋scRNA-seq數(shù)據(jù)的結(jié)果尤揣,標準化是通過調(diào)整由捕獲效率搔啊、測序深度、dropouts和其他技術(shù)影響引起的不想要的偏差來獲得感興趣信號的重要步驟北戏。由于起始量低(starting material)负芋、實驗方案具有挑戰(zhàn)性,scRNA-seq的技術(shù)噪聲是一個明顯的問題。scRNA-seq數(shù)據(jù)的標準化將有利于后續(xù)分析旧蛾,包括細胞亞群識別和差異表達calling莽龟。一般來說,normalization可以分為兩種不同的類型:樣品內(nèi)標準化和樣品間標準化锨天。樣品內(nèi)標準化只要是要消除基因特異性(gene-specific)偏差,這使得基因表達在一個樣本內(nèi)具有可比性(如RPKM/FPKM和TPM)毯盈。而樣本間歸一化則是通過調(diào)整樣本特異性差異(如測序深度和捕獲效率)來實現(xiàn)樣本間基因表達的比較。通常病袄,這些簡單的標準化策略是基于排序深度或上四分位數(shù)搂赋。如果在scRNA-seq protocol中使用了spike-in或UMI,則可以根據(jù)spike-in /UMIs的perfomance 來進行標準化益缠。

bulk RNA-seq 數(shù)據(jù)分析中有很多處理樣品間標準化的方法脑奠。比如DESeq2, TMM(trimmed mean of M values)。DEseq2根據(jù)不同樣本的read count 來計算比例因子(scaling factor)幅慌,而TMM則是刪除極端的對數(shù)倍變化(removes the extreme log fold changes)宋欺。但是這種bulk-based標準化方法不能用于scRNA-seq數(shù)據(jù),因為scRNA-seq產(chǎn)生大量的零表達值胰伍,并且比bulk RNA-seq具有更高的技術(shù)變異水平齿诞。因此使用bulk RNA-seq 標準化方法可能會導(dǎo)致對低表達基因的scRNA-seq過度校正。當然喇辽,有很多處理scRNA-seq標準化的方法掌挚,比如SCnorm,SAMstrt 和最近引入的一種反卷積方法菩咨,它使用跨細胞池的求和表達式值來進行標準化吠式。SCnorm是基于分位數(shù)回歸的方法。SAMstrt是基于spike-in抽米。Bacher等(2017)認為特占,如果用基于bulk RNA-seq開發(fā)的傳統(tǒng)歸一化方法去對scRNA-seq數(shù)據(jù)做標準化可能會引入artifacts,而SCnorm可以有效地對scRNA-seq數(shù)據(jù)進行歸一化云茸,改善主成分分析(PCA)和差異表達基因的識別是目。

scRNA-seq數(shù)據(jù)的imputation(Imputation of ScRNA-Seq Data)

單細胞RNA-seq數(shù)據(jù)通常包含許多由于原始RNA擴增失敗而導(dǎo)致的缺失值或dropouts。(注解: dropouts是指空beats)标捺。dropout 的比例是跟 protocol 相關(guān)懊纳,而且與每個細胞測序的reads數(shù)密切相關(guān)。dropout event 增加了細胞間的變異亡容。導(dǎo)致信號對每個基因的影響嗤疯,模糊了基因-基因關(guān)系的檢測。由于在scRNA-seq中可能無法檢測到大量真正表達的轉(zhuǎn)錄本闺兢,因此dropout 可能會很大程度上影響下游的分析茂缚。Imputation是一種用替代值替換缺失數(shù)據(jù)(dropouts)的有效方法。雖然已經(jīng)提出了一些方法來估算bulk RNA-seq數(shù)據(jù),但它們并不直接適用于scRNA-seq數(shù)據(jù)脚囊。最近為scRNA-seq開發(fā)了幾種imputation方法龟糕。包括SAVER,MAGIC悔耘,ScImpute讲岁,DrImpute,AutoImpute淮逊。SAVER是一個基于貝葉斯的模型催首,用于基于umi的scRNA-seq數(shù)據(jù)來恢復(fù)所有基因的真實表達水平扶踊。MAGIC通過構(gòu)建基于馬爾可夫親和度的基因表達圖來進行基因表達的計算泄鹏。ScImpute的開發(fā)者認為,SAVER和MAGIC可能會導(dǎo)致未受dropout影響的基因的表達變化秧耗,而ScImpute可以利用其他類似細胞中不太可能受dropout影響的相同基因的信息备籽,在不引入新的偏差的情況下計算dropout值。DrImpute是一種基于集群的方法分井,可以有效地將dropout中的零從真正的零中分離出來车猬。AutoImpute是一種基于自編碼的方法,它通過學(xué)習(xí)scRNA-seq數(shù)據(jù)的固有分布來尋找(impute)缺失的值尺锚。最近珠闰,Zhang等人評估了不同的imputation方法,發(fā)現(xiàn)這些方法的性能與其模型假設(shè)和可擴展性相關(guān)瘫辩。

降維和特征選擇(Dimensionality Reduction and Feature Selection)

單細胞RNA測序數(shù)據(jù)具有高維性伏嗜,可能涉及數(shù)千個基因和大量細胞。降維和特征選擇是處理高維數(shù)據(jù)的兩種主要策略伐厌。降維方法通常通過最優(yōu)地保留原始數(shù)據(jù)的一些關(guān)鍵屬性來將數(shù)據(jù)映射到更低的維度空間中承绸。PCA是一種線性降維算法,它假設(shè)數(shù)據(jù)近似正態(tài)分布挣轨。T-SNE(T-distributed stochastic neighbor embedding)是一種主要用于高維數(shù)據(jù)可視化的非線性方法军熏。PCA和t-SNE在大量的的scRNA-seq研究中被廣泛應(yīng)用,以減少數(shù)據(jù)維度卷扮,并將細胞區(qū)分為不同的亞群荡澎。值得注意的是,PCA不能有效地表示scRNA-seq數(shù)據(jù)的復(fù)雜結(jié)構(gòu)晤锹,而t-SNE具有計算時間慢摩幔、多次處理同一數(shù)據(jù)集的嵌入方式不同的局限性。最近抖甘,為了減少scRNA-seq數(shù)據(jù)的維度热鞍,專門開發(fā)了UMAP (uniform manifold approximation and projection) (Becht et al., 2018)和scvis (Ding et al., 2018)。Becht等人指出,與其他降維方法相比薇宠,UMAP提供了最快的運行時間偷办、最高的再現(xiàn)性和最有意義的細胞簇組織(Becht et al., 2018)。

特征選擇刪除了沒有信息的基因澄港,并識別最相關(guān)的特征椒涯,以減少下游分析中使用的維數(shù)。通過執(zhí)行特征選擇來減少基因數(shù)量回梧,可以在很大程度上加速大規(guī)模scRNA-seq數(shù)據(jù)的計算废岂。差異表達在bulk RNA-seq實驗中是一種廣泛應(yīng)用的特征選擇方法,但由于scRNA-seq數(shù)據(jù)差異表達需要預(yù)定的信息或同構(gòu)的亞群信息狱意,因此很難應(yīng)用于scRNA-seq數(shù)據(jù)中湖苞。針對scRNA-seq數(shù)據(jù)設(shè)計的無監(jiān)督特征選擇算法可分為以下幾組:(i)基于高度可變基因(HVG)的;(2)基于spike-in;和(iii)基于dropout。HVG方法依賴于一種假設(shè)详囤,即細胞間表達高度可變的基因是由生物學(xué)效應(yīng)而非技術(shù)噪音造成的财骨。HVG方法包括Brennecke等人(2013)提出的算法,以及Seurat實現(xiàn)的FindVariableGenes (FVG) (Satija等人藏姐,2015)隆箩。基于spikein的方法可以識別出那些表現(xiàn)出顯著高于那些表達水平相似的spikein的基因羔杨。[如scLVM (Buettner et al.捌臊, 2015)和BASiCS (Vallejos et al., 2015)]兜材,兩者對HVG的概念相似理澎。基于Dropout的方法利用scRNA-seq數(shù)據(jù)的Dropout分布進行特征選擇护姆,如M3Drop (Andrews and Hemberg, 2018b)矾端。Andrews和Hemberg證明了他們的M3Drop工具優(yōu)于現(xiàn)有的基于方差(variance-based)的特征選擇方法。

細胞亞群的鑒定(Cell Subpopulation Identification)

scRNA-seq數(shù)據(jù)分析的一個關(guān)鍵目標是鑒定特定條件或組織內(nèi)的細胞亞群(不同的細胞群通常是不同的細胞類型)卵皂,以揭示細胞的異質(zhì)性秩铆。細胞亞群鑒定應(yīng)在質(zhì)量控制和scRNA-seq數(shù)據(jù)標準化后進行,否則可能會引入人為因素灯变∨孤辏基于先驗信息的聚類方法主要可以分為兩類。如果使用一組已知marker進行聚類添祸,則這些方法是基于先驗信息的滚粟。另外,無監(jiān)督聚類方法也可用于用scRNA-seq數(shù)據(jù)重新識別細胞群刃泌。無監(jiān)督聚類算法主要分為以下幾種類型: (i) k-means;(ii) 層次聚類(hierarchical clustering); (iii) density-based clustering; and (iv) graph-based clustering凡壤。K-means是一種快速的方法署尤,它將細胞分配到最近的分類中心,并且它需要預(yù)先確定cluster數(shù)量亚侠。層次聚類可以確定聚類之間的關(guān)系曹体,但是它的工作速度通常比k-means慢∠趵茫基于密度(density-based)的聚類方法需要大量的樣本來精確計算密度箕别,通常假設(shè)所有的聚類具有相同的密度。graph-based聚類可以看作是基于密度的聚類的一種擴展滞谢,可以應(yīng)用于數(shù)百萬個細胞串稀。一些為scRNA-seq數(shù)據(jù)特別設(shè)計的聚類方法,如單細胞一致性聚類(single cell consensus clustering, SC3)狮杨。和Seurat中的聚類方法(Satija et al ., 2015),它可以更容易的鑒定細胞亞群(表3)母截。SC3是一種無監(jiān)督的方法,結(jié)合多個聚類方法,在單細胞cluster中具有精度高和穩(wěn)健性。Seurat主要基于共享近鄰(shared nearest neighbor, SNN)聚類算法來識別細胞簇禾酱。一旦確定了亞群微酬,通常通過差異表達調(diào)用或方差分析(ANOVA)來識別最能區(qū)分不同亞群的marker绘趋。

image.jpeg

ScRNA-Seq數(shù)據(jù)的差異表達分析(Differential Expression Analysis of ScRNA-Seq Data)

差異表達分析對于發(fā)現(xiàn)不同亞群或細胞群間的差異表達基因(DEGs)非常有用颤陶。DEGs對于解釋兩個比較條件之間的生物學(xué)差異至關(guān)重要。scRNA-seq數(shù)據(jù)的技術(shù)可變性陷遮、高噪聲(例如滓走,dropout)和大樣本規(guī)模給差異表達calling帶來了挑戰(zhàn)。此外,在一個細胞群中可能存在多種細胞狀態(tài),導(dǎo)致細胞中基因表達的多樣性巡蘸。最初為 Bulk RNA-seq數(shù)據(jù)開發(fā)的工具已經(jīng)在許多單細胞研究中用于識別DEGs蜗搔,但是這些方法對scRNA-seq數(shù)據(jù)的適用性仍然不清楚。近年來福澡,提出了一些基于scRNA-seq數(shù)據(jù)進行差異表達基因calling的具體方法,如MAST (Finak et al., 2015)、SCDE (Kharchenko et al., 2014)涛漂、DEsingle (Miao et al., 2018)、Census (Qiu et al., 2017)和BCseq (Chen and Zheng,2018)(表4)检诗。MAST是基于線性模型擬合和似然比檢驗匈仗。SCDE是一種貝葉斯方法,使用低強度泊松過程來解釋dropouts逢慌。DEsingle使用了零膨脹負二項模型來估計損失和實際零悠轩。BCseq以數(shù)據(jù)自適應(yīng)的方式降低技術(shù)噪音。Soneson和Robinson最近評估了36種差異表達方法(包括為scRNA-seq設(shè)計的工具和bulk RNA-seq數(shù)據(jù))攻泼,發(fā)現(xiàn)這些方法在DEGs的特征和數(shù)量上存在顯著差異(Soneson and Robinson, 2018)火架。將開發(fā)越來越多的用于scRNA-seq數(shù)據(jù)差異表達分析的工具鉴象,并鼓勵用戶根據(jù)scRNA-seq數(shù)據(jù)的復(fù)雜特性選擇專門為scRNA-seq設(shè)計的工具來識別DEGs。

image.jpeg

細胞譜系和擬時間重構(gòu)(Cell Lineage and Pseudotime Reconstruction)

許多生物系統(tǒng)中的細胞表現(xiàn)出連續(xù)的狀態(tài)譜何鸡,并涉及不同細胞狀態(tài)之間的轉(zhuǎn)換炼列。通過基于scRNA-seq數(shù)據(jù)重建細胞軌跡和偽時間,可以對部分細胞內(nèi)的動態(tài)過程進行計算建模音比。擬時間 (Pseudotime) 是一個系統(tǒng)中沿著一個連續(xù)發(fā)展過程軌跡的細胞排序俭尖,它允許在軌跡的開始、中間和結(jié)束狀態(tài)識別細胞類型(Griffiths et al., 2018)洞翩。除了揭示細胞間的基因表達動態(tài)外稽犁,單細胞軌跡推斷還有助于識別觸發(fā)狀態(tài)轉(zhuǎn)換的因素。一些工具已經(jīng)提出了軌跡推理,例如,Monocle(Trapnell et al., 2014),Waterfall(Shin et al.,2015),Wishbone(Setty et al., 2016), TSCAN (Ji and Ji, 2016),Monocle2(Qiu et al ., 2017),Slingshot(Street et al ., 2018),和CellRouter (Lummertz da Rocha et al., 2018)(表5)骚亿。由此產(chǎn)生的軌跡拓撲可以是線性的,分叉或一棵樹/圖結(jié)構(gòu)已亥。Monocle構(gòu)建了一個最小生成樹(MST),用于細胞基于獨立成分分析(ICA)搜索最長的主鏈来屠。Monocle2使用了一種獨特的方法虑椎,它結(jié)合了非監(jiān)督數(shù)據(jù)驅(qū)動方法和反向圖嵌入(RGE),這比Monocle更健壯俱笛、更快捆姜。Slingshot是一種基于cluster的方法,用于識別具有不同監(jiān)督級別的多個軌跡迎膜。CellRouter利用flow network 來識別細胞狀態(tài)轉(zhuǎn)換軌跡泥技。最近,Saelens等人(2018)評估了一些單細胞軌跡推斷方法(不包括CellRouter)磕仅,發(fā)現(xiàn)Slingshot珊豹、TSCAN和Monocle2優(yōu)于其他方法。

image.jpeg

可變剪切和RNA編輯(Alternative Splicing and RNA Editing Analysis of ScRNA-Seq Data)

大多數(shù)已發(fā)表的單細胞研究主要探討基因水平上單個細胞之間的轉(zhuǎn)錄組差異榕订。在真核生物的基因組中店茶,可變剪切(AS)允許多外顯子基因產(chǎn)生不同的亞型,這在很大程度上增加了蛋白質(zhì)編碼和非編碼rna的多樣性〗俸悖現(xiàn)有公認的5種基本模式:exon-skipping (cassette exon), mutually exclusive exons, alternative donor site, alternative acceptor site, and intron retention.大量研究表明贩幻,AS在哺乳動物中很常見,基于bulk RNA-seq數(shù)據(jù)兼贸,90%以上的人類基因可以發(fā)生AS段直。此外,AS在多種生物學(xué)過程中起著關(guān)鍵作用溶诞,并且AS異逞烀剩可能與癌癥相關(guān)(Sveen et al., 2016)。Bulk RNA-seq 數(shù)據(jù)揭示的結(jié)果只能反映在種群水平上眾多細胞的平均模式螺垢。由于scRNA-seq數(shù)據(jù)的高噪音(如dropouts和不均勻的轉(zhuǎn)錄覆蓋)和低序列覆蓋喧务,最初開發(fā)的用于Bulk RNA-seq數(shù)據(jù)的拼接量化方法不適用于scRNA-seq數(shù)據(jù)赖歌。由于表達動態(tài)是細胞群體的一個關(guān)鍵方面,因此在單細胞分辨率方面進行研究以了解細胞水平亞型的使用是很有希望的功茴。迄今為止庐冯,只有少數(shù)的AS檢測方法被設(shè)計用于scRNA-seq數(shù)據(jù),包括SingleSplice (Welch et al., 2016), Census (Qiu et al., 2017), BRIE (Huang and Sanguinetti, 2017), and Expedition (Song et al., 2017) (Table 6).SingleSplice使用統(tǒng)計模型來檢測具有顯著亞型使用的基因坎穿,而不需要估計全長轉(zhuǎn)錄本的表達水平展父。Census以Dirichlet-multinomial分布的線性模式來模擬每個基因的亞型計數(shù)。BRIE是一個用于微分亞型量化的貝葉斯層次模型玲昧。Expedition包含一套用于識別AS栖茉、分配拼接模式和可視化模式更改的算法。

image.jpeg

另一方面孵延,RNA編輯是一個重要的轉(zhuǎn)錄后處理事件吕漂,會導(dǎo)致RNA分子的序列改變(Gott and Emeson, 2000)。類似地尘应,rna編輯主要使用Bulk RNA-seq技術(shù)進行研究惶凝,但很少在單細胞水平進行研究。目前犬钢,scRNA-seq的局限性很大程度上阻礙了rna編輯檢測在單個細胞中的應(yīng)用苍鲜。因此,隨著scRNA-seq技術(shù)和單細胞編輯檢測算法的發(fā)展娜饵,探索單細胞間的rna編輯動態(tài)將成為可能坡贺。值得注意的是,AS和rna編輯主要適用于scRNA-seq protocol生成的數(shù)據(jù)箱舞,這些protocol可以對Smart-seq2和MATQ-seq等測全長轉(zhuǎn)錄本文庫,而不是3'-end scRNA-seq方法拳亿。

等位基因表達探索(Allelic Expression Exploration with ScRNA-Seq Data)

二倍體物種包含兩套染色體晴股,分別來自它們的父母。等位基因表達分析可以揭示基因在親本和母本基因組中的表達是否相等肺魁。對于常染色體电湘,親代和母體的表達通常是平等的,親代或母體基因組的異常表達可能導(dǎo)致某些疾病(McKean et al., 2016)鹅经。迄今為止寂呛,基于scRNA-seq數(shù)據(jù)的全基因組等位基因表達譜檢測方法較少。等位基因表達calling需要注意scRNA-seq數(shù)據(jù)的高dropouts可能會引入許多假陽性瘾晃。Deng等(2014)在研究小鼠植入前胚胎的等位基因表達譜時贷痪,使用了一系列嚴格的標準來過濾由于scRNA-seq的技術(shù)可變性而可能產(chǎn)生的假等位基因。該策略的穩(wěn)健性在利用小鼠胚胎干細胞分析X染色體失活隨發(fā)育進程的動力學(xué)過程中得到進一步證明(Chen et al.蹦误, 2016a)劫拢。SCALE采用經(jīng)驗貝葉斯方法將基因表達分為沉默狀態(tài)肉津、單等位基因狀態(tài)和雙等位基因狀態(tài)(Jiang et al., 2017)舱沧。我們認為妹沙,單細胞水平的等位基因表達分析可以在很大程度上促進對劑量補償和相關(guān)疾病的潛在機制的理解。值得注意的是熟吏,單細胞水平的等位基因表達研究也需要whole-transcript scRNA-seq距糖,主要適用于有父系和母系單核苷酸多態(tài)性(SNP)信息的生物體。

基因調(diào)控網(wǎng)絡(luò)重構(gòu)(Gene Regulatory Network Reconstruction)

基因調(diào)控網(wǎng)絡(luò)推理(Gene regulatory network inference)已經(jīng)在眾多的Bulk RNA-seq研究中得到廣泛應(yīng)用牵寺,而scRNA-seq也為此類分析提供了巨大的潛力肾筐。對于大量的RNA-seq數(shù)據(jù),通常使用加權(quán)基因共表達網(wǎng)絡(luò)分析(WGCNA)等工具從大量樣本中構(gòu)建網(wǎng)絡(luò)缸剪。一個基本的假設(shè)是吗铐,表達高度相關(guān)的基因可以被共同調(diào)控。由于這樣的分析無法確定監(jiān)管關(guān)系杏节,因此產(chǎn)生的網(wǎng)絡(luò)通常是無定向的唬渗。從理論上講,scRNA-seq細胞可以作為Bulk RNA-seq的樣本奋渔,類似的方法也可以應(yīng)用于scRNA-seq數(shù)據(jù)構(gòu)建基因調(diào)控網(wǎng)絡(luò)镊逝。

對scRNA-seq數(shù)據(jù)的網(wǎng)絡(luò)推斷可以揭示有意義的基因相關(guān)性,并提供bulk RNA-seq的群體水平數(shù)據(jù)所不能揭示的生物學(xué)上的重要見解嫉鲸。然而撑蒜,由于scRNA-seq的技術(shù)噪聲和細胞的不同亞群或狀態(tài),網(wǎng)絡(luò)重構(gòu)應(yīng)該引起大家重視玄渗。為了減少虛假結(jié)果座菠,應(yīng)該對每個亞群或同一階段的細胞進行網(wǎng)絡(luò)推理。最近藤树,Aibar等人(2017)開發(fā)了基于scRNA-seq數(shù)據(jù)重建基因調(diào)控網(wǎng)絡(luò)的SCENIC方法浴滴,結(jié)果表明,該方法能夠較好地預(yù)測轉(zhuǎn)錄因子與靶基因之間的相互作用岁钓。PIDC是另一個利用多元信息論從單細胞數(shù)據(jù)推斷基因調(diào)控網(wǎng)絡(luò)的軟件(Chan et al.升略, 2017)。這些網(wǎng)絡(luò)推理工具有助于從單細胞轉(zhuǎn)錄組數(shù)據(jù)中識別表達調(diào)控網(wǎng)絡(luò)屡限,并為基因間的調(diào)控關(guān)系提供了批判性的生物學(xué)見解品嚣。

Conclusion

在過去的10年中,scRNA-seq取得了很大的進步钧大,開發(fā)了多種scRNA-seq protocol翰撑。scRNA-seq的開發(fā)和創(chuàng)新在很大程度上促進了單細胞轉(zhuǎn)錄組學(xué)的研究,并在細胞表達可變性和動力學(xué)方面產(chǎn)生了深刻的發(fā)現(xiàn)拓型。此外额嘿,隨著細胞條形碼技術(shù)和微流體技術(shù)的發(fā)展瘸恼,scRNA-seq的通量也顯著提高。同時册养,最近也提出了可用于固定和冷凍標本的scRNA-seq方法东帅,這將極大地有利于高異質(zhì)性臨床標本的研究。然而球拦,目前可用的scRNA-seq方法仍然存在高退出問題靠闭,即弱表達基因可能被遺漏。RNA捕獲效率和轉(zhuǎn)錄本覆蓋率的提高必將降低scRNA-seq的技術(shù)噪聲坎炼。此外愧膀,由于目前大多數(shù)scRNA-seq方法主要捕獲polyA rna,因此谣光,開發(fā)能夠同時捕獲polyA和polyA- rna的protocol(如MATQ-seq)將能夠在單細胞分辨率上全面研究蛋白質(zhì)編碼和非編碼基因表達動態(tài)檩淋。

由于scRNA-seq數(shù)據(jù)的噪聲較大,因此采用合適的方法來克服對scRNA-seq數(shù)據(jù)的分析是至關(guān)重要的萄金。質(zhì)量控制是必要的排除那些低質(zhì)量的細胞蟀悦,以避免在數(shù)據(jù)解釋中涉及偽影。此外氧敢,批量效應(yīng)校正(如果需要)日戈、樣本歸一化和歸一化之間也很重要,應(yīng)該在細胞亞群識別孙乖、差異表達調(diào)用和其他下游分析之前進行浙炼。此外,細胞大小和細胞周期狀態(tài)等因素可能在某些類型細胞的細胞變異性中發(fā)揮重要作用唯袄,這些偏差也需要考慮弯屈。雖然越來越多的方法被專門設(shè)計來解釋scRNA-seq數(shù)據(jù),但仍需要先進的新方法來有效地處理技術(shù)噪音和細胞的表達變異性越妈。具體來說季俩,利用scRNA-seq數(shù)據(jù)準確分析AS和rna編輯的方法對闡明單個細胞的轉(zhuǎn)錄后機制非常有用。綜上所述梅掠,scRNA-seq數(shù)據(jù)的生物信息學(xué)分析仍然具有挑戰(zhàn)性,在數(shù)據(jù)解釋方面需要特別重視店归,急需更高效的工具阎抒。

總的來說,scRNA-seq及其相關(guān)的計算方法極大地促進了單細胞轉(zhuǎn)錄組學(xué)的發(fā)展消痛。scRNA-seq技術(shù)的不斷創(chuàng)新和生物信息學(xué)方法的不斷進步且叁,將極大地促進生物學(xué)和臨床研究的發(fā)展,為深入了解細胞的基因表達異質(zhì)性和動態(tài)機制提供重要的理論依據(jù)秩伞。

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
  • 序言:七十年代末逞带,一起剝皮案震驚了整個濱河市欺矫,隨后出現(xiàn)的幾起案子,更是在濱河造成了極大的恐慌展氓,老刑警劉巖穆趴,帶你破解...
    沈念sama閱讀 206,013評論 6 481
  • 序言:濱河連續(xù)發(fā)生了三起死亡事件,死亡現(xiàn)場離奇詭異遇汞,居然都是意外死亡未妹,警方通過查閱死者的電腦和手機,發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 88,205評論 2 382
  • 文/潘曉璐 我一進店門空入,熙熙樓的掌柜王于貴愁眉苦臉地迎上來络它,“玉大人,你說我怎么就攤上這事歪赢』粒” “怎么了?”我有些...
    開封第一講書人閱讀 152,370評論 0 342
  • 文/不壞的土叔 我叫張陵埋凯,是天一觀的道長点楼。 經(jīng)常有香客問我,道長递鹉,這世上最難降的妖魔是什么盟步? 我笑而不...
    開封第一講書人閱讀 55,168評論 1 278
  • 正文 為了忘掉前任,我火速辦了婚禮躏结,結(jié)果婚禮上却盘,老公的妹妹穿的比我還像新娘。我一直安慰自己媳拴,他們只是感情好黄橘,可當我...
    茶點故事閱讀 64,153評論 5 371
  • 文/花漫 我一把揭開白布。 她就那樣靜靜地躺著屈溉,像睡著了一般塞关。 火紅的嫁衣襯著肌膚如雪。 梳的紋絲不亂的頭發(fā)上子巾,一...
    開封第一講書人閱讀 48,954評論 1 283
  • 那天帆赢,我揣著相機與錄音,去河邊找鬼线梗。 笑死椰于,一個胖子當著我的面吹牛,可吹牛的內(nèi)容都是我干的仪搔。 我是一名探鬼主播瘾婿,決...
    沈念sama閱讀 38,271評論 3 399
  • 文/蒼蘭香墨 我猛地睜開眼,長吁一口氣:“原來是場噩夢啊……” “哼!你這毒婦竟也來了偏陪?” 一聲冷哼從身側(cè)響起抢呆,我...
    開封第一講書人閱讀 36,916評論 0 259
  • 序言:老撾萬榮一對情侶失蹤,失蹤者是張志新(化名)和其女友劉穎笛谦,沒想到半個月后抱虐,有當?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體,經(jīng)...
    沈念sama閱讀 43,382評論 1 300
  • 正文 獨居荒郊野嶺守林人離奇死亡揪罕,尸身上長有42處帶血的膿包…… 初始之章·張勛 以下內(nèi)容為張勛視角 年9月15日...
    茶點故事閱讀 35,877評論 2 323
  • 正文 我和宋清朗相戀三年梯码,在試婚紗的時候發(fā)現(xiàn)自己被綠了。 大學(xué)時的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片好啰。...
    茶點故事閱讀 37,989評論 1 333
  • 序言:一個原本活蹦亂跳的男人離奇死亡轩娶,死狀恐怖,靈堂內(nèi)的尸體忽然破棺而出框往,到底是詐尸還是另有隱情鳄抒,我是刑警寧澤,帶...
    沈念sama閱讀 33,624評論 4 322
  • 正文 年R本政府宣布椰弊,位于F島的核電站许溅,受9級特大地震影響,放射性物質(zhì)發(fā)生泄漏秉版。R本人自食惡果不足惜贤重,卻給世界環(huán)境...
    茶點故事閱讀 39,209評論 3 307
  • 文/蒙蒙 一、第九天 我趴在偏房一處隱蔽的房頂上張望清焕。 院中可真熱鬧并蝗,春花似錦、人聲如沸秸妥。這莊子的主人今日做“春日...
    開封第一講書人閱讀 30,199評論 0 19
  • 文/蒼蘭香墨 我抬頭看了看天上的太陽粥惧。三九已至键畴,卻和暖如春,著一層夾襖步出監(jiān)牢的瞬間突雪,已是汗流浹背起惕。 一陣腳步聲響...
    開封第一講書人閱讀 31,418評論 1 260
  • 我被黑心中介騙來泰國打工, 沒想到剛下飛機就差點兒被人妖公主榨干…… 1. 我叫王不留咏删,地道東北人疤祭。 一個月前我還...
    沈念sama閱讀 45,401評論 2 352
  • 正文 我出身青樓,卻偏偏與公主長得像饵婆,于是被迫代替她去往敵國和親。 傳聞我的和親對象是個殘疾皇子,可洞房花燭夜當晚...
    茶點故事閱讀 42,700評論 2 345

推薦閱讀更多精彩內(nèi)容