應審稿人要求| pseudo bulk差異分析

一、寫在前面

最近有粉絲提問，收到了如下的審稿人意見：

1.png

審稿人認為在單細胞測序過程中岖寞，利用findMarker通過Wilcox獲得的差異基因雖然考慮到了不同組別細胞數(shù)量的不同氧映，但是未能考慮到每組樣本數(shù)量的不同窄俏。因此作者希望納入樣本水平的pseudo-bulk分析能夠有助于確認兩種條件下的差異基因。首先我個人覺得審稿人自己的話有些矛盾城须，使用pseudo-bulk計算差異基因蚤认，豈不是無法考慮到不同組別中樣本數(shù)量的差異？另外這有些吹毛求疵糕伐，在Seurat V5之前砰琢，作者甚至沒有在包里集成pseudo-bulk的函數(shù)與算法(當然也可以自己提取矩陣計算)。難道能說作者發(fā)表的這幾篇Cell和Nature給大家推薦的流程不好嗎：
Hao, Hao, et al., Cell 2021 [Seurat v4]
Stuart, Butler, et al., Cell 2019 [Seurat v3]
Butler, et al., Nat Biotechnol 2018 [Seurat v2] Satija, Farrell, et al., Nat Biotechnol 2015 [Seurat v1]再者說赤炒，scRNA-seq比Bulk RNA-Seq更加的稀疏氯析，將前者模擬為后者參與差異計算，其實也沒那么科學莺褒。當然掩缓，審稿人的觀點也不是全無道理，若能夠通過不同的算法得到相同的差異基因結(jié)果遵岩，的確有較高的說服力你辣。

二巡通、`pseudo-bulk`差異分析走起

2.1 數(shù)據(jù)載入

# 加載R包
library(Seurat)``## 載入需要的程序包：SeuratObject``
## 載入需要的程序包：sp``##
## 載入程序包：'SeuratObject'``
## The following objects are masked from 'package:base':
##
##     intersect, t

# 讀取數(shù)據(jù)：
scRNA <- readRDS('test_data/T1D_scRNA.rds')

# 這個數(shù)據(jù)包含24個樣本：
unique(scRNA$sample)``##  [1] "D_503"  "H_120"  "H_630"  "H_3060" "D_609"  "H_727"  "H_4579" "D_504"
##  [9] "H_3128" "H_7108" "D_502"  "D_497"  "D_506"  "H_409"  "H_6625" "D_610"
## [17] "D_501"  "D_500"  "H_4119" "H_1334" "D_498"  "H_2928" "D_644"  "D_505"``# 包含兩個組別的數(shù)據(jù)：
DimPlot(scRNA,split.by = 'Group')`
![2.png](https://upload-images.jianshu.io/upload_images/28196887-db0a2d0cc9c2d667.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/1240)


## 2.2 差異計算

### (1) pseudo-bulk差異計算

`### 生成擬bulk 數(shù)據(jù) **###**
bulk <- AggregateExpression(scRNA, return.seurat = T, slot = "counts", assays = "RNA",
                           group.by = c("cell_type", "sample", "Group")# 分別填寫細胞類型、樣本變量舍哄、分組變量的slot名稱
   )``## Names of identity class contain underscores ('_'), replacing with dashes ('-')
## Centering and scaling data matrix
##
## This message is displayed once every 8 hours.``# 生成的是一個新的Seurat對象
bulk``## An object of class Seurat
## 41056 features across 345 samples within 1 assay
## Active assay: RNA (41056 features, 0 variable features)
##  3 layers present: counts, data, scale.data`

我們可以像普通`scRNA-seq`的`Seurat`對象一樣宴凉，利用`FindMarkers()`進行差異分析，我們這里用`celltype10`做演示表悬。

`# 取出celltype10對應的對象:
ct10.bulk <- subset(bulk, cell_type == "celltype10")
# 改變默認分類變量：
Idents(ct10.bulk) <- "Group"

# 下面的額計算依賴DESeq2弥锄，做過Bulk RNA-Seq的同學都知道：
**if**(!require(DESeq2))BiocManager::install('DESeq2')``## 載入需要的程序包：DESeq2``## 載入需要的程序包：S4Vectors``## 載入需要的程序包：stats4``## 載入需要的程序包：BiocGenerics``##
## 載入程序包：'BiocGenerics'``## The following object is masked from 'package:SeuratObject':
##
##     intersect``## The following objects are masked from 'package:stats':
##
##     IQR, mad, sd, var, xtabs``## The following objects are masked from 'package:base':
##
##     anyDuplicated, aperm, append, as.data.frame, basename, cbind,
##     colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
##     get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
##     match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
##     Position, rank, rbind, Reduce, rownames, sapply, setdiff, table,
##     tapply, union, unique, unsplit, which.max, which.min``##
## 載入程序包：'S4Vectors'``## The following object is masked from 'package:utils':
##
##     findMatches``## The following objects are masked from 'package:base':
##
##     expand.grid, I, unname``## 載入需要的程序包：IRanges``##
## 載入程序包：'IRanges'``## The following object is masked from 'package:sp':
##
##     %over%``## The following object is masked from 'package:grDevices':
##
##     windows``## 載入需要的程序包：GenomicRanges``## 載入需要的程序包：GenomeInfoDb``## 載入需要的程序包：SummarizedExperiment``## 載入需要的程序包：MatrixGenerics``## 載入需要的程序包：matrixStats``##
## 載入程序包：'MatrixGenerics'``## The following objects are masked from 'package:matrixStats':
##
##     colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
##     colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
##     colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
##     colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
##     colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
##     colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
##     colWeightedMeans, colWeightedMedians, colWeightedSds,
##     colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
##     rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
##     rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
##     rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
##     rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
##     rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
##     rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
##     rowWeightedSds, rowWeightedVars``## 載入需要的程序包：Biobase``## Welcome to Bioconductor
##
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.``##
## 載入程序包：'Biobase'``## The following object is masked from 'package:MatrixGenerics':
##
##     rowMedians``## The following objects are masked from 'package:matrixStats':
##
##     anyMissing, rowMedians``##
## 載入程序包：'SummarizedExperiment'``## The following object is masked from 'package:Seurat':
##
##     Assays``## The following object is masked from 'package:SeuratObject':
##
##     Assays``# 差異計算：
bulk_deg <- FindMarkers(ct10.bulk, ident.1 = "D", ident.2 = "H", # 這樣算出來的Fold Change就是D/H
                         slot = "counts", test.use = "DESeq2",# 這里可以選擇其它算法
     verbose = F# 關閉進度提示
     )``## converting counts to integer mode``## gene-wise dispersion estimates``## mean-dispersion relationship``## final dispersion estimates``head(bulk_deg)# 看一下差異列表``##                        p_val avg_log2FC pct.1 pct.2 p_val_adj
## ENSG00000047346 1.900955e-05 -1.1201537 0.917 1.000  0.780456
## ENSG00000168685 7.606182e-05 -0.5817112 1.000 1.000  1.000000
## ENSG00000131759 9.696591e-05 -2.1668759 0.667 1.000  1.000000
## ENSG00000166750 2.136829e-04 -1.5473545 0.750 0.917  1.000000
## ENSG00000163947 3.127113e-04 -0.9136378 0.917 1.000  1.000000
## ENSG00000239713 4.090397e-04 -2.2337008 0.250 0.917  1.000000`

如何寫循環(huán)計算所有細胞類型的差異基因，就留在這里當習題啦蟆沫。

(2)細胞水平的差異計算

# 整理分組變量：
scRNA$CT_Group <- paste(scRNA$cell_type,scRNA$Group,sep = '_')

# 查看新的分組變量：
unique(scRNA$CT_Group)``##  [1] "celltype12_D" "celltype3_H"  "celltype13_H" "celltype6_H"  "celltype0_D"
##  [6] "celltype12_H" "celltype4_H"  "celltype11_D" "celltype14_H" "celltype9_H"
## [11] "celltype11_H" "celltype2_H"  "celltype0_H"  "celltype7_H"  "celltype14_D"
## [16] "celltype1_D"  "celltype4_D"  "celltype1_H"  "celltype8_H"  "celltype3_D"
## [21] "celltype13_D" "celltype8_D"  "celltype7_D"  "celltype5_H"  "celltype6_D"
## [26] "celltype15_H" "celltype2_D"  "celltype5_D"  "celltype10_H" "celltype9_D"
## [31] "celltype10_D" "celltype15_D"``# 差異計算：
cell_deg <- FindMarkers(scRNA,ident.1 = 'celltype10_D',ident.2 = 'celltype10_H' ,group.by = 'CT_Group')# 同樣得到的是celltype10在D組 vs H組的結(jié)果``## For a (much!) faster implementation of the Wilcoxon Rank Sum Test,
## (default method for FindMarkers) please install the presto package
## --------------------------------------------
## install.packages('devtools')
## devtools::install_github('immunogenomics/presto')
## --------------------------------------------
## After installation of presto, Seurat will automatically use the more
## efficient implementation (no further action necessary).
## This message will be shown once per session

(3)兩種算法的對比

library(dplyr)``##
## 載入程序包：'dplyr'``## The following object is masked from 'package:Biobase':
##
##     combine``## The following object is masked from 'package:matrixStats':
##
##     count``## The following objects are masked from 'package:GenomicRanges':
##
##     intersect, setdiff, union``## The following object is masked from 'package:GenomeInfoDb':
##
##     intersect``## The following objects are masked from 'package:IRanges':
##
##     collapse, desc, intersect, setdiff, slice, union``## The following objects are masked from 'package:S4Vectors':
##
##     first, intersect, rename, setdiff, setequal, union``## The following objects are masked from 'package:BiocGenerics':
##
##     combine, intersect, setdiff, union``## The following objects are masked from 'package:stats':
##
##     filter, lag``## The following objects are masked from 'package:base':
##
##     intersect, setdiff, setequal, union``# 先看一下兩種算法的顯著差異基因數(shù)量  
bulk_sig <- filter(bulk_deg,p_val < 0.05)
nrow(bulk_sig)``## [1] 308``cell_sig <- filter(cell_deg,p_val < 0.05)
nrow(cell_sig)``## [1] 1494

可以看出籽暇，pseudo-bulk得到的差異基因數(shù)量要少很多，畫一個韋恩圖看看二者交集

**if**(!require(VennDiagram))install.packages("VennDiagram")``## 載入需要的程序包：VennDiagram``## 載入需要的程序包：grid``## 載入需要的程序包：futile.logger``venn.plot <- venn.diagram(
 x = list(Bulk = rownames(bulk_sig), Cell = rownames(cell_sig)
),
 category.names = c("Bulk DEG", "Single-Cell DEG"),
 filename = NULL,
 output = TRUE,
 main = "Venn Diagram of Significant Genes"
)
grid.draw(venn.plot)

3.png

可以看出包含關系還是挺明顯的饭庞，那我們再用交集基因的avg_log2FC做一個線性回歸看看兩次差異分析的相關性如何：

# 獲得兩次差異分析共同出現(xiàn)的基因：
inter_gene <- intersect(rownames(bulk_sig),rownames(cell_sig))
# 取出avg_log2FC整理為數(shù)據(jù)框
data4plot <- data.frame(Bulk = bulk_sig[inter_gene,'avg_log2FC'],
  Cell = cell_sig[inter_gene,'avg_log2FC']  )

# 線性回歸分析：
lm.model <- lm(Bulk ~ Cell,data = data4plot)
summary(lm.model)#看一下統(tǒng)計學參數(shù)``##
## Call:
## lm(formula = Bulk ~ Cell, data = data4plot)
##
## Residuals:
##      Min       1Q   Median       3Q      Max
## -1.01613 -0.24964 -0.04723  0.17148  2.19351
##
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.17249    0.02757  -6.257 1.58e-09 ***
## Cell         0.70081    0.01959  35.767  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4107 on 263 degrees of freedom
## Multiple R-squared:  0.8295, Adjusted R-squared:  0.8288
## F-statistic:  1279 on 1 and 263 DF,  p-value: < 2.2e-16``mypara <- coefficients(lm.model)#得到截距和斜率
a <- mypara[2]#斜率
b <- mypara[1]#截距
a <- round(a,2)#取兩位有效數(shù)字
b <- round(b,2)

library(ggplot2)
library(ggpubr)``##
## 載入程序包：'ggpubr'``## The following object is masked from 'package:VennDiagram':
##
##     rotate``# 來個散點圖吧~
lmplot <- ggplot( data4plot,aes(x=Bulk, y=Cell))+
 geom_point(color="black")+
 stat_smooth(method="lm",se=TRUE)+stat_cor(data=data4plot, method = "pearson")+#加上置信區(qū)間戒悠、R值、P值
 ggtitle(label = paste(": y = ", a, " * x + ", b, sep = ""))+geom_rug()+#加上線性回歸方程
 labs(x='Bulk DEG',
      y= 'single-cell DEG')

lmplot``## `geom_smooth()` using formula = 'y ~ x'

4.png

R=0.91舟山，那么R^2就是0.83绸狐，可以看出二者的相關性還是不錯的，就看能不能過審稿人這關啦累盗。

大家有什么新的想法寒矿，歡迎在評論區(qū)留言~

環(huán)境信息

sessionInfo()``## R version 4.4.1 (2024-06-14 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 22631)
##
## Matrix products: default
##
##
## locale:
## [1] LC_COLLATE=Chinese (Simplified)_China.utf8
## [2] LC_CTYPE=Chinese (Simplified)_China.utf8  
## [3] LC_MONETARY=Chinese (Simplified)_China.utf8
## [4] LC_NUMERIC=C                              
## [5] LC_TIME=Chinese (Simplified)_China.utf8    
##
## time zone: Asia/Shanghai
## tzcode source: internal
##
## attached base packages:
## [1] grid      stats4    stats     graphics  grDevices utils     datasets
## [8] methods   base    
##
## other attached packages:
##  [1] ggpubr_0.6.0                ggplot2_3.5.1              
##  [3] VennDiagram_1.7.3           futile.logger_1.4.3        
##  [5] dplyr_1.1.4                 DESeq2_1.44.0              
##  [7] SummarizedExperiment_1.34.0 Biobase_2.64.0            
##  [9] MatrixGenerics_1.16.0       matrixStats_1.4.1          
## [11] GenomicRanges_1.56.1        GenomeInfoDb_1.40.1        
## [13] IRanges_2.38.1              S4Vectors_0.42.1          
## [15] BiocGenerics_0.50.0         Seurat_5.1.0              
## [17] SeuratObject_5.0.2          sp_2.1-4                  
##
## loaded via a namespace (and not attached):
##   [1] RcppAnnoy_0.0.22        splines_4.4.1           later_1.3.2            
##   [4] tibble_3.2.1            polyclip_1.10-7         fastDummies_1.7.4      
##   [7] lifecycle_1.0.4         rstatix_0.7.2           globals_0.16.3        
##  [10] lattice_0.22-6          MASS_7.3-60.2           backports_1.5.0        
##  [13] magrittr_2.0.3          plotly_4.10.4           sass_0.4.9            
##  [16] rmarkdown_2.28          jquerylib_0.1.4         yaml_2.3.10            
##  [19] httpuv_1.6.15           sctransform_0.4.1       spam_2.10-0            
##  [22] spatstat.sparse_3.1-0   reticulate_1.39.0       cowplot_1.1.3          
##  [25] pbapply_1.7-2           RColorBrewer_1.1-3      abind_1.4-5            
##  [28] zlibbioc_1.50.0         Rtsne_0.17              purrr_1.0.2            
##  [31] GenomeInfoDbData_1.2.12 ggrepel_0.9.6           irlba_2.3.5.1          
##  [34] listenv_0.9.1           spatstat.utils_3.1-0    openintro_2.5.0        
##  [37] airports_0.1.0          goftest_1.2-3           RSpectra_0.16-2        
##  [40] spatstat.random_3.3-1   fitdistrplus_1.2-1      parallelly_1.38.0      
##  [43] leiden_0.4.3.1          codetools_0.2-20        DelayedArray_0.30.1    
##  [46] tidyselect_1.2.1        UCSC.utils_1.0.0        farver_2.1.2          
##  [49] spatstat.explore_3.3-2  jsonlite_1.8.8          progressr_0.14.0      
##  [52] ggridges_0.5.6          survival_3.6-4          tools_4.4.1            
##  [55] ica_1.0-3               Rcpp_1.0.13             glue_1.7.0            
##  [58] gridExtra_2.3           SparseArray_1.4.8       mgcv_1.9-1            
##  [61] xfun_0.47               withr_3.0.1             formatR_1.14          
##  [64] fastmap_1.2.0           fansi_1.0.6             digest_0.6.37          
##  [67] R6_2.5.1                mime_0.12               colorspace_2.1-1      
##  [70] scattermore_1.2         tensor_1.5              spatstat.data_3.1-2    
##  [73] utf8_1.2.4              tidyr_1.3.1             generics_0.1.3        
##  [76] data.table_1.16.0       usdata_0.3.1            httr_1.4.7            
##  [79] htmlwidgets_1.6.4       S4Arrays_1.4.1          uwot_0.2.2            
##  [82] pkgconfig_2.0.3         gtable_0.3.5            lmtest_0.9-40          
##  [85] XVector_0.44.0          htmltools_0.5.8.1       carData_3.0-5          
##  [88] dotCall64_1.1-1         scales_1.3.0            png_0.1-8              
##  [91] spatstat.univar_3.0-1   knitr_1.48              lambda.r_1.2.4        
##  [94] rstudioapi_0.16.0       tzdb_0.4.0              reshape2_1.4.4        
##  [97] nlme_3.1-164            cachem_1.1.0            zoo_1.8-12            
## [100] stringr_1.5.1           KernSmooth_2.23-24      parallel_4.4.1        
## [103] miniUI_0.1.1.1          pillar_1.9.0            vctrs_0.6.5            
## [106] RANN_2.6.2              promises_1.3.0          car_3.1-2              
## [109] xtable_1.8-4            cluster_2.1.6           evaluate_0.24.0        
## [112] readr_2.1.5             cli_3.6.3               locfit_1.5-9.10        
## [115] compiler_4.4.1          futile.options_1.0.1    rlang_1.1.4            
## [118] crayon_1.5.3            future.apply_1.11.2     ggsignif_0.6.4        
## [121] labeling_0.4.3          plyr_1.8.9              stringi_1.8.4          
## [124] viridisLite_0.4.2       deldir_2.0-4            BiocParallel_1.38.0    
## [127] munsell_0.5.1           lazyeval_0.2.2          spatstat.geom_3.3-2    
## [130] Matrix_1.7-0            RcppHNSW_0.6.0          hms_1.1.3              
## [133] patchwork_1.2.0         future_1.34.0           shiny_1.9.1            
## [136] highr_0.11              ROCR_1.0-11             broom_1.0.6            
## [139] igraph_2.0.3            bslib_0.8.0             cherryblossom_0.1.0

最后編輯于：2024.10.17 09:33:20

?著作權歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者

人面猴
序言：七十年代末，一起剝皮案震驚了整個濱河市幅骄，隨后出現(xiàn)的幾起案子劫窒，更是在濱河造成了極大的恐慌，老刑警劉巖拆座，帶你破解...
沈念sama閱讀 221,548評論 6贊 515
死咒
序言：濱河連續(xù)發(fā)生了三起死亡事件，死亡現(xiàn)場離奇詭異冠息，居然都是意外死亡挪凑，警方通過查閱死者的電腦和手機，發(fā)現(xiàn)死者居然都...
沈念sama閱讀 94,497評論 3贊 399
救了他兩次的神仙讓他今天三更去死
文/潘曉璐我一進店門逛艰，熙熙樓的掌柜王于貴愁眉苦臉地迎上來躏碳，“玉大人，你說我怎么就攤上這事散怖」矫啵” “怎么了？”我有些...
開封第一講書人閱讀 167,990評論 0贊 360
道士緝兇錄：失蹤的賣姜人
文/不壞的土叔我叫張陵镇眷，是天一觀的道長咬最。經(jīng)常有香客問我，道長欠动，這世上最難降的妖魔是什么永乌？我笑而不...
開封第一講書人閱讀 59,618評論 1贊 296
?港島之戀（遺憾婚禮）
正文為了忘掉前任惑申，我火速辦了婚禮，結(jié)果婚禮上翅雏，老公的妹妹穿的比我還像新娘圈驼。我一直安慰自己，他們只是感情好望几，可當我...
茶點故事閱讀 68,618評論 6贊 397
惡毒庶女頂嫁案：這布局不是一般人想出來的
文/花漫我一把揭開白布绩脆。她就那樣靜靜地躺著，像睡著了一般橄抹。火紅的嫁衣襯著肌膚如雪衙伶。梳的紋絲不亂的頭發(fā)上，一...
開封第一講書人閱讀 52,246評論 1贊 308
城市分裂傳說
那天害碾，我揣著相機與錄音矢劲，去河邊找鬼。笑死慌随，一個胖子當著我的面吹牛芬沉，可吹牛的內(nèi)容都是我干的。我是一名探鬼主播阁猜，決...
沈念sama閱讀 40,819評論 3贊 421
雙鴛鴦連環(huán)套：你想象不到人心有多黑
文/蒼蘭香墨我猛地睜開眼丸逸，長吁一口氣：“原來是場噩夢啊……” “哼！你這毒婦竟也來了剃袍？” 一聲冷哼從身側(cè)響起黄刚，我...
開封第一講書人閱讀 39,725評論 0贊 276
萬榮殺人案實錄
序言：老撾萬榮一對情侶失蹤，失蹤者是張志新（化名）和其女友劉穎民效，沒想到半個月后憔维，有當?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體，經(jīng)...
沈念sama閱讀 46,268評論 1贊 320
?護林員之死
正文獨居荒郊野嶺守林人離奇死亡畏邢，尸身上長有42處帶血的膿包…… 初始之章·張勛以下內(nèi)容為張勛視角年9月15日...
茶點故事閱讀 38,356評論 3贊 340
?白月光啟示錄
正文我和宋清朗相戀三年业扒，在試婚紗的時候發(fā)現(xiàn)自己被綠了。大學時的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片舒萎。...
茶點故事閱讀 40,488評論 1贊 352
活死人
序言：一個原本活蹦亂跳的男人離奇死亡程储，死狀恐怖，靈堂內(nèi)的尸體忽然破棺而出臂寝，到底是詐尸還是另有隱情章鲤，我是刑警寧澤，帶...
沈念sama閱讀 36,181評論 5贊 350
?日本核電站爆炸內(nèi)幕
正文年R本政府宣布咆贬，位于F島的核電站败徊，受9級特大地震影響，放射性物質(zhì)發(fā)生泄漏素征。R本人自食惡果不足惜集嵌，卻給世界環(huán)境...
茶點故事閱讀 41,862評論 3贊 333
男人毒藥：我在死后第九天來索命
文/蒙蒙一萝挤、第九天我趴在偏房一處隱蔽的房頂上張望。院中可真熱鬧根欧，春花似錦怜珍、人聲如沸。這莊子的主人今日做“春日...
開封第一講書人閱讀 32,331評論 0贊 24
一樁弒父案酥泛，背后竟有這般陰謀
文/蒼蘭香墨我抬頭看了看天上的太陽。三九已至嫌拣，卻和暖如春柔袁，著一層夾襖步出監(jiān)牢的瞬間，已是汗流浹背异逐。一陣腳步聲響...
開封第一講書人閱讀 33,445評論 1贊 272
情欲美人皮
我被黑心中介騙來泰國打工捶索，沒想到剛下飛機就差點兒被人妖公主榨干…… 1. 我叫王不留，地道東北人灰瞻。一個月前我還...
沈念sama閱讀 48,897評論 3贊 376
代替公主和親
正文我出身青樓腥例，卻偏偏與公主長得像，于是被迫代替她去往敵國和親酝润。傳聞我的和親對象是個殘疾皇子燎竖，可洞房花燭夜當晚...
茶點故事閱讀 45,500評論 2贊 359