Bioconductor 分析基因芯片数据(转帖）

原文地址：https://www.jianshu.com/p/50bcf4ba9d8a

（代码太多，请看原贴）

还有 系统的使用affy处理芯片，看原贴：

https://www.jianshu.com/p/859a7a7d040f

https://www.jianshu.com/p/ccafb8a184c3

https://www.jianshu.com/p/f44dc5d39bed

https://www.jianshu.com/p/e867ed33f43e

https://www.jianshu.com/p/cc171c5faf46

1. 快速入门

安装并加载所需R包

source("http://bioconductor.org/biocLite.R");

biocLite(“CLL”)

# 载入CLL包，CLL包会自动调用affy包，该包含有一系列处理函数

library(CLL)# read example dataset，（CLL包附带的示例数据集）data("CLLbatch")# pre-process using RMA methodCLLrma<- rma(CLLbatch)# read expression value after pre-processinge <- exprs(CLLrma)# 查看部分数据e[1:5,1:5]

image.png

2. 基因芯片数据预处理

2.1 数据输入

# 载入CLL包，CLL包会自动调用affy包，该包含有一系列处理函数library(CLL)# read example dataset，（CLL包附带的示例数据集）data("CLLbatch")# 查看数据类型data.class(CLLbatch)# 读入所有样本的状态信息data(disease)# 查看所有样本的状态信息disease

image.png

# 查看"AffyBatch"的详细介绍help(AffyBatch)phenoData(CLLbatch)featureData(CLLbatch)protocolData(CLLbatch)annotation(CLLbatch)exprs_matrix <- assayData(CLLbatch)[[1]]exprs_matrix[1:5,1:5]exprs(CLLbatch)[1:5,1:5]

image.png

2.2 质量控制

质量控制分三步，直观观察，平均值方法，数据拟合方法。这三个层次的质量控制分别由image 函数、simpleaffy 包和affyPLM包实现

2.2.1 用image包对芯片数据进行质量评估

# 查看第一张芯片的灰度图像image(CLLbatch[,1])

如果图像特别黑，说明型号强度低；如果图像特别亮，说明信号强度很有可能过饱和

image.png

2.2.2 用simpleaffy包对芯片数据进行质量评估

# 安装simpleaffy包source('http://Bioconductor.org/biocLite.R')biocLite("simpleaffy")library(BiocInstaller)biocLite("simpleaffy")library(simpleaffy)library(CLL)data(CLLbatch)# 获取质量分析报告Data.qc <- qc(CLLbatch)# 图型化显示报告plot(Data.qc)

image.png

第一列是所有样品的名称

第二列检出率和平均背景噪声（往往较高的平均背景噪声都伴随着较低的检出率）

第三列

蓝色实现为尺度因子，取值（-3，3）

"o" 不能超过1.25，否则数据质量不好

“三角型“不能超过3，否则数据质量不好

bioB 说明芯片检测没有达标

2.2.3 用affyPLM包对芯片数据进行质量评估

source('http://Bioconductor.org/biocLite.R')biocLite('affyPLM')library(affyPLM)data(CLLbatch)# 对数据集做回归分析，结果为一个PLMset类型的对象Pset <- fitPLM(CLLbatch)image(CLLbatch[,1])# 根据计算结果，画权重图image(Pset,type="weights",which=1, main="Weights")# 根据计算结果，画残差图image(Pset,type="resids",which=1, main="Residuals")# 根据计算结果，画残差符号图image(Pset,type="sign.resids",which=1, main="Residuals.sign")

image.png

image.png

2.2.4 查看总体质量

一个样品的所有探针组的RLE分布可以用一个统计学中常用的箱型图形表示

source('http://Bioconductor.org/biocLite.R')biocLite("RColorBrewer")library(affyPLM)library(RColorBrewer)library(CLL)# 读入数据data("CLLbatch")# 对数据集做回归计算Psel<-fitPLM(CLLbatch)# 载入一组颜色colors<-brewer.pal(12,"Set3")# 绘制RLE箱线图Mbox(Pset, ylim=c(-1,1),col=colors,main="RLE",las=3);# 绘制NUSE箱线图boxplot(Pset,ylim=c(0.95,1.22),col=colors,main="NUSE",las=3)

从以下两幅图中可以看出CLL1，CLL10 的质量明显有别与其他样品，需要去除

image.png

image.png

2.2.5 查看RNA降解曲线

source('http://Bioconductor.org/biocLite.R')biocLite("affy")library(affy)library(RColorBrewer)library(CLL)data("CLLbatch")# 获取降解数据data.deg <- AffyRNAdeg(CLLbatch)# 载入一组颜色colors <- brewer.pal(12,"Set3")# 绘制RNA降解图plotAffyRNAdeg(data.deg, col=colors)# 在左上部位加注图注legend("topleft", rownames(pData(CLLbatch)), col=colors, lwd=1, inset=0.05, cex=0.5)

image.png

去掉三个质量差的

CLLbatch<-CLLbatch[, -match(c("CLL10.CEL","CLL1.CEL","CLL13.CEL"),                              sampleNames(CLLbatch))]

2.2.6 从聚类分析的角度看数据质量

source('http://Bioconductor.org/biocLite.R')biocLite("gcrma")biocLite("graph")biocLite("affycoretools")library(CLL)library(gcrma)library(graph)library(affycoretools)data(CLLbatch)data("disease")# 使用gcrma算法来预处理数据CLLgcrma<-gcrma(CLLbatch)# 提取基因表达矩阵eset<-exprs(CLLgcrma)# 计算样品两两之间的Pearson相关系数pearson_cor<-cor(eset)# 得到Pearson距离的下三角矩阵dist.lower<-as.dist(1-pearson_cor)# 聚类分析hc<-hclust(dist.lower,"ave")plot(hc)# PCA分析samplenames<-sub(pattern="\\.CEL", replacement ="",colnames(eset))groups<-factor(disease[,2])plotPCA(eset,addtext=samplenames,groups=groups,groupnames=levels(groups))

从聚类分析的结果来看，稳定组（红框）和恶化组分不开

image.png

从主成成份分析来看，也分不开

image.png

2.3 背景校正、标准化和汇总

芯片数据通过质量控制，剔除不合格的样品，留下的样品数据往往要通过三步处理（背景校正、标准化和汇总）才能得到下一步分析所需要的基因表达矩阵

去除背景噪声的过程叫背景校正

标准化目的是使各次/组测量或各种实验条件下的测量可以相互比较

使用一定的统计方法将前面得到的荧光强度值从探针水平汇总到探针组水平

> bgcorrect.methods()[1]"bg.correct""mas""none""rma"> normalize.methods(CLLbatch) [1]"constant""contrasts""invariantset""loess""methods""qspline"[7]"quantiles""quantiles.robust""quantiles.probeset""scaling"> pmcorrect.methods()[1]"mas""methods""pmonly""subtractmm"> express.summary.stat.methods()[1]"avgdiff""liwong""mas""medianpolish""playerout"

参数说明

afbatch输入数据的类型必须是AffyBatch对象

bgcorrect.method背景校正方法

bgcorrect.param指定背景校正方法的参数

normalize.method标准化方法

normalize.param指定标准化方法的参数

pmcorrect.methodPM调整方法

pmcorrect.param指定PM方法参数

summary.method汇总方法

summary.param指定汇总方法的参数

芯片内标准化（Loess）前后MA图

# 使用mas方法做背景校正>CLLmas5<- bg.correct(CLLbatch, method ="mas")# 使用constant方法标准化> data_mas5 <- normalize(CLLmas5, method="constant")#　查看每个样品的缩放系数> head(pm(data_mas5)/pm(CLLmas5),5)CLL11.CELCLL12.CELCLL14.CELCLL15.CELCLL16.CELCLL17.CELCLL18.CELCLL19.CELCLL20.CELCLL21.CELCLL22.CELCLL23.CEL17521811.1558491.0238731.4931931.5493692.0002991.4515761.7765010.98251080.70708280.99587331.43236535668911.1558491.0238731.4931931.5493692.0002991.4515761.7765010.98251080.70708280.99587331.43236522769611.1558491.0238731.4931931.5493692.0002991.4515761.7765010.98251080.70708280.99587331.43236523791911.1558491.0238731.4931931.5493692.0002991.4515761.7765010.98251080.70708280.99587331.43236527517311.1558491.0238731.4931931.5493692.0002991.4515761.7765010.98251080.70708280.99587331.432365CLL24.CELCLL2.CELCLL3.CELCLL4.CELCLL5.CELCLL6.CELCLL7.CELCLL8.CELCLL9.CEL1752181.7060261.1563780.84254190.97750820.98165731.1829631.1149761.136390.89392483566891.7060261.1563780.84254190.97750820.98165731.1829631.1149761.136390.89392482276961.7060261.1563780.84254190.97750820.98165731.1829631.1149761.136390.89392482379191.7060261.1563780.84254190.97750820.98165731.1829631.1149761.136390.89392482751731.7060261.1563780.84254190.97750820.98165731.1829631.1149761.136390.8939248# 查看第二个样品的缩放系数是怎么计算来的> mean(intensity(CLLmas5)[,1])/mean(intensity(CLLmas5)[,2])[1]1.155849

2.4 预处理的一体化算法

预处理方法

方法背景校正方法标准化方法汇总方法

MASSmasconstantmas

dChipmasinvariantsetliwong

RMArmaquantilemedianpolish

MAS5 每个芯片可以单独进行标准化；RMA 采用多芯片模型，需要对所有芯片一起进行标准化。

MAS5 利用MM探针的信息去除背景噪声，基本思路是MP-MM；RMA 不使用MM信息，而是基于PM信号分布采用的随机模型来估计表达值。

RMA处理后的数据是经过2为底的对数转换的，而MAS5不是，这一点非常重要，因为大多数芯片分析软件和函数的输入数据必须经过对数变换。

比较不同算法的处理效果

library(affy)library(gcrma)library(affyPLM)library(RColorBrewer)library(CLL)data("CLLbatch")colors <- brewer.pal(12,"Set3")# use MAS5CLLmas5<- mas5(CLLbatch)# use rma CLLrma<- rma(CLLbatch)# use gcrmaCLLgcrma<- gcrma(CLLbatch)## hist plothist(CLLbatch, main="orignal", col=colors)legend("topright", rownames(pData(CLLbatch)), col=colors,      lwd=1, inset=0.05, cex=0.5, ncol=3)hist(CLLmas5, main="MAS 5.0", xlim=c(-150,2^10), col=colors)hist(CLLrma, main="RMA", col=colors)hist(CLLgcrma, main="gcRMA", col=colors)## boxplotboxplot(CLLbatch, col=colors, las=3, main="original")boxplot(CLLmas5, col=colors, las=3, ylim=c(0,1000), main="MAS 5.0")boxplot(CLLrma, col=colors, las=3, main="RMA")boxplot(CLLgcrma, col=colors, las=3, main="gcRMA")

RMA算法将多条曲线重合到了一起，有利于进一步的差异分析，但却出现了双峰现象，不符合高斯正态分布。很显然gcRMA算法在这里表现的更好。

image.png

image.png

image.png

image.png

三个算法处理后的各样品的中值都十分接近。MAS5算法总体而言还不错，有一定拖尾现象；而gcRMA的拖尾现象比RMA要明显得多。这说明针对低表达量的基因，RMA算法比gcRMA算法表现更好一些。

image.png

image.png

image.png

image.png

通过MA图来查看标准化处理的效

library(gcrma)library(RColorBrewer)library(CLL)library(affy)data("CLLbatch")colors<-brewer.pal(12."Set3")CLLgcrma<-gcrma(CLLbatch)MAplot(CLLbatch[,1:4],pairs=TRUE,plot.method="smoothScatter",cex=0.8,main="original MA plot")MAplot(CLLgcrma[,1:4],pairs=TRUE,plot.method="smoothScatter",cex=0.8,main="gcRMA MA plot")

原始数据中，中值（红线）偏离0，经过gcRMA处理后，中值基本保持在零线上。

image.png

image.png

3 基因芯片数据分析

3.1 选取差异表达基因

######## DEG analysis#　如果没有安装limma包，请取消下面两行注释# library(BiocInstaller)# biocLite("limma")# 导入包library(limma)library(gcrma)library(CLL)# 导入CLL数据data("CLLbatch")data(disease)# 移除 CLL1, CLL10, CLL13CLLbatch<-CLLbatch[, -match(c("CLL10.CEL","CLL1.CEL","CLL13.CEL"),                              sampleNames(CLLbatch))]# 用gcRMA算法进行预处理CLLgcrma<- gcrma(CLLbatch)# remove .CEL in sample namessampleNames(CLLgcrma) <- gsub(".CEL$","", sampleNames(CLLgcrma))# remove record in data disease about CLL1, 10 and 13.CELdisease <- disease[match(sampleNames(CLLgcrma), disease[,"SampleID"]), ]# 获得余下21个样品的基因表达矩阵eset <- exprs(CLLgcrma)# 提取实验条件信息disease <- factor(disease[,"Disease"])# 构建实验设计矩阵design <- model.matrix(~-1+disease)# 构建对比模型，比较两个实验条件下表达数据#　这里的.是简写而不是运算符号contrast.matrix <- makeContrasts(contrasts ="diseaseprogres. - diseasestable",                                levels=design)# 线性模型拟合fit <- lmFit(eset, design)# 根据对比模型进行差值计算 fit1 <- contrasts.fit(fit, contrast.matrix)# 贝叶斯检验fit2 <- eBayes(fit1)# 生成所有基因的检验结果报告dif <- topTable(fit2, coef="diseaseprogres. - diseasestable", n=nrow(fit2), lfc=log2(1.5))# 用P.Value进行筛选，得到全部差异表达基因dif <- dif[dif[,"P.Value"]<0.01,]# 显示一部分报告结果head(dif)

>head(dif)logFCAveExprtP.Valueadj.P.ValB39400_at-0.99978505.634004-5.7273291.482860e-050.10345442.44583541303_at-1.34303064.540225-5.5968131.974284e-050.10345442.154635033791_at1.96479626.8379035.4004993.047498e-050.10345441.713558336131_at0.95742149.9453345.3677413.277762e-050.10345441.639622337636_at2.05340935.4786835.1205195.699454e-050.14391121.078831336122_at0.80086047.1462934.7767211.241402e-040.26121180.2922106

下面逐次介绍这个分析过程的六个关键步骤：构建基因表达矩阵、构建实验设计矩阵、构建对比模型（对比矩阵）、线性模型拟合、贝叶斯检验和生成结果报表。

design

image.png

3.2 注释

# 加载注释工具包library(annotate)# 获得基因芯片注释包名称affydb <- annPkgName(CLLbatch@annotation,type="db")affydb# 如果没有安装，先安装biocLite(affydb)#　加载该包library(affydb, character.only = TRUE)# 根据每个探针组的ID获取对应基因Gene Symbol，并作为新的一列dif$symbols<- getSYMBOL(rownames(dif), affydb)# 根据探针ID获取对应基因Entrez IDdif$EntrezID<- getEG(rownames(dif), affydb)# 显示前几行head(dif)

image.png

3.3 统计分析及可视化

# 加载包library(GOstats)# 提取HG_U95Av2芯片中所有探针组对应的EntrezID，注意保证uniqentrezUniverse <- unique(unlist(mget(rownames(eset), hgu95av2ENTREZID)))# 提取差异表达基因及其对应的EntrezIDentrezSelected <- unique(dif[!is.na(dif$EntrezID),"EntrezID"])# 设置GO富集分析的所有参数params <-new("GOHyperGParams", geneIds=entrezSelected, universeGeneIds=entrezUniverse,              annotation=affydb, ontology="BP", pvalueCutoff=0.001, conditional=FALSE,              testDirection="over")# 对所有的GOterm根据params参数做超几何检验hgOver <- hyperGTest(params)# 生成所有GOterm的检验结果报表bp <- summary(hgOver)# 同时生成所有GOterm的检验结果文件，每个GOterm都有指向官方网站的链接，以获得详细信息htmlReport(hgOver, file="ALL_go.html")# 显示结果前几行head(bp)

image.png

# 安装并加载所需包biocLite("GeneAnswers")library(GeneAnswers)# 选取dif中的三列信息构成新的矩阵，新一列必须是EntrezIDhumanGeneInput <- dif[, c("EntrezID","logFC","P.Value")]# 获得humanGeneInput中基因的表达值humanExpr <- eset[match(rownames(dif), rownames(eset)), ]humanExpr <- cbind(humanGeneInput[,"EntrezID"], humanExpr)# 去除NA数据humanGeneInput <- humanGeneInput[!is.na(humanGeneInput[,1]),]humanExpr <- humanExpr[!is.na(humanExpr[,1]),]# KEGG通路的超几何检验y <- geneAnswersBuilder(humanGeneInput,"org.Hs.eg.db", categoryType ="KEGG",                        testType ="hyperG", pvalueT=0.1, geneExpressionProfile=humanExpr,                        verbose = FALSE)getEnrichmentInfo(y)[1:6]

image.png

biocLite("pheatmap")library(pheatmap)# 从基因表达矩阵中，选取差异表达基因对应的数据selected <- eset[rownames(dif), ]# 将selected矩阵每行的名称由探针组ID转换为对应的基因symbolrownames(selected) <- dif$symbols# 考虑显示问题，这里只画前20个基因的热图pheatmap(selected[1:20,], color = colorRampPalette(c("green","black","red"))(100),        fontsize_row = 4, scale ="row", border_color = NA)

image.png

biocLite("Rgraphviz")library(Rgraphviz)# 显著富集的GO term的DAG关系图ghandle <- goDag(hgOver)# 这个图太大，这里只能选一部分数据构建局部图subGHandle <- subGraph(snodes=as.character(summary(hgOver)[,1]), graph = ghandle)plot(subGHandle)

image.png

source("https://bioconductor.org/biocLite.R")biocLite("KEGG.db")yy<- geneAnswersReadable(y)geneAnswersConceptNet(yy, colorValueColumn ="logFC", centroidSize ="pvalue", output="interactive")

这里我报错了，这个网站解决问题http://planspace.org/2013/01/17/fix-r-tcltk-dependency-problem-on-mac/

image.png

yyy<- geneAnswersSort(yy,sortBy ="pvalue")geneAnswersHeatmap(yyy)

image.png

最后输出会话信息：

> sessionInfo()R version3.4.2(2017-09-28)Platform: x86_64-apple-darwin15.6.0(64-bit)Running under: macOS Sierra10.12.6Matrix products:defaultBLAS:/System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylibLAPACK:/Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dyliblocale:[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8attached base packages: [1] grid      stats4    parallel  stats    graphics  grDevices utils    datasets  methods  base    other attached packages: [1] KEGG.db_3.2.3Rgraphviz_2.22.0pheatmap_1.0.8GeneAnswers_2.20.0RColorBrewer_1.1-2[6] Heatplus_2.24.0MASS_7.3-47RSQLite_2.0RCurl_1.95-4.8bitops_1.0-6[11] igraph_1.1.2GO.db_3.5.0GOstats_2.44.0graph_1.56.0Category_2.44.0[16] Matrix_1.2-12hgu95av2.db_3.2.3org.Hs.eg.db_3.5.0annotate_1.56.1XML_3.98-1.9[21] hgu95av2probe_2.18.0AnnotationDbi_1.40.0IRanges_2.12.0S4Vectors_0.16.0hgu95av2cdf_2.18.0[26] CLL_1.18.0gcrma_2.50.0affy_1.56.0Biobase_2.38.0BiocGenerics_0.24.0[31] limma_3.34.1loaded via a namespace (and not attached): [1] matrixStats_0.52.2bit64_0.9-7GenomeInfoDb_1.14.0tools_3.4.2[5] backports_1.1.1affyio_1.48.0rpart_4.1-11Hmisc_4.0-3[9] DBI_0.7lazyeval_0.2.1colorspace_1.3-2nnet_7.3-12[13] gridExtra_2.3DESeq2_1.18.1bit_1.1-12compiler_3.4.2[17] preprocessCore_1.40.0htmlTable_1.9DelayedArray_0.4.1scales_0.5.0[21] checkmate_1.8.5genefilter_1.60.0RBGL_1.54.0stringr_1.2.0[25] digest_0.6.12foreign_0.8-69DOSE_3.4.0AnnotationForge_1.20.0[29] XVector_0.18.0base64enc_0.1-3pkgconfig_2.0.1htmltools_0.3.6[33] htmlwidgets_0.9rlang_0.1.4BiocInstaller_1.28.0BiocParallel_1.12.0[37] acepack_1.4.1GOSemSim_2.4.0magrittr_1.5GenomeInfoDbData_0.99.1[41] Formula_1.2-2Rcpp_0.12.14munsell_0.4.3stringi_1.1.6[45] yaml_2.1.14SummarizedExperiment_1.8.0zlibbioc_1.24.0plyr_1.8.4[49] qvalue_2.10.0blob_1.1.0DO.db_2.9lattice_0.20-35[53] Biostrings_2.46.0splines_3.4.2locfit_1.5-9.1knitr_1.17[57] tcltk_3.4.2fgsea_1.4.0GenomicRanges_1.30.0geneplotter_1.56.0[61] reshape2_1.4.2fastmatch_1.1-0downloader_0.4latticeExtra_0.6-28[65] data.table_1.10.4-3gtable_0.2.0ggplot2_2.2.1xtable_1.8-2[69] survival_2.41-3tibble_1.3.4rvcheck_0.0.9memoise_1.1.0[73] cluster_2.0.6GSEABase_1.40.1

参考文献

http://www.flypeom.site/bioinformatics/r/2017/10/09/microArray-data-analysis/

喜欢就点个赞，留一些宝贵的意见也是最好的赞赏

作者：thinkando

链接：https://www.jianshu.com/p/50bcf4ba9d8a

來源：简书

简书著作权归作者所有，任何形式的转载都请联系作者获得授权并注明出处。

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

推荐阅读更多精彩内容