hdWGCNA:单细胞WGCNA分析方法



0. 数据准备

  • 输入数据集的要求:已经进行了如下分析的Seurat对象
  • 导入演示数据
#官方演示数据集
wget https://swaruplab.bio.uci.edu/public_data/Zhou_2020.rds
seurat_obj <- readRDS('Zhou_2020.rds')

这是一个正常的脑组织数据集,包含了使用Harmony整合的12个样本的Seurat对象。

  • 加载R包和数据集
# single-cell analysis package
library(Seurat)

# plotting and data science packages
library(tidyverse)
library(cowplot)
library(patchwork)

# co-expression network analysis packages:
library(WGCNA)
library(hdWGCNA)

# using the cowplot theme for ggplot
theme_set(theme_cowplot())

# set random seed for reproducibility
set.seed(12345)
#可视化一下看看
p <- DimPlot(seurat_obj, group.by='cell_type', label=TRUE) +
    umap_theme()
p

1. Set up Seurat object for WGCNA

在进行hdWGCNA运算之前,我们首先需要处理Seurat object。绝大多数 hdWGCNA运算所需要的信息都储存在Seurat object’s @misc slot。一个Seurat object can hold multiple hdWGCNA experiments, for example representing different cell types in the same single-cell dataset. Notably, since we consider hdWGCNA to be a downstream data analysis step, we do not support subsetting the Seurat object after SetupForWGCNA has been run.

在这一步我们使用 SetupForWGCNA 来指定 hdWGNCA experiment的名称。这一步指定了用于WGCNA分析的基因 。关于基因的选择,有如下三种方法:

在这个演示里,我们选择了至少在5%的细胞中表达的基因(3857个),并将hdWGCNA experiment 命名为 “tutorial”。

seurat_obj <- SetupForWGCNA(
  seurat_obj,
  gene_select = "fraction", # the gene selection approach
  fraction = 0.05, # fraction of cells that a gene needs to be expressed in order to be included
  wgcna_name = "tutorial" # the name of the hdWGCNA experiment
)
length(seurat_obj@misc$tutorial$wgcna_genes)
# [1] 12639

2. Construct metacells

准备好 Seurat 对象后,首先需要从单细胞数据集中构建metacells。 简而言之,metacells是来自相同生物样品的一小组相似细胞的集合。 k-最近邻 (KNN) 算法可以识别相似细胞组,然后计算这些细胞的平均或总表达量,从而产生metacells基因表达矩阵。 metacells表达矩阵的稀疏性与原始表达矩阵相比大大降低,因此更适用于后续分析(WGCNA 等相关网络方法对数据稀疏性很敏感)。(单细胞表观基因组方法,例如 Cicero,在构建可访问网络之前采用了类似的元细胞聚合方法。)
hdWGCNA 使用 MetacellsByGroups函数来构造给定单细胞数据集的metacells表达矩阵。 此函数为存储在 hdWGCNA experiment中的metacells数据集构造一个新的 Seurat 对象。 group.by 参数确定将在哪些组中构建metacells。我们只想从来自相同生物样本的细胞构建metacells,因此通过 group.by 参数将该信息传递给 hdWGCNA 至关重要。 此外,我们通常为每种细胞类型分别构建metacells。 因此,在这个演示中,我们按 Sample 和 cell_type 进行分组以达到预期的结果。
k值(The number of cells to be aggregated)应根据输入数据集的大小进行调整,通常较小的 k 数可用于小数据集。 我们一般使用 20 到 75 之间的 k 值。本教程使用的数据集有 40,039 个细胞,每个生物样本中的细胞数从 890 到 8,188 不等,这里我们使用 k=25。 可以使用 max_shared 参数调整元单元之间允许的重叠量。

注意:metacells聚合方法对extremely underrepresented cell types效果不好。 例如,在这个数据集中,脑血管细胞(周细胞和内皮细胞)是the least represented,因此将它们排除在此分析之外。 MetacellsByGroups 有一个参数 min_cells 来排除小于指定单元数的组。

# construct metacells  in each group
seurat_obj <- MetacellsByGroups(
  seurat_obj = seurat_obj,
  group.by = c("cell_type", "Sample"), # specify the columns in seurat_obj@meta.data to group by
  k = 25, # nearest-neighbors parameter
  max_shared = 10, # maximum number of shared cells between two metacells
  ident.group = 'cell_type' # set the Idents of the metacell seurat object
)

# normalize metacell expression matrix:
seurat_obj <- NormalizeMetacells(seurat_obj)

Optional: Process the Metacell Seurat Object

由于我们将 Metacell 表达信息存储为单独的 Seurat 对象,因此我们可以对 Metacell 数据运行 Seurat 函数。 我们可以使用 GetMetacellObject 从 hdWGCNA 实验中获取元单元对象。

metacell_obj <- GetMetacellObject(seurat_obj)

Additionally, we have included a few wrapper functions to apply the Seurat workflow to the metacell object within the hdWGCNA experiment. Here we apply these wrapper functions to process the metacell object and visualize the aggregated expression profiles in two dimensions with UMAP.

seurat_obj <- NormalizeMetacells(seurat_obj)
seurat_obj <- ScaleMetacells(seurat_obj, features=VariableFeatures(seurat_obj))
seurat_obj <- RunPCAMetacells(seurat_obj, features=VariableFeatures(seurat_obj))
seurat_obj <- RunHarmonyMetacells(seurat_obj, group.by.vars='Sample')
seurat_obj <- RunUMAPMetacells(seurat_obj, reduction='harmony', dims=1:15)

这些结果储存在Seurat object@misc slot

p1 <- DimPlotMetacells(seurat_obj, group.by='cell_type') + umap_theme() + ggtitle("Cell Type")
p2 <- DimPlotMetacells(seurat_obj, group.by='Sample') + umap_theme() + ggtitle("Sample")
p1 | p2

3. 共表达网络分析

3.1 Set up the expression matrix

这一步将指定将用于网络分析的表达矩阵,使用SetDatExpr 函数,用于存储将用于下游网络分析的给定单元组的转置表达式矩阵。这一步默认情况下使用元单元表达式矩阵 (use_metacells=TRUE)。

由于我们只想分析抑制性神经元 (INH) ,因此使用group_name进行了指定。如果想包含多个细胞类型,可以用group_name = c("INH", "EX")

seurat_obj <- SetDatExpr(
  seurat_obj,
  group_name = "INH", # the name of the group of interest in the group.by column
  group.by='cell_type', # the metadata column containing the cell type info. This same column should have also been used in MetacellsByGroups
  assay = 'RNA', # using RNA assay
  slot = 'data' # using normalized data
)
3.2 挑选软阈值

接下来我们将选择“软阈值”。这是 hdWGNCA 分析(以及普通 WGCNA)中极其重要的一步。 hdWGCNA 构建基因-基因相关邻接矩阵来推断基因之间的共表达关系。将相关性提高到一个幂以减少相关矩阵中存在的噪声,从而保留强连接并去除弱连接。因此,确定软功率阈值的适当值至关重要。

我们使用 TestSoftPowers 函数来为不同的软阈值进行parameter sweep。此功能通过检查不同power值的结果网络拓扑结构,帮助我们指导我们选择构建共表达网络的软阈值。共表达网络应该具有无标度拓扑,因此 TestSoftPowers 函数模拟了共表达网络在不同软功率阈值下与无标度图的相似程度。此外,我们可以使用函数 PlotSoftPowers 来可视化parameter sweep的结果。

# Test different soft powers:
seurat_obj <- TestSoftPowers(
  seurat_obj,
  networkType = 'signed' # you can also use "unsigned" or "signed hybrid"
)

# plot the results:
plot_list <- PlotSoftPowers(seurat_obj)

# assemble with patchwork
wrap_plots(plot_list, ncol=2)

WGCNA 和 hdWGCNA 的一般要求是选择无标度拓扑模型拟合大于或等于 0.8 的最低软阈值,因此在这种情况下,我们选择9为软阈值。如果不设置软阈值,ConstructNetwork 将 自动选择软阈值。

Parameter sweep的输出结果存储在 hdWGCNA experiment中,可以使用 GetPowerTable 函数访问以进行进一步检查:

power_table <- GetPowerTable(seurat_obj)
head(power_table)
# Power    SFT.R.sq      slope truncated.R.sq   mean.k. median.k.    max.k.
# 1     1 0.192964167 10.1742363      0.9508966 5887.9290 5901.7852 6501.0283
# 2     2 0.001621749  0.4801145      0.9828730 3092.6379 3092.7878 3835.2202
# 3     3 0.150812707 -3.2345892      0.9918393 1657.6198 1646.6169 2349.6967
# 4     4 0.407513393 -4.5383131      0.9898440  905.8402  890.9431 1486.7695
# 5     5 0.585949838 -4.8303680      0.9917528  504.4288  490.2683  968.1157
# 6     6 0.677668481 -4.6749295      0.9935437  286.1585  274.3233  646.6741
3.3 构建共表达网络

随后我们就可以使用ConstructNetwork函数构建共表达网络了。(The parameters for blockwiseConsensusModules can be passed directly to ConstructNetwork with the same parameter names.)

# construct co-expression network:
seurat_obj <- ConstructNetwork(
  seurat_obj, soft_power=9,
  setDatExpr=FALSE,
  tom_name = 'INH' # name of the topoligical overlap matrix written to disk
)

共表达网络的可视化

PlotDendrogram(seurat_obj, main='INH hdWGCNA Dendrogram')

模块组成基因

seurat_obj@misc$tutorial$wgcna_modules %>% head
#             gene_name module color   kME_grey   kME_black  kME_green  kME_grey60    kME_red kME_lightcyan      kME_tan kME_royalblue  kME_brown kME_turquoise kME_lightyellow
# AL627309.1 AL627309.1   grey  grey 0.06467554  0.01169450 0.01346971 -0.03798701 0.06389461 -0.0234879179  0.024303002   0.038220179 0.04439502    0.06965388     0.032234528
# LINC01409   LINC01409  black black 0.09562911  0.14258271 0.04532450 -0.03838167 0.07259315  0.0088322467  0.089630925   0.090604216 0.05090351    0.08046802     0.020385146
# LINC01128   LINC01128   grey  grey 0.15025185  0.01456062 0.08874896  0.05263830 0.10848233  0.0233505152  0.036381392   0.003267192 0.03539629    0.13056570     0.034278458
# NOC2L           NOC2L   grey  grey 0.16386360 -0.03203850 0.03588840  0.02706705 0.09822920  0.0223949415  0.044098042   0.013722673 0.02077281    0.15718774     0.007470014
# AGRN             AGRN   grey  grey 0.06157292  0.05102518 0.08148903 -0.01092934 0.07061436 -0.0000888185 -0.011453627   0.028133157 0.06998927    0.04722421     0.043001474
# C1orf159     C1orf159   grey  grey 0.03297087  0.02592184 0.02742151 -0.01461799 0.01558220 -0.0035206063  0.001851281   0.010550009 0.01190812    0.02965008     0.025278788
#            kME_lightgreen kME_yellow kME_salmon    kME_pink kME_magenta kME_purple     kME_cyan kME_greenyellow kME_midnightblue     kME_blue kME_darkred kME_darkgreen
# AL627309.1     0.04998851 0.03999694 0.05183194 0.057363638  0.02197360 0.07381381  0.022132715    0.0309476545      0.020243486  0.039376677  0.04946240  -0.041431912
# LINC01409      0.06062908 0.12675277 0.06172935 0.086179445  0.03545920 0.10836689 -0.010795315    0.1022885505     -0.013422943  0.008179931  0.12100397  -0.005445180
# LINC01128      0.12536877 0.04878530 0.09510907 0.037100032  0.06866021 0.07657537  0.108189685    0.0433716590      0.045969133  0.112251246  0.08318723  -0.045492473
# NOC2L          0.07886350 0.01475673 0.13338651 0.049125991  0.07227435 0.07841803  0.157399769    0.0232836501      0.025693226  0.166539003  0.03711394  -0.050635030
# AGRN           0.05589207 0.11086059 0.02339966 0.032878117 -0.03134549 0.05606373  0.036350264    0.0422577304      0.009418002  0.034792567  0.01463038   0.012060664
# C1orf159       0.01596504 0.03185948 0.02135547 0.004147534  0.00511424 0.02773708 -0.003693451    0.0004414052      0.005552397 -0.002849447  0.03548955   0.002986991

table(seurat_obj@misc$tutorial$wgcna_modules$module)

#         grey        black        green       grey60          red    lightcyan          tan    royalblue        brown    turquoise  lightyellow   lightgreen 
#         9597          154          188           87          166           88          129           67          230          281           76           86 
#       yellow       salmon         pink      magenta       purple         cyan  greenyellow midnightblue         blue      darkred    darkgreen 
#          203          123          151          150          149          123          131          112          232           60           56
3.4 Optional: inspect the topoligcal overlap matrix (TOM)

hdWGCNA represents the co-expression network as a topoligcal overlap matrix (TOM). This is a square matrix of genes by genes, where each value is the topoligcal overlap between the genes. The TOM is written to the disk when running ConstructNetwork, and we can load it into R using the GetTOM function. Advanced users may wish to inspect the TOM for custom downstream analyses.

TOM <- GetTOM(seurat_obj)

4. Module Eigengenes and Connectivity

在这一部分,我们将介绍如何计算单个细胞中的模块特征基因,以及如何计算每个基因的基于特征基因的连通性。

4.1 Compute harmonized module eigengenes

Module Eigengenes (MEs) 是一种常用的指标,用于总结整个共表达模块的基因表达谱。简而言之,模块特征基因是通过对包含每个模块的基因表达矩阵的子集执行主成分分析 (PCA) 来计算的。这些 PCA 矩阵中的第一个PC 就是 ME。

hdWGCNA 包含一个函数 ModuleEigengenes 来计算单个单元格中的模块特征基因。此外,我们允许用户对 ME 应用 Harmony 批量校正,从而产生协调的模块特征基因 (hME)。以下代码使用 group.by.vars 参数执行由原始样本协调的模块特征基因计算。

# need to run ScaleData first or else harmony throws an error:
seurat_obj <- ScaleData(seurat_obj, features=VariableFeatures(seurat_obj))

# compute all MEs in the full single-cell dataset
seurat_obj <- ModuleEigengenes(
 seurat_obj,
 group.by.vars="Sample"
) #时间较长

ME矩阵存储为一个矩阵,其中每一行是一个细胞,每一列是一个模块。 可以使用 GetMEs 函数从 Seurat 对象中提取此矩阵,该函数默认检索 hME。

# harmonized module eigengenes:
hMEs <- GetMEs(seurat_obj) #每个细胞对于每个模块的特征值
head(hMEs)
#                     lightcyan darkgreen       red turquoise   salmon   purple    black  darkred lightgreen royalblue      tan greenyellow
# AGTCTTTGTTGATTCG-11 -3.728266 -2.410602 3.8018569  5.540111 6.002367 4.829816 5.308432 4.066452   3.087318  5.677753 6.637234    7.003169
# AGCAGCCTCCAGATCA-11 -3.242143 -2.725894 3.1256861  5.243863 5.191766 4.053458 5.157324 3.919190   4.276658  5.821838 6.023351    7.905414
# CTTAGGATCTCATTCA-11 -3.697165 -1.910883 4.2795263  5.803626 5.108230 3.798885 5.144346 3.710556   3.451550  5.056269 6.701767    8.761117
# AACTCAGAGGAATGGA-11 -3.438431 -3.141482 0.9050427  2.689078 2.172585 2.790341 7.279000 4.412214   2.103700  4.536073 5.719890    6.730259
# ACGAGGAGTCTCCACT-11 -3.752210 -2.695012 3.8369253  1.600519 2.965712 3.075344 5.888933 2.554901   2.433528  5.539103 8.523871    7.847494
# ACTTACTTCGGAGCAA-11 -3.193025 -1.252224 3.3207978  4.391861 4.744844 2.290026 4.946544 2.965951   4.713082  4.360662 6.201866    6.620763
#                     midnightblue     pink  magenta    brown lightyellow   yellow     grey    green   grey60       cyan       blue
# AGTCTTTGTTGATTCG-11     6.576508 5.369716 6.578784 7.616044    4.476879 7.459703 37.92130 7.224577 4.975482 13.5867856 14.6903483
# AGCAGCCTCCAGATCA-11     5.517671 4.143108 5.390045 8.251945    4.413418 7.395514 38.43114 8.578032 6.748117 13.6829303 12.4873768
# CTTAGGATCTCATTCA-11     5.156724 5.510369 4.768236 8.243442    4.932943 8.103645 36.48388 7.867925 7.018552 12.0219104 10.7403398
# AACTCAGAGGAATGGA-11     5.603454 6.998764 6.598624 9.449461    6.600172 7.963021 29.11238 8.023318 1.401920  0.4015516  0.9184507
# ACGAGGAGTCTCCACT-11     7.024574 7.239748 7.934361 9.010840    5.841042 8.041680 32.85802 7.300763 4.790319  4.6499158  3.8687891
# ACTTACTTCGGAGCAA-11     5.045800 4.115106 5.537796 8.347864    6.416239 6.914440 30.72694 7.490030 5.289845  9.2345490 10.2069815
# module eigengenes:
MEs <- GetMEs(seurat_obj, harmonized=FALSE)
4.2 Compute module connectivity

在共表达网络分析中,我们经常希望关注“hub gene”,即在每个模块内高度连接的那些基因。 因此,我们希望确定每个基因的eigengene-based connectivity,也称为 kME。 hdWGCNA 使用 ModuleConnectivity 函数来计算完整单细胞数据集中的 kME 值,而不是元细胞数据集。 该函数本质上计算基因和模块特征基因之间的成对相关性。 kME可被用于该数据集中的所有细胞,但我们建议在用于计算 ConstructNetwork 的细胞类型或分组中计算 kME。

# compute eigengene-based connectivity (kME):
seurat_obj <- ModuleConnectivity(
  seurat_obj,
  group.by = 'cell_type', group_name = 'INH'
)

We can visualize the genes in each module ranked by kME using the PlotKMEs function.

# plot genes ranked by kME for each module
p <- PlotKMEs(seurat_obj, ncol=5)

p
4.3 Getting the module assignment table

hdWGCNA 允许使用 GetModules 函数轻松访问module assignment table。 该表由三列组成:gene_name: 存储基因的符号或 ID,module: 存储基因的模块分配,color: 存储每个模块的颜色映射,用于许多下游绘图步骤。 如果在此 hdWGCNA 实验中调用了 ModuleConnectivity,则此表将具有每个模块的 kME 的附加列。

# get the module assignment table:
modules <- GetModules(seurat_obj)

# show the first 6 columns:
head(modules[,1:6])
# gene_name  module color   kME_INH-M1   kME_INH-M2    kME_INH-M3
# AL627309.1 AL627309.1    grey  grey -0.032349090  0.029917426  0.0379323320
# LINC01409   LINC01409 INH-M19 brown -0.045140924 -0.013473381  0.0102194168
# LINC01128   LINC01128    grey  grey  0.050793988  0.109578251  0.1153173093
# NOC2L           NOC2L    grey  grey  0.032490535  0.164524557  0.1699131451
# AGRN             AGRN INH-M19 brown -0.008488577  0.035558532  0.0379326966
# C1orf159     C1orf159    grey  grey -0.015618737  0.002235229 -0.0003824554

鉴定模块内hub基因

# get hub genes
hub_df <- GetHubGenes(seurat_obj, n_hubs = 10)

head(hub_df)
# gene_name module       kME
# 1   SEPTIN2 INH-M1 0.3019739
# 2     MICU2 INH-M1 0.3140922
# 3     WDR37 INH-M1 0.3204428
# 4   UBE2Q2L INH-M1 0.3297472
# 5   HSD17B4 INH-M1 0.3399253
# 6      ASB3 INH-M1 0.3417400

This wraps up the critical analysis steps for hdWGCNA, so remember to save your output.

saveRDS(seurat_obj, file='hdWGCNA_object.rds')
  • 计算每个细胞对于每个模块hub基因的表达活性(module score)
    可使用seurat的AddModuleScore 函数或者Ucell算法
# compute gene scoring for the top 25 hub genes by kME for each module
# with Seurat method
seurat_obj <- ModuleExprScore(
  seurat_obj,
  n_genes = 25,
  method='Seurat'
)

# compute gene scoring for the top 25 hub genes by kME for each module
# with UCell method
library(UCell)
seurat_obj <- ModuleExprScore(
  seurat_obj,
  n_genes = 25,
  method='UCell'
)

5. Basic Visualization

5.1 Module Feature Plots

hdWGCNA includes the ModuleFeaturePlot function to consruct FeaturePlots for each co-expression module colored by each module’s uniquely assigned color.

# make a featureplot of hMEs for each module
plot_list <- ModuleFeaturePlot(
  seurat_obj,
  features='hMEs', # plot the hMEs
  order=TRUE # order so the points with highest hMEs are on top
)

# stitch together with patchwork
wrap_plots(plot_list, ncol=6)

We can also plot the hub gene signature score using the same function:

# make a featureplot of hub scores for each module
plot_list <- ModuleFeaturePlot(
  seurat_obj,
  features='scores', # plot the hub gene scores
  order='shuffle', # order so cells are shuffled
  ucell = TRUE # depending on Seurat vs UCell for gene scoring
)

# stitch together with patchwork
wrap_plots(plot_list, ncol=6)
5.2 Module Correlations

hdWGCNA includes the ModuleCorrelogram function to visualize the correlation between each module based on their hMEs, MEs, or hub gene scores using the R package corrplot.

ModuleCorrelogram(seurat_obj)
5.3 Seurat plotting functions

The base Seurat plotting functions are also great for visualizing hdWGCNA outputs. Here we demonstrate plotting hMEs using DotPlot and VlnPlot. The key to using Seurat’s plotting functions to visualize the hdWGCNA data is to add it into the Seurat object’s @meta.data slot:

# get hMEs from seurat object
MEs <- GetMEs(seurat_obj, harmonized=TRUE)
mods <- colnames(MEs); mods <- mods[mods != 'grey']

# add hMEs to Seurat meta-data:
seurat_obj@meta.data <- cbind(seurat_obj@meta.data, MEs)

Now we can easily use Seurat’s DotPlot function:

# plot with Seurat's DotPlot function
p <- DotPlot(seurat_obj, features=mods, group.by = 'cell_type')

# flip the x/y axes, rotate the axis labels, and change color scheme:
p <- p +
  coord_flip() +
  RotatedAxis() +
  scale_color_gradient2(high='red', mid='grey95', low='blue')

# plot output
p
# plot with Seurat's DotPlot function
p <- DotPlot(seurat_obj, features=mods, group.by = 'cell_type')

# flip the x/y axes, rotate the axis labels, and change color scheme:
p <- p +
  coord_flip() +
  RotatedAxis() +
  scale_color_gradient2(high='red', mid='grey95', low='blue')

# plot output
p

6. 相关文献


参考:https://smorabit.github.io/hdWGCNA/articles/basic_tutorial.html

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

推荐阅读更多精彩内容