10X空间转录组分析之转录组信息 & 空间位置的联合分析(SEDR)

hello,大家好,今天给大家继续分享10X空间转录组的分析,其实之前分享的内容已经多次强调过空间转录组一定要转录组信息和空间位置联合进行分析,大家可以参考之前我的文章10X空间转录组空间高变基因联合组织区域识别之SpatialDE210X空间转录组聚类分析之图卷积网络(graph convolutional network)10X空间转录组聚类分析之BayesSpace算法聚类等,就不给大家一一列举了,大家做研究的一定好多多学习。

今天给大家分享的联合方法的文章在Unsupervised Spatially Embedded Deep Representation of Spatial Transcriptomics,下面放一张分析的原理图

图片.png

原理图大家也可以看到,需要计算临近矩阵,其实就是联合空间位置进行下游分析,非常值得大家的借鉴,好了,我们开始分享今天我们的分析内容,最后看一看示例代码。

Abstract

空间转录组学使我们能够剖析组织异质性并绘制出细胞间通讯(这里的通讯主要就是临近通讯)。转录组学数据和相关空间信息的最佳整合对于充分利用数据至关重要(转录组 & 空间位置)。作者开发了 SEDR,这是一种转录本和空间信息的无监督空间嵌入深度表示(转录组 & 空间位置联合分析)。 SEDR pipeline使用深度自动编码器构建基因表达的低维潜在表示,然后通过变分图自动编码器同时嵌入相应的空间信息。将 SEDR 应用于人类背外侧前额叶皮层数据并获得了更好的聚类精度,并通过轨迹分析正确追溯了产前皮层发育顺序。还发现 SEDR 表示非常适合批量集成。将 SEDR 应用于人类乳腺癌数据,在视觉上同质的肿瘤区域内识别出异质子区域,识别具有促炎微环境的肿瘤核心和富含肿瘤相关巨噬细胞的外环区域,该区域驱动免疫抑制微环境。

Introduction

单细胞组学技术能够以单细胞分辨率进行测量,并导致在健康和患病状态下的各种组织中发现新的亚群。然而,在高通量组学数据采集之前将组织分解为单个细胞会导致细胞空间信息丢失,从而阻碍我们剖析单个细胞的空间组织和细胞间相互作用的能力。虽然已经开发出计算工具来从配体和受体表达预测细胞间相互作用,但它们需要使用免疫组织化学 (IHC) 或免疫荧光 (IF) 进行验证。新兴的空间组学技术通过同时测量基因/蛋白质表达和细胞的空间位置来克服这些限制。这种组织学组织的空间解析转录组能够重建组织结构和细胞间相互作用。 这种方法已被证明在许多应用中很有价值,包括对脑部疾病、肿瘤微环境和胚胎发育的研究。
在当前可用的空间转录组学方法中,基于原位捕获的技术,例如 10x Genomics Visium 和 Nanostring GeoMX DSP,由于它们的accessibility和在每个spot内分析大量 mRNA 目标的能力而受到欢迎。原则上,来自组织样本的组织切片被透化,释放的 mRNA 被载玻片表面上空间排列的寡核苷酸或手动定义的感兴趣区域 (ROI) 中的预杂交 RNA 目标条形码捕获。然而,这两种技术都受到 mRNA 捕获区域的限制,最小直径通常为~50μm,比单个细胞大。为了克服这个问题,已经开发了几种计算方法来对空间点的细胞混合物进行去卷积。最近,mRNA 捕获方法的改进导致直径约 1-10μm 的亚细胞捕获区域更小。这些高分辨率空间转录组学方法可以获得具有更高空间保真度的空间解析转录组,而不会影响捕获的基因数量。它们包括 Slide-seq、DBiT-seq、Stereo-seq、PIXEL-seq 和 Seq-Scope,后三者迄今为止获得的最高分辨率(~1μm)。这些亚微米分辨率的方法通常需要voxel binning或细胞分割,以产生用于下游分析的基因-细胞表达矩阵。捕获区域的大小也得到了改善,从而提高了整体细胞通量,需要可以处理大空间数据的新计算方法。
在分析空间转录组学数据时,结合基因表达和空间信息来学习每个细胞或spot的判别式表示至关重要。然而,已建立的工作流程,例如 Seurat,仍然采用专为单细胞 RNA-seq 分析设计的管道,主要关注基因表达数据,而忽略了空间邻域的结构关系。最近,已经开发了几种新的空间转录组学方法来克服这一限制。例如,BayesSpace 从马尔可夫随机场 (MRF) 先验开始,该先验假设属于相同细胞类型的点应该彼此更接近,并使用贝叶斯方法更新模型。 Giotto 实施了一个隐马尔可夫随机场 (HMRF) 模型,通过比较细胞及其邻居之间的基因表达来检测具有连贯模式的域。 SpaGCN 结合空间距离和组织学差异来构建spot的加权图,然后使用图卷积网络 (GCN) 将该图与基因表达相结合以对spot进行聚类。 stLearn 利用对图像的深度学习模型来提取形态特征,并在其上计算形态距离。然后, stLearn使用形态距离和空间邻域信息根据每个spot的识别邻居来规范化每个点的基因表达。然后将归一化的基因表达用作线性主成分分析 (PCA) 的输入,然后是均匀流形近似和投影 (UMAP) 以及空间聚类。值得注意的是,这些方法主要依靠 PCA 来提取基因表达数据的高度可变特征,其中涉及线性变换,因此无法对复杂的非线性关系进行建模。虽然 stLearn确实利用了深度学习,但它仅应用于图像模态,并且该模型仍然依赖于线性 PCA 从空间归一化的基因表达数据中提取特征。此外,这些方法不会产生联合嵌入的基因表达和空间信息的低维表示。基因表达和空间信息的联合嵌入对于有效整合两种模式以实现更好的可视化、聚类、轨迹推断和批量集成至关重要。
在研究中,作者开发了一种无监督的空间嵌入深度表示 (SEDR) 方法,用于学习嵌入空间信息的基因表达的低维潜在表示。 SEDR 模型由两个主要部分组成,一个是用于学习基因表示的深度自编码器网络,另一个是用于嵌入空间信息的变分图自编码器网络。 这两个组件联合优化以生成用于空间转录组学数据分析的潜在表示。 在 10x Genomics Visium 空间转录组学和 Stereo-seq 数据集上应用了 SEDR,并证明了它能够为各种后续分析任务实现更好的表示,包括聚类、可视化、轨迹推断和批量效应校正。

Results

Overview of SEDR.

SEDR 在具有联合嵌入空间信息的低维潜在空间中学习基因表示(下图)
图片.png
给定空间转录组学数据,SEDR 首先使用深度自动编码器网络学习从基因表达空间到低维特征空间的非线性映射。 同时,利用变分图自动编码器将基因表示与相应的空间邻域关系聚合以产生空间嵌入。 然后,将基因表示和空间嵌入连接起来,形成用于重建基因表达的最终潜在表示。 此后,采用无监督的深度聚类方法来增强学习到的潜在表示的紧凑性。 这种迭代深度聚类生成一种软聚类形式,利用特定于聚类和特定于细胞的表征学习之间的推论,为每个细胞分配特定于聚类的概率。 最后,学习到的潜在表示可以应用于各种分析内容。
为了对 SEDR 与其他方法进行定量比较,下载了 10x Genomics Visium 空间转录组学数据和 LIBD 人类背外侧前额叶皮层 (DLPFC) 数据的手动注释层。 LIBD 数据包括来自人类 DLPFC 的 12 个切片,跨越 六个皮质层加上白质。 选择这个数据集是因为人类 DLPFC 具有明确和建立的形态边界,可以作为基本事实。 首先应用标准 Seurat pipeline仅使用表达谱来处理和聚类细胞,并将结果设置为基准基线,以研究空间信息改善细胞聚类的程度。 由于 Giotto、stLearn、SpaGCN 和 BayesSpace 集成了空间信息和 RNA-seq 数据以进行聚类,还将它们与推荐的默认参数一起应用于同一数据集,以针对 SEDR 进行基准测试。
在具有 3,639 个spot和 33,538 个基因的脑切片 151673(下图)中,SEDR 和 BayesSpace 在层边界和调整后的Rand index (ARI,关于Rand index,大家可以参考文章机器学习的评价指标-Rand index) 方面均取得了最佳性能。比较所有 12 个 DLPFC 样本的结果时,SEDR 的平均 ARI (0.427) 位居第二(下图,右下角),但 SEDR 与表现最佳的 BayesSpace (0.457) 之间的差异并不显着(Mann-Whitney U Test : pvalue=0.55)。
图片.png
需要注意的是,BayesSpace 的聚类算法针对空间组学进行了优化,而 SEDR 是一种降维方法,其目标是找到最佳的latent representation。 SEDR 之后是 Leiden 聚类,它不是专门为聚类空间组学而设计或优化的,它实现了与 BayesSpace 相当的聚类性能。这表明 SEDR latent representation有效地整合了基因表达和空间信息以捕获cluster间差异。将 SEDR 与更适合空间组学的聚类算法相结合,有望进一步提高聚类精度。此外,与不产生latent representation的 BayesSpace 相比,SEDR 衍生的嵌入不仅可以用于聚类,还可以用于各种下游分析任务,例如 UMAP 可视化、轨迹推断和批量效应校正,从而提供更大的灵活性和实用性。与 SEDR 类似,SpaGCN 也使用 GCN 来处理空间转录组学数据。此外,它还包含了 SEDR 中未包含的组织学信息。但是,SEDR 的聚类性能优于SpaGCN (Mann-WhitneyU Test,p-value < 0.05)。 stLearn 也集成了组织学数据,但性能同样较差。这可能表明SpaGCNstLearn当前用于合并组织学数据的方法不是最佳的。为了充分利用组织学信息,可能需要将其视为单独的数据模态,并使用专用的多视图算法进行集成。
SEDR 生成一组低维表示特征,可用于各种下游分析,例如轨迹推断。在实验中,我们使用 Monocle3 对样本 151673 执行轨迹推断,其中包含 Seurat 输出(仅限 RNA)和低维 SEDR 表示特征。发现 SEDR 显示出比 Seurat 显着提高的性能(下图)。
图片.png
在 SEDR 输出的 UMAP 图中,属于不同层的细胞被很好的表征,当选择白质 (WM) 作为root时,伪时间反映了皮层层的正确“由内而外”的发育顺序(上图)。这表明,与仅 RNA 分析相比,结合空间信息使 SEDR 能够生成更好的潜在表征,总结空间转录组学数据。我们使用另一种名为基于分区的图形抽象 (PAGA) 的轨迹推理方法进一步证实了观察结果,使用 SEDR 派生的潜在空间嵌入而不是 UMAP 坐标(下图)。
图片.png
PAGA 结果表明,相邻的皮质层往往具有更大的相似性,这表明空间相邻性与转录组学甚至功能相似性有关。值得注意的是,轨迹与皮层发育的时间顺序一致。然后比较了使用 Seurat 衍生的主成分和 SEDR 嵌入生成的 PAGA 图。对于 12 个 DLPFC 切片中的每一个,计算了相邻皮质层之间的边缘权重与所有边缘权重总和的比率。我们发现 SEDR 的比率明显高于 Seurat(Mann-Whitney U 检验 p 值 < 0.05)。

SEDR corrects for batch effects(注意这里是空间数据矫正批次效应)

空间组学应用的激增正在不同实验室产生越来越多的空间解析组学数据。然而,在尝试就空间解析的组织图谱达成共识时,protocols和技术的差异使比较和数据整合变得复杂。与单细胞 RNA-seq (scRNA-Seq) 一样,消除空间组学数据集中的批次效应是一项重大挑战。迄今为止,还没有可用的方法。在这里,作者证明 SEDR 可以学习跨多个批次的联合嵌入并将它们投影到共享的潜在空间中。此外,SEDR 采用深度嵌入聚类 (DEC) 损失函数,使其能够保留生物变异,同时减少技术变异。评估了 SEDR 在 DLPFC 数据集上的批量校正性能。我们首先评估了十二个数据集之间的批次变化,并选择了三组(151507、151672、151673),它们表现出显着的批次效应。如 UMAP 图(下图)所示,不同批次的共同皮质层被分开。
图片.png
由于其在 scRNA-seq 数据集成方面的卓越性能,首先应用 Harmony 来消除批次效应。Harmony能够在保持不同层分开的同时混合批次。然而,当放大各个层时,发现了不同的批次特定子集群,这表明批次效应并未完全消除(下图)。
图片.png
接下来,测试了 SEDR,发现批次效应显着降低(下图)。
图片.png
批次间的公共层非常接近并且对齐良好,而不同层的混合最少。Harmony在 SEDR 嵌入上的进一步应用将批次均匀混合,同时保持层之间的分离(上图D)。值得注意的是,特定于批次的cluster不再存在于各个层中。这表明 SEDR 与Harmony 的结合有效地消除了批次效应。在其他空间组学分析方法中,只有 stLearn能够产生潜在空间嵌入,可以将其输入 Harmony 进行批次校正。因此,针对stLearn对 SEDR 进行了基准测试。由于stLearn 需要将组织学图像作为输入,因此无法将不同批次联合投影到共享的潜在空间,因此从每个数据集生成潜在空间嵌入,然后将它们连接起来进行 Harmony 集成。结果表明批次没有很好地混合并且层分离很差(上图E)。综上所述,SEDR结合Harmony优于单独使用HarmonystLearn with Harmony,可以作为空间组学数据批量校正的有效方法。

Dissecting tumor heterogeneity and immune microenvironments using SEDR.

癌症中的肿瘤内异质性使有效的治疗方案复杂化,并且与较差的生存前景相关。空间转录组学是剖析和表征肿瘤内异质性和肿瘤免疫串扰的有效工具。在人类乳腺癌的 10x Visium 空间转录组学数据上测试了 SEDR,该数据以其高的瘤内和瘤间差异而闻名。为了帮助解释 SEDR 结果,根据 H&E 染色进行了手动病理标记。需要注意的是,与具有明确和既定形态边界的大脑皮层不同,肿瘤组织具有高度异质性并包含复杂的微环境。仅基于肿瘤形态的手动标记不足以表征这种复杂性。根据病理特征,手动将组织学图像分割为 20 个区域,然后将其分为四种主要形态类型:导管原位癌/小叶原位癌 (DCIS/LCIS)、健康组织 (Healthy)、浸润性导管癌 (IDC) ),以及具有低恶性特征的肿瘤周围区域(肿瘤边缘)(下图左上角)。在视觉上,所有五种方法都在宏观层面上与手动注释一致(下图)。尽管如此,与其他方法相比,SEDR 聚类呈现出更平滑的分割,而由 Seurat、stLearn 和 SpaGCN 派生的那些聚类出现了不规则边界的碎片。值得注意的是,SEDR 在肿瘤区域内发现了更多的subcluster,而其他方法倾向于将健康区域划分为subcluster,因为所有方法都设置为生成相同数量的cluster。具体而言,在看似同质的肿瘤区域 DCIS/LCIS_3 内,SEDR 将外“环”(下图,SEDR cluster 7)与肿瘤核心(下图,SEDR cluster 3)分开。这些 SEDR 簇在视觉上同质的肿瘤区域内表明转录和空间上不同的区室。除了聚类分析之外,我们还采用了基于 Seurat3“anchor”的集成工作流程,将注释从人类乳腺癌的 scRNA-seq 参考数据概率转移到空间数据。对于每个点,获得了每个 scRNA-seq 派生类的概率分类(下图B)。转移的类别概率能够描绘出肿瘤区域和存在免疫细胞或成纤维细胞的区域,这有助于进一步剖析肿瘤微环境。
图片.png
为了进一步表征 DCIS/LCIS_3 区域的 SEDR cluster 3(肿瘤核心)和cluster 7(肿瘤边缘)之间的转录差异,进行了差异表达分析,然后进行了通路富集分析(下图)。
图片.png
在第 3 组中,观察到干扰素信号通路(IFIT1、IFITM1、IFITM3 和 TAP1)和 NK 或中性粒细胞活性(FCGR3B 和 TNFSF10)的上调。此外,该区域的 RHOB 上调,表明转移潜力降低。第 3 组代表了促炎免疫反应限制癌症生长的区域。另一方面,在第 7 组中,观察到 TAM、记忆 B 细胞(IGHG1、IGHG3、IGHG4、IGLC2 和 IGLC3)和成纤维细胞(COL1A1、COL1A2、COL3A1、COL5A1、COL6A1)的存在、COL6A2 和 FN1)。上调的组织蛋白酶活性(CTSB、CTSD 和 CTSZ)和补体途径(C1QA、C1S)表明该区域中的 TAM 具有促肿瘤活性。肌动蛋白细胞骨架信号的上调也表明cluster 7 具有更高的转移潜力。此外,上调的组织蛋白酶活性和金属蛋白酶抑制剂(TIMP1 和 TIMP3)也表明细胞外基质完整性受到干扰。总体而言,cluster 7代表了一个具有免疫抑制的促肿瘤微环境和高癌症转移潜力的区域。
许多驱动力被假设为导致肿瘤细胞从浸润前状态转移到浸润癌的原因,包括促肿瘤免疫微环境和肿瘤内细胞间相互作用的减少。在这项研究中,我们使用 PAGA 来推断手动注释的肿瘤区域之间的相互关联性,以追踪转移转移过程。使用 SEDR 嵌入生成的 PAGA 图表明 DCIS_LCIS_3 与相邻的 IDC_6 区域相关(下图)。
图片.png
与所有其他 DCIS_LCIS 区域相比,DCIS_LCIS_3 的差异表达基因 (DEG) 和富集途径表明 DCIS_LCIS_3 具有更多的免疫浸润,尤其是肿瘤相关巨噬细胞 (TAM),而其他 DCIS_LCIS 区域主要由活跃的分裂/循环上皮细胞具有上调的糖酵解和代谢过程。由于 TAM 浸润促进肿瘤血管生成并诱导肿瘤迁移、侵袭和转移,因此已知 TAM 浸润与实体瘤患者的不良生存率密切相关。因此,进行了 Monocle3 分析以推断从 DCIS_LCIS_3 到 IDC_6 过渡的伪时间。由于 DCIS_LCIS_3 和 IDC_6 与 SEDR cluster 3、7 和 11重合,我们在这三个cluster上应用 Monocle3,并将cluster 3(肿瘤核心)设置为起点。 Monocle3 分析表明,与来自 Seurat PCA 嵌入的伪时间相比,来自 SEDR 嵌入的伪时间更好地追踪了由内而外的肿瘤进展。我们随后确定了沿 Monocle3 伪时间改变表达的基因,并揭示了沿轨迹的基因调控的连续波(下图)。
图片.png
总之,SEDR 分析揭示了视觉同质肿瘤区域内的肿瘤内异质性,并揭示了具有 TAM 浸润和癌症相关成纤维细胞 (CAF) 的肿瘤外环(cluster 7),SEDR 还启用了映射 从肿瘤核心到其邻近侵袭区域的分子轨迹,证明了从促炎到免疫抑制微环境的转变,这可能有助于肿瘤转移。

SEDR can handle spatial transcriptomics of high resolution.

当前可用的空间组学技术,包括 10x Visium Spatial Omics、Nanostring GeoMX DSP、SLIDE-seq 和 DBIT-seq,不提供单细胞分辨率,每个捕获点包含 1 到 10 个细胞。然而,Stereo-seq、PIXEL-Seq 和 Seq-Scope 等新兴方法可以实现亚微米分辨率,从而实现亚细胞分辨率。随着技术的不断进步,每个组织检测到的空间分辨率和细胞数量将显着提高,从而产生具有高吞吐量的大型数据集。

Method

图片.png

图片.png

图片.png

图片.png

最后了,我们来看看代码,网址在SEDR

安装
conda create -n SEDR_Env python=3.8 pip
conda activate SEDR_Env
pip install -r requirements.txt
SDER 使用 anndata(基于 Scanpy)作为输入,并输出一个潜在的表示,保存在 SED_result.npz 中。 用户可以将 Python 中的 SEDR 功能提取为:
saved_result = np.load(os.path.join(save_path, "SED_result.npz"), allow_pickle=True)
sed_feat = saved_result["sed_feat"].item()
...
or in R with 'reticulate' package as:
library(reticulate)
np <- import("numpy")
saved_result <- np$load("SED_result.npz")
sed_feat = saved_result$f[["sed_feat"]]
...
We also provide three examples:
  1. run_SEDR_DLPFC_all_data.py: The demo code for 'LIBD human dorsolateral prefrontal cortex (DLPFC, http://research.libd.org/spatialLIBD/)' data. The raw data are publicly available from the Globus endpoint ‘jhpce#HumanPilot10x’ that is also listed at http://research.libd.org/globus. It will take around 40 minutes to run all DLPFC samples on GPU.
    图片.png

    run_SEDR_10x_Genomics_Visium.py: The demo code for 'Visium 10x Genomics' data (https://support.10xgenomics.com/spatial-gene-expression/datasets/1.1.0/V1_Breast_Cancer_Block_A_Section_1 ).
    图片.png

    run_UBC_DLPFC_data.py: The demo code for unsupervised batch correction for multiple DLPFC data.
    图片.png

生活很好,有你更好

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

推荐阅读更多精彩内容