hdWGCNA：單細胞WGCNA分析方法

0. 數(shù)據準備

  • 輸入數(shù)據集的要求：已經進行了如下分析的Seurat對象
  • 導入演示數(shù)據
#官方演示數(shù)據集
wget https://swaruplab.bio.uci.edu/public_data/Zhou_2020.rds
seurat_obj <- readRDS('Zhou_2020.rds')

這是一個正常的腦組織數(shù)據集，包含了使用Harmony整合的12個樣本的Seurat對象日川。

  • 加載R包和數(shù)據集
# 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赖条。絕大多數(shù) 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 對象后，首先需要從單細胞數(shù)據集中構建metacells化戳。 簡而言之单料，metacells是來自相同生物樣品的一小組相似細胞的集合。 k-最近鄰 (KNN) 算法可以識別相似細胞組点楼，然后計算這些細胞的平均或總表達量扫尖，從而產生metacells基因表達矩陣。 metacells表達矩陣的稀疏性與原始表達矩陣相比大大降低掠廓，因此更適用于后續(xù)分析（WGCNA 等相關網絡方法對數(shù)據稀疏性很敏感）换怖。（單細胞表觀基因組方法，例如 Cicero蟀瞧，在構建可訪問網絡之前采用了類似的元細胞聚合方法沉颂。）
hdWGCNA 使用 MetacellsByGroups函數(shù)來構造給定單細胞數(shù)據集的metacells表達矩陣。 此函數(shù)為存儲在 hdWGCNA experiment中的metacells數(shù)據集構造一個新的 Seurat 對象悦污。 group.by 參數(shù)確定將在哪些組中構建metacells铸屉。我們只想從來自相同生物樣本的細胞構建metacells，因此通過 group.by 參數(shù)將該信息傳遞給 hdWGCNA 至關重要切端。 此外彻坛，我們通常為每種細胞類型分別構建metacells。 因此踏枣，在這個演示中昌屉，我們按 Sample 和 cell_type 進行分組以達到預期的結果。
k值（The number of cells to be aggregated）應根據輸入數(shù)據集的大小進行調整茵瀑，通常較小的 k 數(shù)可用于小數(shù)據集间驮。 我們一般使用 20 到 75 之間的 k 值。本教程使用的數(shù)據集有 40,039 個細胞瘾婿，每個生物樣本中的細胞數(shù)從 890 到 8,188 不等蜻牢，這里我們使用 k=25烤咧。 可以使用 max_shared 參數(shù)調整元單元之間允許的重疊量偏陪。

注意：metacells聚合方法對extremely underrepresented cell types效果不好抢呆。 例如，在這個數(shù)據集中笛谦，腦血管細胞（周細胞和內皮細胞）是the least represented抱虐，因此將它們排除在此分析之外。 MetacellsByGroups 有一個參數(shù) min_cells 來排除小于指定單元數(shù)的組饥脑。

# 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 數(shù)據運行 Seurat 函數(shù)。 我們可以使用 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 函數(shù)，用于存儲將用于下游網絡分析的給定單元組的轉置表達式矩陣笋颤。這一步默認情況下使用元單元表達式矩陣 (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 函數(shù)來為不同的軟閾值進行parameter sweep喉悴。此功能通過檢查不同power值的結果網絡拓撲結構棱貌，幫助我們指導我們選擇構建共表達網絡的軟閾值。共表達網絡應該具有無標度拓撲粥惧，因此 TestSoftPowers 函數(shù)模擬了共表達網絡在不同軟功率閾值下與無標度圖的相似程度键畴。此外，我們可以使用函數(shù) 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 函數(shù)訪問以進行進一步檢查：

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函數(shù)構建共表達網絡了嘀粱。（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) 是一種常用的指標垄分，用于總結整個共表達模塊的基因表達譜。簡而言之娃磺，模塊特征基因是通過對包含每個模塊的基因表達矩陣的子集執(zhí)行主成分分析 (PCA) 來計算的薄湿。這些 PCA 矩陣中的第一個PC 就是 ME。

hdWGCNA 包含一個函數(shù) ModuleEigengenes 來計算單個單元格中的模塊特征基因偷卧。此外豺瘤，我們允許用戶對 ME 應用 Harmony 批量校正，從而產生協(xié)調的模塊特征基因 (hME)听诸。以下代碼使用 group.by.vars 參數(shù)執(zhí)行由原始樣本協(xié)調的模塊特征基因計算坐求。

# 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 函數(shù)從 Seurat 對象中提取此矩陣静尼，該函數(shù)默認檢索 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 函數(shù)來計算完整單細胞數(shù)據集中的 kME 值豆混，而不是元細胞數(shù)據集篓像。 該函數(shù)本質上計算基因和模塊特征基因之間的成對相關性。 kME可被用于該數(shù)據集中的所有細胞皿伺，但我們建議在用于計算 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 函數(shù)輕松訪問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 函數(shù)或者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

最后編輯于 ：
?著作權歸作者所有,轉載或內容合作請聯(lián)系作者
  • 人面猴
    序言：七十年代末讼油，一起剝皮案震驚了整個濱河市杰赛，隨后出現(xiàn)的幾起案子，更是在濱河造成了極大的恐慌矮台，老刑警劉巖乏屯，帶你破解...
    沈念sama閱讀 217,185評論 6贊 503
  • 死咒
    序言：濱河連續(xù)發(fā)生了三起死亡事件根时，死亡現(xiàn)場離奇詭異，居然都是意外死亡辰晕，警方通過查閱死者的電腦和手機蛤迎，發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 92,652評論 3贊 393
  • 救了他兩次的神仙讓他今天三更去死
    文/潘曉璐 我一進店門，熙熙樓的掌柜王于貴愁眉苦臉地迎上來伞芹，“玉大人忘苛，你說我怎么就攤上這事蝉娜〕希” “怎么了？”我有些...
    開封第一講書人閱讀 163,524評論 0贊 353
  • 道士緝兇錄：失蹤的賣姜人
    文/不壞的土叔 我叫張陵召川，是天一觀的道長南缓。 經常有香客問我，道長荧呐，這世上最難降的妖魔是什么汉形？ 我笑而不...
    開封第一講書人閱讀 58,339評論 1贊 293
  • ?港島之戀（遺憾婚禮）
    正文 為了忘掉前任，我火速辦了婚禮倍阐，結果婚禮上概疆，老公的妹妹穿的比我還像新娘。我一直安慰自己峰搪，他們只是感情好岔冀，可當我...
    茶點故事閱讀 67,387評論 6贊 391
  • 惡毒庶女頂嫁案：這布局不是一般人想出來的
    文/花漫 我一把揭開白布。 她就那樣靜靜地躺著概耻，像睡著了一般使套。 火紅的嫁衣襯著肌膚如雪。 梳的紋絲不亂的頭發(fā)上鞠柄，一...
    開封第一講書人閱讀 51,287評論 1贊 301
  • 城市分裂傳說
    那天侦高，我揣著相機與錄音，去河邊找鬼厌杜。 笑死奉呛，一個胖子當著我的面吹牛，可吹牛的內容都是我干的夯尽。 我是一名探鬼主播瞧壮，決...
    沈念sama閱讀 40,130評論 3贊 418
  • 雙鴛鴦連環(huán)套：你想象不到人心有多黑
    文/蒼蘭香墨 我猛地睜開眼，長吁一口氣：“原來是場噩夢啊……” “哼呐萌！你這毒婦竟也來了馁痴？” 一聲冷哼從身側響起，我...
    開封第一講書人閱讀 38,985評論 0贊 275
  • 萬榮殺人案實錄
    序言：老撾萬榮一對情侶失蹤肺孤，失蹤者是張志新（化名）和其女友劉穎罗晕，沒想到半個月后济欢，有當?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體，經...
    沈念sama閱讀 45,420評論 1贊 313
  • ?護林員之死
    正文 獨居荒郊野嶺守林人離奇死亡小渊，尸身上長有42處帶血的膿包…… 初始之章·張勛 以下內容為張勛視角 年9月15日...
    茶點故事閱讀 37,617評論 3贊 334
  • ?白月光啟示錄
    正文 我和宋清朗相戀三年法褥，在試婚紗的時候發(fā)現(xiàn)自己被綠了。 大學時的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片酬屉。...
    茶點故事閱讀 39,779評論 1贊 348
  • 活死人
    序言：一個原本活蹦亂跳的男人離奇死亡半等，死狀恐怖，靈堂內的尸體忽然破棺而出呐萨，到底是詐尸還是另有隱情杀饵，我是刑警寧澤，帶...
    沈念sama閱讀 35,477評論 5贊 345
  • ?日本核電站爆炸內幕
    正文 年R本政府宣布谬擦，位于F島的核電站切距，受9級特大地震影響，放射性物質發(fā)生泄漏惨远。R本人自食惡果不足惜谜悟，卻給世界環(huán)境...
    茶點故事閱讀 41,088評論 3贊 328
  • 男人毒藥：我在死后第九天來索命
    文/蒙蒙 一、第九天 我趴在偏房一處隱蔽的房頂上張望北秽。 院中可真熱鬧葡幸，春花似錦、人聲如沸贺氓。這莊子的主人今日做“春日...
    開封第一講書人閱讀 31,716評論 0贊 22
  • 一樁弒父案，背后竟有這般陰謀
    文/蒼蘭香墨 我抬頭看了看天上的太陽掠归。三九已至缅叠，卻和暖如春，著一層夾襖步出監(jiān)牢的瞬間虏冻，已是汗流浹背肤粱。 一陣腳步聲響...
    開封第一講書人閱讀 32,857評論 1贊 269
  • 情欲美人皮
    我被黑心中介騙來泰國打工， 沒想到剛下飛機就差點兒被人妖公主榨干…… 1. 我叫王不留厨相，地道東北人领曼。 一個月前我還...
    沈念sama閱讀 47,876評論 2贊 370
  • 代替公主和親
    正文 我出身青樓，卻偏偏與公主長得像蛮穿，于是被迫代替她去往敵國和親庶骄。 傳聞我的和親對象是個殘疾皇子，可洞房花燭夜當晚...
    茶點故事閱讀 44,700評論 2贊 354

推薦閱讀更多精彩內容