GSEA背后的統(tǒng)計(jì)學(xué)原理

GSEA簡介

基因集富集分析(GSEA)是一種計(jì)算方法,可確定先驗(yàn)定義的基因集是否在統(tǒng)計(jì)上顯示兩種生物學(xué)狀態(tài)之間存在顯著一致的差異琢唾。
來自官網(wǎng):GSEA

GSEA首先由Mootha等人描述。 (Mootha 2003)試圖闡明2型糖尿病的機(jī)理滞伟。 他們認(rèn)為與疾病相關(guān)的基因表達(dá)改變可以在生物學(xué)途徑或共同調(diào)控的基因組而非個體基因的水平上體現(xiàn)出來媚值。 缺乏在基因水平上檢測真實(shí)信號的能力可能是高通量表達(dá)測量的結(jié)果宴树,該測量涉及異類樣品红省,適度的樣品量以及微妙但有意義的改變表達(dá)水平额各。 最終,這些混淆的嘗試得出了具有感興趣的生物學(xué)狀態(tài)的可再現(xiàn)的關(guān)聯(lián)吧恃。

在他們的研究中虾啦,使用微陣列比較了糖尿病(DM2)和正常葡萄糖耐量(NGT)對照小鼠中的基因表達(dá)痕寓。 一方面傲醉,兩種統(tǒng)計(jì)分析方法未能在這兩組小鼠之間鑒定單個差異表達(dá)的基因。 另一方面厂抽,GSEA將氧化磷酸化(OXPHOS)確定為在DM2中被下調(diào)的最高得分基因需频。 接下來的四個得分最高的基因組與OXPHOS很大程度上重疊。

所以GSEA的輸入是一個基因表達(dá)量矩陣筷凤,其中的樣本分成了A和B兩組,首先對所有基因進(jìn)行排序苞七,簡單理解就是根據(jù)處理后的差異倍數(shù)值進(jìn)行從大到小排序

1.Local statistic

在這一步中藐守,我們描述了一種局部或基因水平的測量方法,該方法用于按照GSEA對基因進(jìn)行排名蹂风。 之前卢厂,我們描述了如何獲取和處理RNA-seq數(shù)據(jù)集,使其成為單個基因列表惠啄,該列表根據(jù)每個基因的p值(作為差異表達(dá)測試的一部分)的函數(shù)進(jìn)行排序慎恒。 在這種情況下,分配給基因的p值可以解釋為各組之間基因表達(dá)差異的概率至少與在沒有組間差異的情況下觀察到的差異一樣大撵渡。 我們只需聲明此p值函數(shù)的 rank metric(圖1)

圖1

圖1.導(dǎo)出GSEA Local statistic信息:rank metric 融柬,描述了總共m個樣品的基因表達(dá)的成對比較。 確定n個基因中每個基因的RNA水平趋距。 差異表達(dá)測試為每個基因分配一個p值(P)粒氧,并用于導(dǎo)出表示GSEA rank metric的 Local statistic。 基因列表(L)根據(jù)等級排序

排名指標(biāo)的一個示例是表達(dá)式更改中“方向”的符號(即“上調(diào)”為正节腐,而“下調(diào)”為負(fù))與p值(P)的乘積外盯。



在此rank metric下摘盆,具有相對較小p值的上調(diào)基因出現(xiàn)在列表的頂部,而具有較小p值的下調(diào)基因出現(xiàn)在列表的底部

第一步先根據(jù)表達(dá)矩陣(根據(jù)表達(dá)矩陣不同處理組)計(jì)算每個基因的 Fold change 和 P-value饱苟,然后根據(jù) P-value(Rank Metric)進(jìn)行排名孩擂,Rank Metric 小的在上面,大的在下面

2.Global statistic

Global statistic是GSEA的核心箱熬,用于計(jì)算它的相當(dāng)簡單的算法掩蓋了一種判斷每個候選基因集的相當(dāng)深刻的統(tǒng)計(jì)方法肋殴。 我們將暫時擱置技術(shù)細(xì)節(jié),以了解GSEA如何計(jì)算給定途徑的富集得分坦弟。
讓我們按照圖1所示的過程進(jìn)行操作护锤,這里有一個排名基因的列表。 出于說明目的酿傍,假設(shè)我們希望測試列表中細(xì)胞周期途徑的富集情況(圖2)烙懦。

圖2

圖2.全局統(tǒng)計(jì)的樣本計(jì)算:GSEA富集得分。 該過程需要根據(jù)排名指標(biāo)排序的排序基因列表(L)以及候選基因集(G)赤炒。 在這種情況下氯析,候選者是哺乳動物細(xì)胞周期途徑。 通過從排名列表的頂部開始并依次考慮每個基因來計(jì)算運(yùn)行總和:如果該基因存在于基因集中(紅色莺褒; +)掩缓,則將其加和,否則將其總和(-)減一遵岩。 GSEA富集分?jǐn)?shù)(S)是列表中任何一點(diǎn)的總和的最大值你辣。 盡管未示出,但是行駛總和可能沿負(fù)方向偏離尘执,因此舍哄,S實(shí)際上是行駛總和的最大絕對值
GSEA一次考慮一個候選基因集版仔。為了計(jì)算富集得分些举,GSEA從排名靠前的基因列表的頂部開始。如果一個基因是候選基因集的成員沼侣,那么它會加上一個連續(xù)的總和丧靡,否則就減去它蟆沫。對排名列表中的每個基因重復(fù)此過程,該基因集的富集得分等于運(yùn)行總和所達(dá)到的最大絕對值温治。
我們避開的該算法的一個方面是增加或減少的值饭庞。在GSEA的原始版本中(Mootha 2003),這些值是經(jīng)過特別選擇的罐盔,以使所有基因的總和為零但绕。在圖2中,這相當(dāng)于在列表末尾滿足水平軸的運(yùn)行總和。在這種情況下捏顺,濃縮分?jǐn)?shù)是Kolmogorov-Smirnov(K-S)統(tǒng)計(jì)信息六孵,用于確定濃縮分?jǐn)?shù)是否在統(tǒng)計(jì)上有意義。 Subramanian等人(Subramanian 2005)描述的更新程序使用了該程序的“加權(quán)”版本幅骄,由此運(yùn)行總和的增量與該基因的rank metric成正比劫窒。這些選擇的原因是相當(dāng)技術(shù)性的,我們在下一節(jié)中將其保留給那些更傾向于數(shù)學(xué)的人拆座。

富集分?jǐn)?shù):


注明: n指的是所有基因數(shù)量
這一小節(jié)單獨(dú)討論了對于單個基因集Gk(其中k為索引)主巍,該單基因集Gk表示的是你感興趣的功能基因集,而非你自己的data,那么一共有k個單基因集挪凑;對于每一個基因集Gk來說孕索,里面一共有nk個基因(每個基因用 gkj 表示)
而不存在于基因集Gk的,我們將其表示為:

并且有:


式子中 st 表示 rank metric, nk表示的是Gk基因集的基因數(shù), gene(t)表示你自己的數(shù)據(jù): Ranked gene list(L)
其中1是指定基因集中成員的指標(biāo)函數(shù)躏碳。 我們將在此處注意到搞旭,富集得分是基因集大小的函數(shù),當(dāng)我們在下面討論顯著性檢驗(yàn)時菇绵,它將發(fā)揮作用肄渗。 為了更好地理解這些方程式的含義,讓我們看看改變指數(shù)α?xí)r會發(fā)生什么咬最。
意味著 Sk(GSEA) 表示你自己的data(L 中的)中g(shù)ene(t) 存在于某功能的單基因集Gk中,那么就得分,否則減分

Equal weights: The ‘classic’ Kolmogorov-Smirnov statistic:

考慮當(dāng)α= 0時的簡單情況翎嫡。 然后,來自基因集中基因的所有的貢獻(xiàn)都不會考慮(2)中的 rank metric 永乌。實(shí)際上惑申,這賦予所有基因相等的權(quán)重。


注明: nk指的是Gk集中的基因數(shù)量
如果上升到連續(xù)型變量:

其中X(i)表示累計(jì)值, X(1) = 1 , X(2) = 2 , X(i) = i , X(n) = n ; 以此類推铆遭,并且Xt的值要小于等于你給定的x值
也就是說當(dāng)x = 5時 , 基因list(L) 只能取前五個基因
當(dāng) L 中的基因存在于Gk中時,采用(4)式 ; 當(dāng) L 中的基因不存在于Gk中時,采用(5)式 , 因此我們可以根據(jù)我們的 L 基因集和感興趣通路的基因集Gk , 做一個累計(jì)分布曲線 :

圖3.經(jīng)驗(yàn)累積分布函數(shù)硝桩。 上圖 : 假設(shè)的有序基因列表,其中每個豎線表示基因集中的一個基因(紅色)或不代表(藍(lán)色); 下圖 : 基因組中的基因(紅色)和基因組外的基因(藍(lán)色)的經(jīng)驗(yàn)累積分布函數(shù)枚荣。 ecdfs之間的最大偏差是綠色的虛線
經(jīng)過仔細(xì)檢查,將累計(jì)總和(圖2底部)與ecdfs(圖3)聯(lián)系起來很簡單:累計(jì)總和的增加對應(yīng)于基因集中基因的 ecdf 的增加(圖3啼肩,紅色)橄妆,同樣地 ,總和的減少對應(yīng)于外部基因(藍(lán)色)的 ecdf 增加(實(shí)際上是減分)祈坠。 然后害碾,富集得分(運(yùn)行總和的最大偏差)對應(yīng)于ecdf圖之間最大的垂直距離(綠色虛線)。

GSEA軟件計(jì)算最終的富集得分即為上圖對應(yīng)于ecdf圖之間最大的垂直距離(綠色虛線)赦拘,也稱為峰值

如上圖所示慌随,如果我們將某一個通路(功能)相關(guān)的基因集合定義為 P,自己的數(shù)據(jù)集基因(按照Foldchange或者相關(guān)性排過序)集合定義為 G,GSEA的本質(zhì)是如果某基因在 G與P基因集 中都有阁猜,那么就加分丸逸;否則減分

最終的富集得分(NES)為正,說明相關(guān)基因在 G 的頂部富集剃袍,如果是按Foldchange從高到低排序黄刚,那么說明該通路 (集合P表示為與該通路相關(guān)的基因集合)對應(yīng)存在于基因集 G 中的基因在 G 里面顯示為上調(diào)的居多(紅框所示基因);最終的富集得分(NES)為負(fù)民效,說明相關(guān)基因在 G 的底部富集憔维,如果是按Foldchange從高到低排序,那么說明該通路 (集合P表示為與該通路相關(guān)的基因集合)對應(yīng)存在于基因集 G 中的基因在 G 里面顯示為下調(diào)的居多(紫框所示基因)

這里對應(yīng)存在的意思是某基因在 G與 P基因集 中都有

Kolmogorov-Smirnov擬合優(yōu)度檢驗(yàn)

我們已經(jīng)了解到

可以表示為圖3中紅色ecdfs與藍(lán)色ecdfs之間的最大距離,那么給定的一個顯著的



是否可以理解為兩個不同的 ecdfs 代表不同的樣本 , 這里面是否存在抽樣誤差

單樣本Kolmogorov-Smirnov擬合優(yōu)度檢驗(yàn)
該檢驗(yàn)的目的是檢測實(shí)際樣本的分布是否與理論分布是否相契合

基本原理 : 將理論分布下的累計(jì)頻數(shù)分布與觀測到的累計(jì)頻數(shù)分布相比較,找出它們之間的最大差異點(diǎn)(即圖中綠色的線所在的點(diǎn)) , 并參照抽樣分布,判斷最大差異是否來自偶然

其中:


代表觀測到的累計(jì)頻數(shù)分布 , 而:

代表理論分布下的累計(jì)頻數(shù)分布 , 而它們的原假設(shè)(H0)即為判斷它們是否相等
那么我們構(gòu)造:

Dn代表它們之間的誤差 , 構(gòu)造統(tǒng)計(jì)量 , 判斷Dn是否小于某個臨界值
參考資料傳送門: 傳送門

雙樣本Kolmogorov-Smirnov擬合優(yōu)度檢驗(yàn)

對于雙樣本的檢驗(yàn)來說 :



對于兩個樣本

來說,我們要尋找這兩個分布的累計(jì)曲線的差異 , 即圖3中紅色曲線和藍(lán)色曲線的差異(最大差異為綠色虛線)

我們將這兩者的差異表示為Dab , 構(gòu)造統(tǒng)計(jì)量 , 判斷Dab是否小于某個臨界值 , 是否顯著

3. 顯著性檢驗(yàn)

回顧一下畏邢,GSEA使用了一組用于基因列表的 rank metrics 來計(jì)算候選基因集的富集得分业扒。在這一點(diǎn)上,主要問題是哪些分?jǐn)?shù)表明是富集的舒萎?在假設(shè)檢驗(yàn)行話中程储,我們希望確定每個 global statistic 的統(tǒng)計(jì)意義。我們通過得出一個p值來實(shí)現(xiàn)這一點(diǎn)逆甜,該p值表示觀察到給定富集得分或其他更極端的概率虱肄。為此,我們需要對統(tǒng)計(jì)信息的分布方式有所了解交煞。

零分布
從我們對 global statistic 的討論中咏窿,使用加權(quán)的富集得分使我們沒有對其空分布的解析描述。 也就是說素征,用本地統(tǒng)計(jì)量加權(quán)濃縮分?jǐn)?shù)SGSEAk不同于經(jīng)典的Kolmogorov-Smirnov統(tǒng)計(jì)量集嵌,該統(tǒng)計(jì)量通常遵循類似K-S的分布。
GSEA采用“重采樣”或“自舉”方法來為每個基因集的富集得分得出空分布的經(jīng)驗(yàn)樣本御毅。 GSEA軟件提供了兩種排列方法的選擇根欧,這些方法是空分布計(jì)算的基礎(chǔ)(圖4)。


然后用你自己的數(shù)據(jù)基因集 L 與每一個(一共Nk個)功能基因集計(jì)算富集分?jǐn)?shù) Sk(Nk)GSEA {一共Nk個富集分?jǐn)?shù)}
其中Nk表示被置換的次數(shù) , Sk(Nk)GSEA表示富集分?jǐn)?shù)

圖4. GSEA使用置換方法為每個基因集生成零分布端蛆。 為了簡潔起見凤粗,我們描述了單個基因集的排列方法示意圖。 在GSEA中今豆,對每個基因集分別重復(fù)此過程嫌拣。 A.表型排列。 B.基因組排列呆躲。 C. p值的計(jì)算
圖4 C 中紅色部分表示未經(jīng)過置換的富集分?jǐn)?shù)值
經(jīng)過Nk次置換后,我們一共可以得到 Nk 個富集分?jǐn)?shù), 那么不同的富集分?jǐn)?shù)擬合除一共頻率分布直方圖,作為判斷顯著性的標(biāo)準(zhǔn), 并將頻率分布的期望記為:

表型(功能)置換

對于Phenotype permutation來說,這里的Phenotype gene list(表型基因集)表示的其實(shí)是功能基因集


GSEA描述了利用表型置換的方法對零分布進(jìn)行采樣 , 對于給定的某功能的基因集Gk , 利用表型置換的方法將它們的功能基因集Gk的 "某功能" 標(biāo)簽進(jìn)行置換,從而得到置換分?jǐn)?shù) (圖4 B):

從直觀的角度來看异逐,這在基因等級排名與表型之間沒有特定關(guān)聯(lián)的假設(shè)下生成了一個富集得分分布的樣本。 換句話說插掂,我們了解了富集分?jǐn)?shù)變化的范圍以及兩組有效地相同的頻率灰瞻。
從統(tǒng)計(jì)的角度來看腥例,作者聲稱這為空模型提供了更準(zhǔn)確的描述。
重要的是酝润,類別標(biāo)記的排列保留了基因與基因之間的相關(guān)性燎竖,因此,與通過排列基因獲得的意義相比袍祖,提供了更生物學(xué)上合理的重要性評估底瓣。
確實(shí),GSEA的各種變種蕉陋,其聲稱只是簡單的方法論捐凭,都依賴于基因-基因獨(dú)立性的假設(shè)(Irizarry 2009)。 然而凳鬓,經(jīng)驗(yàn)比較表明茁肠,忽略這些相關(guān)性會導(dǎo)致方差膨脹-更寬和更平坦的零分布-導(dǎo)致基因集錯誤發(fā)現(xiàn)的風(fēng)險更高(Tamayo 2016)。

P-value 計(jì)算

根據(jù) P-value 的定義我們不難發(fā)現(xiàn), P-value的計(jì)算可以是:



根據(jù)定義缩举,P是在原假設(shè)下觀察統(tǒng)計(jì)或更極端的概率

4. 多次校正

當(dāng)我們檢驗(yàn)一個假設(shè)時垦梆,觀察具有較小p值的統(tǒng)計(jì)信息的機(jī)會會增加。如果小于顯著性水平仅孩,則可以將它們錯誤地分類為發(fā)現(xiàn)或 I 類錯誤托猩。多種校正程序試圖對這些進(jìn)行量化和控制。
在GSEA中辽慕,觀察到的數(shù)據(jù)分布是否等于理論數(shù)據(jù)分布這是一個假設(shè)京腥。量化 I 型錯誤的指標(biāo)是錯誤發(fā)現(xiàn)率(FDR)。 FDR被定義為實(shí)際上是被拒絕的零假設(shè)的一部分的期望值溅蛉。實(shí)際上公浪,GSEA憑經(jīng)驗(yàn)確定了這一比例。
通常船侧,給定 Local statistic 的指定閾值 T欠气,F(xiàn)DR是大于T的真實(shí)無效假設(shè)的數(shù)量除以大于T的真實(shí)和假無效假設(shè)的總和。對于GSEA镜撩,T是富集的某個值分?jǐn)?shù)和真實(shí)的無效假設(shè)將代表實(shí)際上未富集的基因集预柒。實(shí)際上,我們不會直接觀察后者袁梗,因此可以根據(jù)超過T的經(jīng)驗(yàn)空值進(jìn)行估算(圖6)卫旱。


圖6.經(jīng)驗(yàn)錯誤發(fā)現(xiàn)率A.觀察到的統(tǒng)計(jì)信息的分布和通過置換方法得出的零分布。 T代表閾值围段。 B. 假設(shè)觀察到n(obs)和經(jīng)驗(yàn)空分布為n(null)的數(shù)據(jù)點(diǎn)。A中右側(cè)分布尾部的放大視圖投放。超出閾值的空分布值的數(shù)量用于估計(jì)真實(shí)的空假設(shè)奈泪。 錯誤拒絕的比例估計(jì)為s_null(紅色)與s_obs(藍(lán)色)之比,并在數(shù)據(jù)點(diǎn)數(shù)不相等時引入校正

如果超過閾值T的觀察到的統(tǒng)計(jì)信息和無效統(tǒng)計(jì)信息的數(shù)目分別為s(obs)和s(null),則經(jīng)驗(yàn)錯誤發(fā)現(xiàn)率為:


標(biāo)準(zhǔn)化的富集得分:考慮基因組大小

從等式(1)-(3)的定義可以看出富集得分取決于基因組大小涝桅。因此拜姿,除非我們處于不可能的情況下,即我們測試的所有基因集大小均相等冯遂,否則富集分?jǐn)?shù)將不會均勻分布蕊肥,因此無法直接進(jìn)行比較。這排除了經(jīng)驗(yàn)錯誤發(fā)現(xiàn)率的計(jì)算蛤肌。
GSEA通過對所計(jì)算的富集得分進(jìn)行轉(zhuǎn)換以使其處于可比較的規(guī)模來解決此問題壁却。特別地,該歸一化的富集得分(NES)是富集得分除以相應(yīng)的空分布的期望值(即裸准,平均值)展东。具體而言,對于每個基因集炒俱,我們通過基因集置換得出富集得分:



和相應(yīng)的空分布:



并且定義, 標(biāo)準(zhǔn)化富集得分:

未經(jīng)過置換后GSEA富集分?jǐn)?shù)的標(biāo)準(zhǔn)化值:
實(shí)際情況中肯存在正負(fù)值, 那么對正負(fù)值分別進(jìn)行標(biāo)準(zhǔn)化:


經(jīng)過置換后GSEA富集分?jǐn)?shù)的標(biāo)準(zhǔn)化值:
實(shí)際情況中肯存在正負(fù)值, 那么對正負(fù)值分別進(jìn)行標(biāo)準(zhǔn)化:

翻譯自:
Gene Set Enrichment Analysis

最后編輯于
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
  • 序言:七十年代末盐肃,一起剝皮案震驚了整個濱河市,隨后出現(xiàn)的幾起案子权悟,更是在濱河造成了極大的恐慌砸王,老刑警劉巖,帶你破解...
    沈念sama閱讀 219,366評論 6 508
  • 序言:濱河連續(xù)發(fā)生了三起死亡事件峦阁,死亡現(xiàn)場離奇詭異谦铃,居然都是意外死亡,警方通過查閱死者的電腦和手機(jī)拇派,發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 93,521評論 3 395
  • 文/潘曉璐 我一進(jìn)店門荷辕,熙熙樓的掌柜王于貴愁眉苦臉地迎上來,“玉大人件豌,你說我怎么就攤上這事疮方。” “怎么了茧彤?”我有些...
    開封第一講書人閱讀 165,689評論 0 356
  • 文/不壞的土叔 我叫張陵骡显,是天一觀的道長。 經(jīng)常有香客問我曾掂,道長惫谤,這世上最難降的妖魔是什么? 我笑而不...
    開封第一講書人閱讀 58,925評論 1 295
  • 正文 為了忘掉前任珠洗,我火速辦了婚禮溜歪,結(jié)果婚禮上,老公的妹妹穿的比我還像新娘许蓖。我一直安慰自己蝴猪,他們只是感情好调衰,可當(dāng)我...
    茶點(diǎn)故事閱讀 67,942評論 6 392
  • 文/花漫 我一把揭開白布。 她就那樣靜靜地躺著自阱,像睡著了一般嚎莉。 火紅的嫁衣襯著肌膚如雪。 梳的紋絲不亂的頭發(fā)上沛豌,一...
    開封第一講書人閱讀 51,727評論 1 305
  • 那天趋箩,我揣著相機(jī)與錄音,去河邊找鬼加派。 笑死叫确,一個胖子當(dāng)著我的面吹牛,可吹牛的內(nèi)容都是我干的哼丈。 我是一名探鬼主播启妹,決...
    沈念sama閱讀 40,447評論 3 420
  • 文/蒼蘭香墨 我猛地睜開眼,長吁一口氣:“原來是場噩夢啊……” “哼醉旦!你這毒婦竟也來了饶米?” 一聲冷哼從身側(cè)響起,我...
    開封第一講書人閱讀 39,349評論 0 276
  • 序言:老撾萬榮一對情侶失蹤车胡,失蹤者是張志新(化名)和其女友劉穎檬输,沒想到半個月后,有當(dāng)?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體匈棘,經(jīng)...
    沈念sama閱讀 45,820評論 1 317
  • 正文 獨(dú)居荒郊野嶺守林人離奇死亡丧慈,尸身上長有42處帶血的膿包…… 初始之章·張勛 以下內(nèi)容為張勛視角 年9月15日...
    茶點(diǎn)故事閱讀 37,990評論 3 337
  • 正文 我和宋清朗相戀三年,在試婚紗的時候發(fā)現(xiàn)自己被綠了主卫。 大學(xué)時的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片逃默。...
    茶點(diǎn)故事閱讀 40,127評論 1 351
  • 序言:一個原本活蹦亂跳的男人離奇死亡,死狀恐怖簇搅,靈堂內(nèi)的尸體忽然破棺而出完域,到底是詐尸還是另有隱情,我是刑警寧澤瘩将,帶...
    沈念sama閱讀 35,812評論 5 346
  • 正文 年R本政府宣布吟税,位于F島的核電站,受9級特大地震影響姿现,放射性物質(zhì)發(fā)生泄漏肠仪。R本人自食惡果不足惜,卻給世界環(huán)境...
    茶點(diǎn)故事閱讀 41,471評論 3 331
  • 文/蒙蒙 一备典、第九天 我趴在偏房一處隱蔽的房頂上張望异旧。 院中可真熱鬧,春花似錦提佣、人聲如沸泽艘。這莊子的主人今日做“春日...
    開封第一講書人閱讀 32,017評論 0 22
  • 文/蒼蘭香墨 我抬頭看了看天上的太陽匹涮。三九已至,卻和暖如春槐壳,著一層夾襖步出監(jiān)牢的瞬間然低,已是汗流浹背。 一陣腳步聲響...
    開封第一講書人閱讀 33,142評論 1 272
  • 我被黑心中介騙來泰國打工务唐, 沒想到剛下飛機(jī)就差點(diǎn)兒被人妖公主榨干…… 1. 我叫王不留雳攘,地道東北人。 一個月前我還...
    沈念sama閱讀 48,388評論 3 373
  • 正文 我出身青樓枫笛,卻偏偏與公主長得像吨灭,于是被迫代替她去往敵國和親。 傳聞我的和親對象是個殘疾皇子刑巧,可洞房花燭夜當(dāng)晚...
    茶點(diǎn)故事閱讀 45,066評論 2 355