ArchR官网教程学习笔记13:ChromVAR偏差富集

如前几章所示，TF motif富集可以帮助我们预测在我们感兴趣的细胞类型中哪些转录因子最活跃。然而，这些富集并不以每个细胞为基础计算，也没有考虑到Tn5转座酶的插入序列偏倚。chromVAR，是一个Greenlead实验室构建的R包，就是为了解决这些问题而创建的。chromVAR是设计用于在稀疏染色质可接近数据里预测富集的转录因子活性（在每个细胞的基础上）。chromVAR的两个主要输出是:

1.“偏差”：偏差是对给定特征(即motif)的每个细胞的可接近性与基于所有细胞或样本平均值的预期可接近性的偏差校正测量。

2.“z分数”：也称为“偏差分数”，是所有细胞中每个偏差校正后的z分数。偏差分数的绝对值与每个细胞的读取深度相关。这是因为，读得越多，你就越有信心认为给定特征(即motif)的每个细胞可接近性与预期的差异大于偶然发生的差异。

chromVAR的主要限制之一是，它是在scATAC-seq数据生成早期设计的，当时的实验只有几百个细胞。在这个实验尺度下，chromVAR可以很容易地将整个“细胞－峰”矩阵读取到内存中，以快速计算TF偏差。然而，目前的实验方法使用成千上万的细胞，产生“细胞－峰”矩阵，这很难读取到内存中。即使对于50,000个细胞的中等大小的数据集，这也会导致运行时间和内存使用的显著增加。

为了规避这些限制，ArchR通过独立分析样本子矩阵（sub-matrices）实现了相同的chromVAR分析工作流程。

首先，ArchR读取每个子样本中所有细胞的每个峰的全局可接近性。其次，对于每个峰，ArchR识别出一组与GC含量和可接近性相匹配的背景峰。第三，ArchR使用这一背景峰集和全局可接近性，独立地使用chromVAR计算每个样本的偏置校正偏差。实现在任何给定时间只需要将来自5,000-10,000个细胞的数据加载到内存中，从而最小化内存需求，能够进行可扩展分析，并提高运行时间性能。

（一）Motif偏差

首先，确定我们在ArchRProject中添加了motif注释。

> if("Motif" %ni% names(projHeme5@peakAnnotation)){
    projHeme5 <- addMotifAnnotations(ArchRProj = projHeme5, motifSet = "cisbp", name = "Motif")
}

我们还需要添加一组用于计算偏差的背景峰。使用chromVAR::getBackgroundPeaks()函数来选择背景峰，该函数基于GC含量的相似性和使用马氏距离在所有样本中的片段数量来选择峰。

> projHeme5 <- addBgdPeaks(projHeme5)
Identifying Background Peaks!

现在我们可以使用addDeviationsMatrix()函数计算在所有motif注释里每个细胞的偏差。这个函数有一个名为matrixName的可选参数，它允许我们定义将存储在Arrow files中的偏差矩阵的名称。如果我们不为该参数提供值，如下面的示例所示，则该函数通过在peakAnnotation的名称中添加单词“matrix”来创建一个矩阵名称。下面的例子在每个Arrow files中创建了一个偏差矩阵，称为“MotifMatrix”:

> projHeme5 <- addDeviationsMatrix(
  ArchRProj = projHeme5, 
  peakAnnotation = "Motif",
  force = TRUE
)

要访问这些偏差，我们使用getVarDeviations() 函数。如果我们希望这个函数返回一个ggplot对象，我们设置plot = TRUE，否则，这个函数将返回DataFrame对象。该DataFrame对象的前几行在函数运行时默认显示：

> plotVarDev <- getVarDeviations(projHeme5, name = "MotifMatrix", plot = TRUE)
DataFrame with 6 rows and 6 columns
     seqnames     idx        name     combinedVars       combinedMeans      rank
        <Rle> <array>     <array>        <numeric>           <numeric> <integer>
f388        z     388   GATA2_388 11.8942823585118  -0.037174877691045         1
f336        z     336    SPIB_336 11.6325682721513 -0.0831267955064216         2
f383        z     383   GATA1_383 11.6253116037962 -0.0402170665586582         3
f155        z     155   CEBPA_155 11.3334664092501  -0.179614402768299         4
f385        z     385   GATA5_385 10.3039978837462 -0.0327070456187338         5
f651        z     651 SMARCC1_651 10.1573838523962  -0.122959684614117         6

从上面的数据中，你可以看到MotifMatrix的seqname不是染色体。通常，在TileMatrix、PeakMatrix和GeneScoreMatrix等矩阵中，我们将染色体信息存储在seqnames中。MotifMatrix没有对应的位置信息，但相反，从chromVAR里来的“devations”和“z-scores”chromVAR到相同的矩阵使用两种不同的seqnames——偏差和z。在之后的分析中很重要，如果你尝试使用MotifMatrix(类Sparse.Assays.Matrix)等getMarkerFeatures()函数，在这些类型的操作中，ArchR会期望你将MotifMatrix子集划分为两个seqname中的一个(即选择z-scores或偏差来执行计算)。

然后我们可以画出这些可变偏差：

> plotVarDev
> plotPDF(plotVarDev, name = "Variable-Motif-Deviation-Scores", width = 5, height = 5, ArchRProj = projHeme5, addDOC = FALSE)

如果我们想要为下游分析提取一个motif子集该怎么办呢?我们可以使用getFeatures()函数来实现这一点。下面的paste(motifs, collapse="|")语句创建了一个连接的or语句，该语句允许选择所有的motifs:

> motifs <- c("GATA1", "CEBPA", "EBF1", "IRF4", "TBX21", "PAX5")
> markerMotifs <- getFeatures(projHeme5, select = paste(motifs, collapse="|"), useMatrix = "MotifMatrix")
> markerMotifs
 [1] "z:TBX21_780"          "z:PAX5_709"           "z:IRF4_632"          
 [4] "z:GATA1_383"          "z:CEBPA_155"          "z:EBF1_67"           
 [7] "z:SREBF1_22"          "deviations:TBX21_780" "deviations:PAX5_709" 
[10] "deviations:IRF4_632"  "deviations:GATA1_383" "deviations:CEBPA_155"
[13] "deviations:EBF1_67"   "deviations:SREBF1_22"

如上所述，MotifMatrix包含了z-score和偏差的seqname，如上面的“z:”和“deviations:”所示。要获得与z-scores对应的特性，我们可以使用grep。不幸的是，在上面的例子中，你可以看到除了“EBF1”之外，我们还选择了我们不想分析的“SREBF1”。因此，我们使用%ni%表达式将其删除，该表达式是ArchR辅助函数，提供了以R为基础的%in%的反义词。

> markerMotifs <- grep("z:", markerMotifs, value = TRUE)
> markerMotifs <- markerMotifs[markerMotifs %ni% "z:SREBF1_22"]
> markerMotifs
[1] "z:TBX21_780" "z:PAX5_709"  "z:IRF4_632"  "z:GATA1_383" "z:CEBPA_155"
[6] "z:EBF1_67"

现在我们已经有了我们感兴趣的特征的名称，我们可以绘制每个cluster的chromVAR偏差分数的分布。请注意，我们提供了之前在基因评分分析中计算的impute weights。提醒一下，这些权重允许我们使附近的细胞信号变得平滑，这对我们稀疏的scATAC-seq数据的背景是有帮助的。

> p <- plotGroups(ArchRProj = projHeme5, 
  groupBy = "Clusters2", 
  colorBy = "MotifMatrix", 
  name = markerMotifs,
  imputeWeights = getImputeWeights(projHeme5)
)

在单个图形中绘制所有motif的分布：

> p2 <- lapply(seq_along(p), function(x){
  if(x != 1){
    p[[x]] + guides(color = FALSE, fill = FALSE) + 
    theme_ArchR(baseSize = 6) +
    theme(plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "cm")) +
    theme(
        axis.text.y=element_blank(), 
        axis.ticks.y=element_blank(),
        axis.title.y=element_blank()
    ) + ylab("")
  }else{
    p[[x]] + guides(color = FALSE, fill = FALSE) + 
    theme_ArchR(baseSize = 6) +
    theme(plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "cm")) +
    theme(
        axis.ticks.y=element_blank(),
        axis.title.y=element_blank()
    ) + ylab("")
  }
})
> do.call(cowplot::plot_grid, c(list(nrow = 1, rel_widths = c(2, rep(1, length(p2) - 1))),p2))

> plotPDF(p, name = "Plot-Groups-Deviations-w-Imputation", width = 5, height = 5, ArchRProj = projHeme5, addDOC = FALSE)

除了看这些z-score的分布，我们还可以将z-score覆盖在UMAP中，就像我们之前对基因评分所做的那样：

> p <- plotEmbedding(
    ArchRProj = projHeme5, 
    colorBy = "MotifMatrix", 
    name = sort(markerMotifs), 
    embedding = "UMAP",
    imputeWeights = getImputeWeights(projHeme5)
)

画这些motif的UMAP图：

> p2 <- lapply(p, function(x){
    x + guides(color = FALSE, fill = FALSE) + 
    theme_ArchR(baseSize = 6.5) +
    theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) +
    theme(
        axis.text.x=element_blank(), 
        axis.ticks.x=element_blank(), 
        axis.text.y=element_blank(), 
        axis.ticks.y=element_blank()
    )
})
> do.call(cowplot::plot_grid, c(list(ncol = 3),p2))

为了了解这些TF偏差z-scores与推断的相应TF基因的基因评分之间的关系，我们可以将每个TFs的基因评分叠加在UMAP图上：

> markerRNA <- getFeatures(projHeme5, select = paste(motifs, collapse="|"), useMatrix = "GeneScoreMatrix")
> markerRNA <- markerRNA[markerRNA %ni% c("SREBF1","CEBPA-DT")]
> markerRNA
[1] "TBX21" "CEBPA" "EBF1"  "IRF4"  "PAX5"  "GATA1"

> p <- plotEmbedding(
    ArchRProj = projHeme5, 
    colorBy = "GeneScoreMatrix", 
    name = sort(markerRNA), 
    embedding = "UMAP",
    imputeWeights = getImputeWeights(projHeme5)
)

> p2 <- lapply(p, function(x){
    x + guides(color = FALSE, fill = FALSE) + 
    theme_ArchR(baseSize = 6.5) +
    theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) +
    theme(
        axis.text.x=element_blank(), 
        axis.ticks.x=element_blank(), 
        axis.text.y=element_blank(), 
        axis.ticks.y=element_blank()
    )
})
> do.call(cowplot::plot_grid, c(list(ncol = 3),p2))

类似地，因为我们之前将scATAC-seq数据与相应的scRNA-seq数据相连接，所以我们可以在UMAP中绘制每个TFs的连接的基因表达：

> markerRNA <- getFeatures(projHeme5, select = paste(motifs, collapse="|"), useMatrix = "GeneIntegrationMatrix")
> markerRNA <- markerRNA[markerRNA %ni% c("SREBF1","CEBPA-DT")]
> markerRNA
[1] "TBX21" "CEBPA" "EBF1"  "IRF4"  "PAX5"  "GATA1"

> p <- plotEmbedding(
    ArchRProj = projHeme5, 
    colorBy = "GeneIntegrationMatrix", 
    name = sort(markerRNA), 
    embedding = "UMAP",
    continuousSet = "blueYellow",
    imputeWeights = getImputeWeights(projHeme5)
)

> p2 <- lapply(p, function(x){
    x + guides(color = FALSE, fill = FALSE) + 
    theme_ArchR(baseSize = 6.5) +
    theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) +
    theme(
        axis.text.x=element_blank(), 
        axis.ticks.x=element_blank(), 
        axis.text.y=element_blank(), 
        axis.ticks.y=element_blank()
    )
})
> do.call(cowplot::plot_grid, c(list(ncol = 3),p2))

（二）自定义偏差

在峰注释富集一章中，我们介绍了如何为任意一组基因组区域创建峰注释。这包括(i) ArchR支持的区域集，例如从ENCODE和从bulk ATAC-seq的峰集挑选的TF结合位点，以及(ii)用户提供的自定义区域集。如果你还没有阅读本节，我们建议你阅读一下，以便更好地理解peak annotations是如何工作的。

这些峰注释可以像motif一样用于偏差计算。这里我们提供了如何运行这些分析的示例，但请注意，下游分析与上一节关于motifs的分析相同，因此我们不提供关于每一步代码的详细信息。在Arrow files中创建了偏差矩阵之后，其余的操作都是相同的。

（１）Encode TFBS

如果你还没有为“EncodeTFBS”区域集添加注释矩阵，现在让我们做一下：

> if("EncodeTFBS" %ni% names(projHeme5@peakAnnotation)){
    projHeme5 <- addArchRAnnotations(ArchRProj = projHeme5, collection = "EncodeTFBS")
}

然后，我们创建一个偏差矩阵，将这个峰注释提供给peakAnnotation参数:

> projHeme5 <- addDeviationsMatrix(
  ArchRProj = projHeme5, 
  peakAnnotation = "EncodeTFBS",
  force = TRUE
)

画一个偏差排序的图：

> plotVarDev <- getVarDeviations(projHeme5, plot = TRUE, name = "EncodeTFBSMatrix")
DataFrame with 6 rows and 6 columns
     seqnames     idx                    name     combinedVars        combinedMeans      rank
        <Rle> <array>                 <array>        <numeric>            <numeric> <integer>
f222        z     222     222.GATA2_S-K562... 13.6462717266687  -0.0347371238702175         1
f41         z      41      41.EZH2_39-NHEK... 12.1902835430786    0.111737526644295         2
f584        z     584 584.GATA_1-PBDEFetal... 11.4020981531359 -0.00726058996329407         3
f542        z     542     542.TAL1_SC-K562... 11.3641140640469  -0.0178812008386299         4
f498        z     498      498.GATA_2-K562... 10.3108802557775  -0.0304219839395099         5
f44         z      44      44.EZH2_39-NHLF... 8.98198697694539   0.0672812498713452         6
> plotVarDev
> plotPDF(plotVarDev, name = "Variable-EncodeTFBS-Deviation-Scores", width = 5, height = 5, ArchRProj = projHeme5, addDOC = FALSE)

或者我们可以将这些TF结合位点划分为我们感兴趣的特定motif，然后在我们的UMAP嵌入中绘制每个细胞的偏差z-scores。

> tfs <- c("GATA_1", "CEBPB", "EBF1", "IRF4", "TBX21", "PAX5")
> getFeatures(projHeme5, select = paste(tfs, collapse="|"), useMatrix = "EncodeTFBSMatrix")
 [1] "z:584.GATA_1-PBDEFetal..."         
 [2] "z:582.GATA_1-PBDE..."              
 [3] "z:497.GATA_1-K562..."              
......
[39] "deviations:87.EBF1_SC-GM12878..."  
[40] "deviations:86.CEBPB_S-GM12878..."

> markerTFs <- getFeatures(projHeme5, select = paste(tfs, collapse="|"), useMatrix = "EncodeTFBSMatrix")
> markerTFs <- sort(grep("z:", markerTFs, value = TRUE))
> TFnames <- stringr::str_split(stringr::str_split(markerTFs, pattern = "\\.", simplify=TRUE)[,2], pattern = "-", simplify = TRUE)[,1]
> markerTFs <- markerTFs[!duplicated(TFnames)]
> markerTFs
[1] "z:101.PAX5_C2-GM12878..." "z:102.PAX5_N1-GM12878..."
[3] "z:173.CEBPB_S-HepG2..."   "z:278.CEBPB-A549..."     
[5] "z:293.EBF1_SC-GM12878..." "z:497.GATA_1-K562..."    
[7] "z:93.IRF4_SC-GM12878..."

> p <- plotEmbedding(
    ArchRProj = projHeme5, 
    colorBy = "EncodeTFBSMatrix", 
    name = markerTFs, 
    embedding = "UMAP",
    imputeWeights = getImputeWeights(projHeme5)
)

> p2 <- lapply(p, function(x){
    x + guides(color = FALSE, fill = FALSE) + 
    theme_ArchR(baseSize = 6.5) +
    theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) +
    theme(
        axis.text.x=element_blank(), 
        axis.ticks.x=element_blank(), 
        axis.text.y=element_blank(), 
        axis.ticks.y=element_blank()
    )
})
> do.call(cowplot::plot_grid, c(list(ncol = 3),p2))

（２）Bulk ATAC-seq

同样，我们可以使用ArchR-curated bulk ATAC-seq峰集进行motif偏差计算。

> if("ATAC" %ni% names(projHeme5@peakAnnotation)){
    projHeme5 <- addArchRAnnotations(ArchRProj = projHeme5, collection = "ATAC")
}

然后，我们创建一个偏差矩阵，将这个峰注释提供给peakAnnotation参数。

> projHeme5 <- addDeviationsMatrix(
  ArchRProj = projHeme5, 
  peakAnnotation = "ATAC",
  force = TRUE
)

> plotVarDev <- getVarDeviations(projHeme5, plot = TRUE, name = "ATACMatrix")
DataFrame with 6 rows and 6 columns
    seqnames     idx                  name     combinedVars       combinedMeans      rank
       <Rle> <array>               <array>        <numeric>           <numeric> <integer>
f22        z      22 IAtlas_T_CD8posCenMem 12.8801483758475 -0.0929263764838128         1
f86        z      86              Heme_CD8 12.5648193703514 -0.0758332391096771         2
f85        z      85              Heme_CD4  12.371741410312 -0.0543437044277343         3
f23        z      23 IAtlas_T_CD8posEffMem 12.3442474876336 -0.0977163360337115         4
f33        z      33 IAtlas_T_Th1Precursor  12.231771365977 -0.0834317302855867         5
f21        z      21       IAtlas_T_CD8pos 12.1835121500956  -0.090039542520433         6
> plotVarDev
> plotPDF(plotVarDev, name = "Variable-ATAC-Deviation-Scores", width = 5, height = 5, ArchRProj = projHeme5, addDOC = FALSE)

我们也可以绘制每个细胞里这些峰集的偏差z-score（在UMAP里）：

> ATACPeaks <- c("Heme_HSC", "Heme_LMPP", "Heme_Ery", "Heme_Mono", "Heme_CD4", "Heme_CD8", "Heme_B", "Heme_NK", "IAtlas_DC_Plasmacytoid")
> getFeatures(projHeme5, select = paste(ATACPeaks, collapse="|"), useMatrix = "ATACMatrix")
[1] "z:Heme_NK"                        
 [2] "z:Heme_Mono"                      
 [3] "z:Heme_LMPP"                      
......            
[17] "deviations:Heme_B"                
[18] "deviations:IAtlas_DC_Plasmacytoid"

> markerATAC <- getFeatures(projHeme5, select = paste(ATACPeaks, collapse="|"), useMatrix = "ATACMatrix")
> markerATAC <- sort(grep("z:", markerATAC, value = TRUE))
> markerATAC
[1] "z:Heme_B"                 "z:Heme_CD4"              
[3] "z:Heme_CD8"               "z:Heme_Ery"              
[5] "z:Heme_HSC"               "z:Heme_LMPP"             
[7] "z:Heme_Mono"              "z:Heme_NK"               
[9] "z:IAtlas_DC_Plasmacytoid"

> p <- plotEmbedding(
    ArchRProj = projHeme5, 
    colorBy = "ATACMatrix", 
    name = markerATAC, 
    embedding = "UMAP",
    imputeWeights = getImputeWeights(projHeme5)
)

> p2 <- lapply(p, function(x){
    x + guides(color = FALSE, fill = FALSE) + 
    theme_ArchR(baseSize = 6.5) +
    theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) +
    theme(
        axis.text.x=element_blank(), 
        axis.ticks.x=element_blank(), 
        axis.text.y=element_blank(), 
        axis.ticks.y=element_blank()
    )
})
> do.call(cowplot::plot_grid, c(list(ncol = 3),p2))

（３）自定义偏差

除了使用上面描述的ArchR管理的区域集，我们可以提供自己的自定义区域集作为peak注释。这些自定义注释的使用方式与ArchR管理的注释完全相同。

首先，如果您在前一章中还没有创建这个“EncodePeaks”注释，可以通过下载一些ENCODE峰集并调用addPeakAnnotations()来创建它。

> EncodePeaks <- c(
  Encode_K562_GATA1 = "https://www.encodeproject.org/files/ENCFF632NQI/@@download/ENCFF632NQI.bed.gz",
  Encode_GM12878_CEBPB = "https://www.encodeproject.org/files/ENCFF761MGJ/@@download/ENCFF761MGJ.bed.gz",
  Encode_K562_Ebf1 = "https://www.encodeproject.org/files/ENCFF868VSY/@@download/ENCFF868VSY.bed.gz",
  Encode_K562_Pax5 = "https://www.encodeproject.org/files/ENCFF339KUO/@@download/ENCFF339KUO.bed.gz"
)
> if("ChIP" %ni% names(projHeme5@peakAnnotation)){
    projHeme5 <- addPeakAnnotations(ArchRProj = projHeme5, regions = EncodePeaks, name = "ChIP")
}

创建偏差矩阵：

> projHeme5 <- addDeviationsMatrix(
  ArchRProj = projHeme5, 
  peakAnnotation = "ChIP",
  force = TRUE
)

分析工作流程的其余部分与前面多次提到的内容相同。
画出排列好的偏差：

> plotVarDev <- getVarDeviations(projHeme5, plot = TRUE, name = "ChIPMatrix")
DataFrame with 4 rows and 6 columns
   seqnames     idx                 name      combinedVars        combinedMeans      rank
      <Rle> <array>              <array>         <numeric>            <numeric> <integer>
f1        z       1    Encode_K562_GATA1  6.35199363570919  -0.0161048051589612         1
f3        z       3     Encode_K562_Ebf1  3.43122738670321   0.0364557074214317         2
f4        z       4     Encode_K562_Pax5  3.16880877416097 -0.00697454330040878         3
f2        z       2 Encode_GM12878_CEBPB 0.806426078308586  0.00560392567474442         4
> plotVarDev
> plotPDF(plotVarDev, name = "Variable-ChIP-Deviation-Scores", width = 5, height = 5, ArchRProj = projHeme5, addDOC = FALSE)

> plotPDF(plotVarDev, name = "Variable-ChIP-Deviation-Scores", width = 5, height = 5, ArchRProj = projHeme5, addDOC = FALSE)

在UMAP上绘制偏差z-score：

> markerChIP <- getFeatures(projHeme5, useMatrix = "ChIPMatrix")
> markerChIP <- sort(grep("z:", markerChIP, value = TRUE))
> markerChIP
[1] "z:Encode_GM12878_CEBPB" "z:Encode_K562_Ebf1"    
[3] "z:Encode_K562_GATA1"    "z:Encode_K562_Pax5" 

> p <- plotEmbedding(
    ArchRProj = projHeme5, 
    colorBy = "ChIPMatrix", 
    name = markerChIP, 
    embedding = "UMAP",
    imputeWeights = getImputeWeights(projHeme5)
)

> p2 <- lapply(p, function(x){
    x + guides(color = FALSE, fill = FALSE) + 
    theme_ArchR(baseSize = 6.5) +
    theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) +
    theme(
        axis.text.x=element_blank(), 
        axis.ticks.x=element_blank(), 
        axis.text.y=element_blank(), 
        axis.ticks.y=element_blank()
    )
})
> do.call(cowplot::plot_grid, c(list(ncol = 2),p2))

最后编辑于：2020.12.01 08:52:58

禁止转载，如需转载请通过简信或评论联系作者。

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 206,723评论 6赞 481
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 88,485评论 2赞 382
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 152,998评论 0赞 344
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 55,323评论 1赞 279
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 64,355评论 5赞 374
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 49,079评论 1赞 285
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 38,389评论 3赞 400
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 37,019评论 0赞 259
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 43,519评论 1赞 300
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 35,971评论 2赞 325
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 38,100评论 1赞 333
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 33,738评论 4赞 324
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 39,293评论 3赞 307
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 30,289评论 0赞 19
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 31,517评论 1赞 262
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 45,547评论 2赞 354
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 42,834评论 2赞 345