10X单细胞(10X空间转录组)之基因网络推断方法之间的优劣势

大家应该经常用SCENIC或者pySCENIC推断单细胞数据的基因网络吧,但是大家知道网络推断方法之间的算法和优劣势么,这一篇我们来分享一下推断基因网络调控的多方法比较,参考文章在Identifying strengths and weaknesses of methods for computational network inference from single cell RNA-seq data,分享一下如何准确推断细胞调控网络。

图片.png

先说一下最终的结论

  • top(PIDC、MERLIN、SCENIC、PEARSON)
  • mid(Inferelator、 SCODE、LEAP、Scribe)
  • bottom(knnDREMI、SILGGM)

Abstract

单细胞 RNA 测序 (scRNA-seq) 通过测量数千个单个细胞的转录组,为不同细胞状态的转录程序提供了无与伦比的洞察力scRNA-seq 分析中的一个新问题是转录基因调控网络的推断,并且已经开发了许多具有不同学习框架的方法。在这里,在考虑不同类型的黄金标准网络和评估指标的情况下,针对人类、小鼠和酵母中 6 个已发布的单细胞 RNA 测序数据集,介绍了最近 11 种网络推理方法的扩展基准研究根据它们的计算要求以及它们恢复网络结构的能力来评估方法。分析发现,虽然没有一种方法是万能的,而且大多数方法都可以适度恢复基于 AUPR 等全局指标的实验衍生交互,但方法能够捕获与所研究系统相关的regulators的目标。基于整体性能,将这些方法分为三个主要类别,并发现information-theoretic and regression-based methods的组合具有普遍较高的性能。分析还评估了基因调控网络推断插补的效用,并发现少数方法受益于插补,这进一步取决于数据集。最后,在可比较的bulk条件下与推断网络的比较表明,从 scRNA-seq 数据集推断的网络通常更好或与bulk推断的网络相当这表明 scRNA-seq 数据集可以成为基因调控网络推断的一种具有cost-effective的方式。分析应该有利于选择执行网络推理的算法,但也主张改进方法和更好的黄金标准,以准确评估哺乳动物系统的regulatory network推理方法。

Introduction

Inference of genome-scale regulatory networks from global mRNA profiles has been a long-standing problems in systems biology and gene regulation.这些网络对于理解细胞类型规范和疾病的机制很重要。 单细胞组学数据的可用性为 reverse engineer transcriptional regulatory networks开辟了新的机会,使我们能够研究多种细胞类型的调控网络。 单细胞 RNA-seq (scRNA-seq) 测量数万个细胞的全基因组表达谱,这大大降低了生成计算网络推理所需的大样本数据集的成本。 此外,它提供了在同一实验中为已知和新细胞群推断细胞类型特异性调节网络的潜力,提供了一种剖析异质细胞群的有效方法 。

scRNA-seq 数据集的可用性推动了许多从这些数据进行网络推理的方法的发展,这些方法使用不同类型的模型,包括高斯图模型、信息论方法、随机森林、常微分方程和布尔网络。方法在包含伪时间或估算的去噪信号方面也有所不同。其中一些方法专门对 scRNAseq 数据的统计特性进行建模,而另一些方法则是对现有方法的bulk数据的改编。在方法开发的同时,有两个基准研究来比较现有的网络推理算法。一项研究表明,现有方法在模拟数据和真实数据上的表现都很差,而另一项研究表明,在模拟数据上表现良好的方法在真实数据上可能相对较差。虽然这些研究很有用,但仍有许多未解决的问题需要解决,以全面了解来自 scRNA-seq 数据集的现有网络推理算法的优缺点。现有的基准测试工作还研究了来自真实数据(<2000 个基因)的相对较少数量的网络,而实际上基因组规模的调控网络可能有 5-10k 个基因。另一个未知数是不同的基准(例如来自 ChIP 与监管机构扰动的基准)会在多大程度上影响性能,以及方法在多大程度上彼此一致以及在黄金标准网络的specific components上。此外,有人建议对从 scRNA-seq 获得的稀疏计数矩阵进行插补,以改进 scRNA-seq 数据集的下游分析,包括识别基因之间的功能关系。然而,目前尚不清楚这是否更普遍地有益于网络推理方法。最后,目前尚不清楚从同一系统的单个细胞与bulk数据集推断出的网络相互比较的程度 。

图片.png

在这里,在人类、小鼠和酵母样本的已发表 scRNA-seq 实验的六个数据集上比较了 11 种网络推理方法。对不同大小数据集的时间和内存消耗计算要求的算法进行基准测试,确定了几种不太可能扩展到全基因组基因调控网络的算法。分析使用不同的全局指标(例如精确召回曲线下的面积和 F 分数)以及局部指标(例如可以准确预测目标的监管机构的数量)来比较算法性能。基于简单相关性的方法在这些指标上表现得同样好或更好。基于全局指标,网络推理方法与随机方法相比性能适中,但局部指标表明该方法可以预测多个相关监管机构的目标。通常,与其他算法相比表现良好的算法往往彼此一致,尽管这取决于数据集的深度。 scRNA-seq 数据集的插补并没有为大多数网络推理方法提供实质性的好处。最后,我们将这些算法在单细胞数据集上的性能与bulk数据进行了比较,并发现了几种情况,其中 scRNA-seq 数据集的推理与改进的推理相关。总体而言,研究对来自 scRNA-seq 数据集的网络推理的最新算法进行了扩展分析比较,确定了不同类型黄金标准的共同优势和劣势,这些优势和劣势应该有利于网络推理社区以改进网络推理

Results

Computing requirements of regulatory network inference algorithms from single cell RNA-seq data

从 11 种专门设计用于从单细胞测序数据推断监管网络的算法以及bulk推断方法开始。除了这些已发表的算法之外,分析还构建了一个网络,其中基因对之间的边缘通过它们在实验范围内的表达相关性进行加权。首先使用从 10 到 8000 个基因的不同数量的基因对这些算法进行了基准测试。估计了每个网络推理算法的内存和运行时间。在比较的方法中,SCHiRM 和 BTR 没有在合理的时间内完成,因此被排除在进一步分析之外。 SCRIBE 和 HurdleNormal 最多可以运行 2000 个基因,但超出这个范围需要更长的时间来运行。大多数可针对多达 8k 个基因运行的算法消耗了多达 15G 的 RAM。例外情况是 SCHiRM、Inferelator 和 PIDC,它们占用了更高的内存。 SCHiRM 内存消耗呈指数增长,被认为对于更大的基因集是不可行的。对于随后的分析,我们从我们的分析中排除了 SCHiRM、BTR 和 HurdleNormal。

图片.png

图片.png

Assessing the performance of scRNA-seq regulatory network inference algorithms based on global metrics

接下来比较了算法在来自三个物种(酵母、小鼠和人类)的真实黄金标准网络上的性能。这些数据集包括两个酵母应激相关数据集、四个小鼠数据集、一个用于树突细胞、三个从细胞重编程到不同培养基产生的胚胎干细胞,以及一个用于人类胚胎干细胞分化的人类数据集。数据集有一系列细胞,从 Gasch 等人的 163 个细胞到 Zhao 等人的 36,199 个细胞的最大数据集。对于每个物种数据集,我们有三种黄金标准,那些来自敲除或敲除调节器,然后进行全局 mRNA 分析,来自 ChIP-chip 或 ChIP-seq 实验的那些,以及来自这两者交集的那些黄金标准。我们使用精确召回曲线下的标准面积和 F 分数评估算法的性能,考虑每个数据集上的每个算法。虽然精确召回曲线 (AUPR) 下的面积考虑了由网络推理算法估计的所有边缘,但 F 分数是使用顶部边缘计算的。为了选择 F-score 计算的边数,考虑了不同的边数并选择了 5k 边,因为这导致了跨算法最稳定的结果。

图片.png

图片.png

图片.png

图片.png

图片.png

图片.png

首先检查了每种算法对来自不同生物体的单个数据集的性能,并根据其 F-score 对每种方法进行排名。 大多数方法在他们的排名中表现出一些变化,但是除了少数例外,性能总体上是一致的。 例如,SILGGM 在大多数数据集上的性能相对较低,但 Sridharan A2S/FBS Perturb+ChIP 除外,它的排名很高。 SCODE 在酵母数据集上表现良好,但在哺乳动物数据集上表现相对较差。 LEAP 在 Shalek 数据集上表现不佳,但通常属于中等排名方法。 算法的相对性能在两个指标 F-score 和 AUPR 之间基本一致。 在不同的方法中,PIDC、SCENIC、MERLIN 和 Pearson 相关在基于 F-score 和 AUPR 的数据集上的性能最稳定。

根据不同类型的黄金标准检查了算法的性能。 与 ChIP 黄金标准相比,Perturb 黄金标准数据集的 F 分数和 AUPR 通常高于随机数。 例外是 SCODE,与 Knockdown 相比,它更适合 ChIP。 接下来,我们根据每个黄金标准的中位数排名汇总所有数据集的算法,并根据其在两种黄金标准中的整体中位数排名对每种方法进行排名。 根据跨数据集的所有算法的中位数排名,相关性、SCENIC、MERLIN 和 PIDC 在使用 F-score 或 AUPR 时也名列前四。 当使用 F-score 和 AUPR 检查随机网络性能的individual fold improvements时,这种趋势很明显。

Assessing the performance of scRNA-seq regulatory network inference algorithms based on subnetwork metrics

AUPR 和 F-score 提供了监管网络全局准确性的量化,并一次比较一个边缘。然而,从生物学的角度来看,了解是否有一些转录因子 (TF) 的目标与其他转录因子相比更容易预测可能是有益的。此外,跨方法的比较可以确定特定方法是否对某些 TF 更敏感。因此,我们接下来通过测量预测单个转录因子的目标的能力,根据 Siahpirani 等人描述的预测集中 TF 真实目标的倍数富集来关注推断网络的更细粒度的视图。我们使用 FDR 校正的超几何测试为重叠分配了 p 值,并使用显着可预测的 TF 的数量作为网络性能的另一种量化。与算法性能通常接近随机的 AUPR 和 F-score 不同,随机基线的可预测 TF 的数量不超过 1,通常为 0。在数据集上,这些方法排名一致,除了 LEAP在 Shalek 中表现更差,而在 Gasch 中表现更好的 SCODE。当使用中位数对数据集进行聚合时,我们发现 SCENIC 和 Pearson 是排名靠前的方法,其次是 MERLIN 和 PIDC。综合考虑 AUPR、F-score 和可预测 TF 的所有三个指标,表现最好的方法是 SCENIC、MERLIN 和 Correlation。

图片.png

图片.png

接下来检查了每个数据集哪些特定的 TF 是通过不同的方法预测的。此信息可以帮助确定网络的某些部分是否可以通过特定类别的方法或所有方法进行预测。对于在胚胎干细胞状态谱系规范过程中描述转录程序的 Han 数据集,发现几个 ESC 特异性调节因子,如 POU5F1、SOX2 和 NANOG,在使用 ChIP 黄金标准的多种方法中作为可预测的 TF。这些是负责建立此状态的关键 ESC TF。还发现了一些谱系特异性 TF,如 CDX2 和 TBX3,它们是通过多种方法预测的不同谱系的主要调节因子。在比较 Perturb 金标准时,发现了与 ChIP 金标准类似的regulators,但另外还有用于视网膜和造血谱系的 OTX2 和 GATA3 等regulators。对于小鼠细胞重编程数据集,发现许多方法都具有类似的行为,这些方法能够在 ChIP 和 Perturb 黄金标准中识别多种发育调节因子,例如 Pou5f1、Esrrb 和 Sox2。 SCODE 在 Han 数据集中有许多可预测的 TF,但是,其中许多是general regulators,例如 SP1、YY1、POL2RA、ATF2。重要的是,与树突细胞数据集 (Shalek) 相比,发现不同方法的一组不同的regulators一致地over-representation。特别是,这包括与免疫反应相关的调节因子,如 Rel、Nfkb、Stat1 和 Stat3,并通过多种方法进行鉴定,包括 SCENIC、MERLIN、PIDC、Scribe、Inferelator 和 Pearson。发现酵母数据集有类似的行为,这些方法一致地能够恢复与应激反应相关的关键调节因子,如 HAP4(氧化应激)和 GCN4(氨基酸饥饿)。与哺乳动物数据集相比,酵母数据集的可预测 TF 的总体富集倍数更高。在两个数据集之间,与 Gasch 相比,Jackson 表现出更丰富的目标,这可能是因为 Gasch 数据集中的细胞数量较少。总的来说,这表明虽然这些方法的 AUPR 不大,但这些方法能够始终如一地恢复特定系统的相关regulators

图片.png

图片.png

图片.png

图片.png

图片.png

图片.png

最后,为了更深入地了解每种方法恢复哺乳动物调节网络的能力,从文献中报道的调节相互作用中策划了small黄金标准网络。调节网络有 35 个调节子和 90 个目标基因(其中一些也包括调节子),由 267 个边连接。接下来,使用了从 Tran A2S 数据集推断出的网络中的前 5000 条边,并在 267 次交互中的 199 次中描绘了真阳性(通过一种方法发现)和假阴性(通过一种方法遗漏),由于缺少表达,在去除边缘后仍然存在.大多数方法恢复了 Nanog、Sox2 和 Pou5f1 之间的交叉调节相互作用,它们是建立 ESC 状态的主要调节器。恢复最少真阳性数的方法包括 SCODE、knnDREMI 和 PIDC。有趣的是,像 PIDC 一样基于信息理论的 SCRIBE 与 PIDC 相比能够推断出更多的真阳性。总体而言,基于不是真正正边缘恢复的性能与分析的全局和可预测的 TF 指标一致,但也突出了网络推理方法在策划监管交互恢复方面的其他属性

Defining common and method-specific network components

接下来query算法在他们的预测中的一致程度。为此,测量了每对推断网络的前 500、5,000 和 50,000 条边之间的 Jaccard 相似度和 F-score。发现无论包含多少顶部边缘和相似性度量,相似性模式通常都是相似的。随着考虑更多边缘,相似性的大小通常会增加。专注于在前 5k 边缘获得的 Jaccard 分数,以与我们的其他结果保持一致。我们使用两种顺序对方法进行分组。首先,我们根据为每个数据集推断的网络的中值相似度对方法进行聚类。这使我们能够跨数据集将所有方法相互比较。与其他数据集相比,一些数据集产生了更多相似的网络。例如,在 Han hESC 数据集上,与具有较低相似性的 Shalek 和 Gasch 相比,最大 Jaccard 系数最高。这两个数据集也是两个数据集中最小的,这可以解释较低的相似性。与其他方法相比,SILGGM 和 knnDREMI 通常学习不同的网络,包括基于同一类模型的网络,例如knnDREMI (PIDC, Scribe) 的信息论和 SILGGM (MERLIN, Inferelator) 的概率模型。 Scribe 与其他方法的相似性取决于数据集,然而,在大多数情况下它与 PIDC 最相似,这与依赖信息论度量的两种方法一致。 LEAP 与 Correlation 最相似,但通常识别不同的网络。其中相似度最高的方法是 Pearson、SCENIC、Inferelator 和 PIDC。

图片.png

图片.png

图片.png

图片.png

图片.png

图片.png

接下来根据每个数据集上的网络 Jaccard 相似度对方法进行聚类。 后者揭示了更多特定于数据集的分组。 特别是,发现 LEAP 和 Pearson 是 Gasch 和 Zhao 数据集中最相似的方法,而 PIDC 在七个数据集中的四个(Jackson、Tran (A2S)、Tran (FBS) 和 Han)中与 Pearson 最相似。 还发现 LEAP 和 SILGGM 在 Shalek 上形成了一个集群,而另一组包括 SCRIBE、PIDC、SCENIC、Pearson、MERLIN 和 Inferelator。 总之,这些比较表明,性能最好的算法倾向于学习相似的网络,然而,PIDC 和 Pearson 等方法经常学习相似的网络。 此外,方法之间的一致性和差异可以由数据集的性质决定,包括样本大小和实验条件

Examining the impact of imputation on scRNA-seq network inference

单细胞表达数据的特点是有大量的零表达计数,这可能是由于技术(低深度)或生物学原因而产生的,这是基因在细胞中真正不表达。为了解决这个问题,已经开发了几种方法来估算缺失表达计数的值。我们将 MAGIC 算法应用于每个数据集,这是最近基准研究中最重要的插补方法之一。 MAGIC 计算单元格之间的成对相似性,并基于这些相似性创建马尔可夫转移图。然后它基于这个马尔可夫转换图在相似的单元格之间“扩散”表达计数。我们使用来自每个实验的估算数据来推断网络,并将它们的 AUPR、F 分数和可预测的 TF 指标与从稀疏数据中推断出的网络进行比较。基于 F 分数,大多数方法并没有从跨数据集的插补中受益。一个例外是 Shalek 和 Gasch 数据集,正如我们之前提到的,它们属于较小的数据集。我们汇总了不同黄金标准和数据集的 F 分数值,并比较了与插补前的比率,发现使用插补数据时 SCODE 和 kNN-DREMI 似乎有所改善。然而,我们发现插补通常不利于网络推理过程。我们根据 AUPR 和可预测的 TF 重复了这些比较。 AUPR 在有和没有插补的情况下通常没有太大变化,但是,我们注意到 Shalek 数据集的性能有类似的提升。最后,基于可预测的 TF,我们发现大多数算法并没有从跨数据集的平滑中受益,Tran 数据集为 SCODE 和 knnDREMI 方法提供了适度的好处。综上所述,我们在不同实验数据集和黄金标准上的实验表明,插补可以在少数算法和数据集的情况下受益,但通常对大多数方法的网络推理没有帮助。


图片.png

图片.png

Discussion

单细胞 RNA-seq 数据集的快速增长为基于表达的基因调控网络推断领域开辟了巨大的机遇。 因此,已经投入了大量努力来开发和应用监管网络推理算法到单细胞数据集。 在这里,我们对跨越不同物种和模型细胞系的大量数据集的计算要求和网络结构恢复方法的整体性能进行了基准测试。 我们的努力通过考虑更多的数据集、额外的方法类别和多个本地和全球网络评估指标来扩展先前的努力,这些指标可以为不同方法所做的预测提供更多的生物学洞察力,并揭示方法学和资源生成进步的潜在方向,以改进网络恢复。

根据计算要求以及在不同黄金标准网络上的性能评估方法。由于过多的计算要求,无法对几种专门用于处理单细胞数据集统计特性的方法进行排名,例如 SCHiRM、HurdleNormal 和 BTR。随着单细胞数据集的增长,算法的有效实施将变得很重要。我们的评估指标既包括 AUPR 和 F-score 等全球指标,也包括可显着预测目标的数量regulators等本地指标。基于 AUPR 和 F 分数,方法的整体性能仍然略好于随机,但是,发现可预测的 TF 作为更敏感的指标,突出了网络推理方法的优势。重要的是,不同的方法能够概括感兴趣系统的相关regulators,例如发育数据集中的关键干细胞regulators和树突细胞数据集中的免疫反应regulators。没有一种方法能在所有数据集和黄金标准中获胜,但是,根据方法的整体性能和计算要求,将它们分为三大类:顶级(PIDC、MERLIN、SCENIC、PEARSON)、中(Inferelator、 SCODE、LEAP、Scribe)和底部(knnDREMI、SILGGM)的一组方法

scRNA-seq 数据集的一个挑战是零的高比例,这可能是由于生物学和技术原因造成的。 由于已经提出插补作为处理高度稀疏数据集(例如 scRNA-seq 数据集)的方法,因此研究了插补对推断网络质量的影响。 由于这种插补并没有提高大多数方法的性能,但是,它确实受益的数据集往往是那些细胞数量相对较少的数据集。 我们分析中的一个警告是,我们只考虑了一种插补方法 MAGIC,它被证明是最重要的插补方法之一。 未来工作的一个方向是考虑额外的插补方法并进行额外的实验来检查插补对网络推理的影响。

单细胞转录组数据集的优势在于,单个实验可以产生大量样本,这些样本与已用于网络推理的现有bulk数据集相当或更大。 因此,我们使用bulk和单细胞数据集的方法比较了酵母、小鼠和人类的bulk和单细胞 RNA-seq 数据集的推断网络的质量。 发现 scRNA-seq 数据集尽管很稀疏,但在使用bulk RNA-seq 数据集时能够捕获足够的有意义的生物变异并表现出同等水平。 然而,我们的研究并不完美,因为bulk和单细胞数据集是从不同来源收集的。 生成的受控数据集捕获同一系统的bulk和单细胞profiles,可以进一步了解 scRNA-seq 数据集在推断基因调控网络方面的相对优势

目前的研究比较了仅依赖表达的方法。来自bulk数据的实验表明,结合先验信息来约束网络结构以及估计隐藏的监管活动可以有益于网络推理。未来工作的一个方向是开发将先验知识纳入推断网络估计的基准方法。在工作中收集的黄金标准和数据集应该有利于这些未来的研究。未来工作的另一个方向是利用 scRNA-seq 数据集的固有异质性和种群结构。基于多任务学习的方法是一个很有前景的框架来模拟人口和网络的异质性。发现的一个令人惊讶的发现是一个简单的基于 Pearson 相关性的指标的相对较好的表现。这可能是我们公认不完美的黄金标准的产物。基于更新的高通量扰动研究(如 Perturb-seq 和 Perturb-ATAC)产生改进的黄金标准,特别是对于哺乳动物系统,可以显着提高我们推断基因组规模基因调控网络的能力

Method

Regulatory network inference algorithms for single cell RNA-seq datasets

图片.png

图片.png

图片.png

生活很好,有你更好

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
禁止转载,如需转载请通过简信或评论联系作者。
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 212,185评论 6 493
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 90,445评论 3 385
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 157,684评论 0 348
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 56,564评论 1 284
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 65,681评论 6 386
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 49,874评论 1 290
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 39,025评论 3 408
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 37,761评论 0 268
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,217评论 1 303
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,545评论 2 327
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,694评论 1 341
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 34,351评论 4 332
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 39,988评论 3 315
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,778评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,007评论 1 266
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 46,427评论 2 360
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 43,580评论 2 349

推荐阅读更多精彩内容