mRNA-seq转录组二代测序从raw reads到表达矩阵：上中游分析pipeline

mRNA-seq上中游分析需要用到Linux系统的terminal！

尽管TCGA数据库中已经提供了大量的处理好的表达矩阵数据，但是无论是自身实验需要，还是从NCBI的GEO数据库（储存了许多发表文章样本的测序原始数据）下载文件，都会经常有处理转录组测序下机原始数据的需求。因此我尝试总结出一个具有一定普适性的mRNA-seq上中游pipeline，可以供有需求的朋友进行参考、学习和交流。第一次写shell脚本，望包涵！

转录组应该大家都很熟悉了，就不多介绍了，通常我们可以使用针对RNA的NGS（next generation sequencing）来高通量检测不同条件下细胞组织内的基因（/non coding RNA）的表达情况，此外更高深度的RNA测序还可以用来检测融合基因、可变剪接等。本教程着眼于适用面最广的表达量检测，尝试对GEO数据库中的一批二代RNA-seq数据进行下载和自己重新处理。


从原始文件到表达矩阵！

mRNA-seq的上中游数据处理流程大致可以分为如下步骤：
原始数据质检→低质量碱基去除+低质量reads过滤→（过滤数据的再次质检）→reads比对到参考基因组→格式转换→去除多重比对的低质量比对结果和去重复reads→counts计算（或RPKM/FPKM/TPM计算）

至于这些流程为何被称为上中游分析，是相对以表达矩阵为起点的差异分析、通路富集分析等后续分析流程命名的~

软件准备

由于生信的迅速发展，现在每一步流程都有不少已经开发好的软件可供选择，我这里都是使用了较多人使用的经典软件。

（1）原始数据质检：FastQC

FastQC在reads质检界的地位已经是无可撼动的了。关于FastQC的具体参数和报告解读，字数原因不多赘述，推荐简书教程：利用fastqc检测原始序列的质量

安装

Anaconda依赖python，其中的bioconda子库包含了很多常见的生信软件，并且通过它下载可以直接添加到环境变量（调用时不需要输入绝对路径），非常方便！那么我们需要下载anaconda，添加bioconda查找路径，然后就可以愉快地下载这些软件了。

  1. 下载并配置anaconda3.
#download and install the software.
wget https://repo.anaconda.com/archive/Anaconda3-2020.02-Linux-x86_64.sh -P ~/Downloads
bash ~/Downloads/Anaconda3-2020.02-Linux-x86_64.sh

#Then we have to go through a long agreement, and enter yes and ENTER to complete the installation.
#Then anaconda3 has been installed at /home/xuzihan/anaconda3

#We need to add the command to the environment.
vi ~/.bashrc
#Go to the bottom of this file, press "i" to start text insertion. And add this command:
export PATH="$PATH:/home/xuzihan/anaconda3/bin"
#Press esc to exit the texting mode, then press ":", and enter "wq!" to save the modulation and exit the vim texter.

#reload the environment setting.
source ~/.bashrc

#add some databases to anaconda3 for bioinfo software searching and downloading.
conda config --add channels defaults
conda config --add channels bioconda
conda config --add channels conda-forge
anaconda安装过程的截图
  1. fastqc下载
conda install fastqc
#then press "y" to confirm the installation.

然后就安装好并已经加入环境变量了！是不是很方便！

使用

普通用法很简单，输入fastq文件，指定输出目录即可。我们可以得到每个样本的html网页式的QC报告，里面对这一样本的测序数据各项参数进行了解读。

fastqc raw.reads.fastq -o ./fastqc.output.directory
（2）低质量碱基去除+低质量reads过滤：trimmomatic

Trimmomatic是基于java的使用非常广泛的去除接头和质量过滤的软件，通过它的reads去除了部分reads中也测到的测序建库时连接的外源接头，和可能影响reads序列准确性的低质量碱基，使数据质量更好。关于它的工作原理和模式可参考简述教程：NGS 数据过滤之 Trimmomatic 详细说明

安装

从开发者网站上下载压缩包解压后就可以直接得到jar文件了。

wget http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/Trimmomatic-0.39.zip -P /home/xuzihan/
unzip Trimmomatic-0.39.zip
mv ~/Trimmomatic-0.39 ~/trimmomatic
#Then we will get a folder containing adapters and trimmomatic java package.
使用

本教程使用的测序数据是单端测序（Single End），所以使用SE模式就好。

java -jar /home/xuzihan/trimmomatic/trimmomatic-0.39.jar SE -threads 8 -phred33 input_reads.fastq output_reads.fastq ILLUMINACLIP:/home/xuzihan/trimmomatic/adapters/TruSeq3-SE.fa:4:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:50

参数含义：用java执行jar文件，开启8线程，使用单端模式，测序质量分数体系为phred33（现在的illumina测序都是），输入input_reads.fastq，输出output_reads.fastq。过滤参数包括：使用trimmomatic自带的TruSeq3-SE.fa作为测序接头序列进行剪切；对reads头尾均以碱基质量分数为3作为阈值进行剪切；进行滑窗质检，通过标准为4个碱基构成的滑窗内平均质量分数大于15；剪切后可接受的read最小长度为50.

（3）参考基因组比对：bowtie2

对于有参转录组测序，我们需要将测到的reads回贴到参考基因组上，从而鉴定其所属的基因组位置。这一步的选择其实挺多的，比如适合短序列的bowtie，相对更适合50bp以上长序列的bowtie2，还有BWA，HISAT2，Tophat等。Bowtie2在运行时间上比较有优势，还是很受欢迎的。关于bowtie的详细解释和使用参数，请移步bowtie2 manual：bowtie2 manual

安装

基于Anaconda3就可以直接安装了，一步到位。

conda install bowtie2
#Then press "y" to confirm the installation.
使用

第一次进行参考基因组比对前，需要对各数据库（ENSEMBL/NCBI/UCSC等）发布的标准人类基因组数据进行index。完成index后我们会获得6个以bt2为后缀的index文件，一劳永逸，以后就可以直接使用了。

我个人比较喜欢ENSEMBL数据库的版本，会进行定期更新，此外版本号比较规整，每个版本的参考基因组（fasta文件）和基因注释（gtf文件）会同步发布，不会怕版本间出现错误。

本教程会使用到的参考基因组是刚刚（4.29.20）发布的release100中的Homo_sapiens.GRCh38.dna.primary_assembly.fa，以及对应的release100中的Homo_sapiens.GRCh38.100.gtf注释文件（在最后数counts时才会用到）。（ftp://ftp.ensembl.org/pub/release-100/）

bowtie2-build /The_directory_containig_downloaded_reference_genome_fasta/Homo_sapiens.GRCh38.dna.primary_assembly.fa The_unified_index_prefix
#Then we can get 6 index files: The_unified_index_prefix.1.bt2 ~.2.bt2 ~3.bt2 ~4.bt2 The_unified_index_prefix.rev.1.bt2 ~.rev.2.bt2

然后就是用测序得到的reads进行比对了。由于要下载的是单端测序数据，所以使用-U参数。-p为线程数，-x定义参考基因组index，-S定义输出sam文件名。其他更复杂和相对不常用的参数使用参见manual。

bowtie2 -p 8 -x /directory/The_unified_index_prefix -U input_reads.fastq -S output.sam
（4）格式转化：samtools

比对后输出的格式为sam文件，我们选择将其压缩到二进制的bam文件中。我们还需要对比对上的bam格式的reads按照基因组的位置进行排序，方便进行基因组注释。samtools绝对是这一过程的不二选择。此外，在有其他需求的时候，它也含有其他功能。详细参数解读可右转至：samtools常用命令详解

安装
conda install samtools
#Then press "y" to confirm the installation.
使用

-S参数指输入格式为.sam，-b参数指输出.bam。-o参数指输出文件。

#format transformation.
samtools view -bS input.sam > output.bam
#sort the reads by genomic sites.
samtools sort input.bam -o sorted_output.bam
（5）过滤低质量比对结果：sambamba

相比picard，好像sambamba这个软件小众很多了。不过它的语法挺有意思的，在基于关键词的基础上支持一点自然语言。sambamba的处理速度确实是挺快的！这也是它的一个亮点。具体参数解读见：Sambamba: process your BAM data faster!

安装
conda install sambamba
#Then press "y" to confirm the installation.
使用

-t参数指线程；-h指保留bam文件的header；-f指输入文件类型；-F配合后面的语法表示过滤的条件。这里是表示过滤多重比对/未比对上/鉴定为PCR复制的reads。为了减少假阳性，本教程就将条件设置地严苛一些，仅仅只保留比对到参考基因组上一次的非PCR重复reads。（sam文件中，XS参数是指reads比对到多个参考基因组时，除了最佳比对位置外的第二好得分，它存在也就意味着read存在多次比对）

sambamba view -h -t 8 -f bam -F "[XS] == null and not unmapped and not duplicate" ordered.bam > ordered_filtered.bam
（6）从比对后bam文件得到各个样本的表达矩阵：htseq

TCGA数据库可供下载的开放数据包括htseq-counts，就可以说明这一软件的使用广泛程度了。过滤得到的bam通过htseq处理后得到的counts矩阵就是可以进行差异分析的标准形式了！

安装
conda install htseq
#Then press "y" to confirm the installation.
使用

最常用的就是通过htseq获得各基因组注释下的count数了。-s规定是否为链特异性建库（默认为yes，不懂为啥要用这个默认...）；-f指输入文件形式；输入已排序和过滤的bam文件，以及能和参考基因组对应的基因组注释（这里就是ENSEMBL-release100的Homo_sapiens.GRCh38.100.gtf）；最后输出该样本的counts文件。

htseq-count -f bam -s no ordered_filtered.bam /directory/Homo_sapiens.GRCh38.100.gtf > count.txt
（7）从GEO数据库下载：sra-toolkit

通过sra-toolkit中的prefetch指令，我们只需要知道GEO数据库中想下载样本的SRR序列号，就可以直接下载数据了。GEO数据库的结构和使用不是本教程的重点，所以大家可以自行补充学习！

虽然sra-toolkit可以通过apt-get install命令下载，也可以通过conda install下载，但它在安装后常需要configure（比如重新设置下载文件的所属路径），因此我们需要明确sra-toolkit的安装路径。然而不巧的是好像这俩命令的安装路径不是很好找，所以就还是使用麻烦一点但是更透明的办法：自行下载安装包-设置环境变量。

安装
#download the software.
mkdir ~/sratoolkit && cd ~/sratoolkit
wget http://ftp-trace.ncbi.nlm.nih.gov/sra/sdk/2.10.5/sratoolkit.2.10.5-ubuntu64.tar.gz -P ~/sratoolkit
tar -xzvf sratoolkit.2.10.5-ubuntu64.tar.gz

#modulate the environment.
vi ~/.bashrc
#add the directory containing commands of sratoolkit into the environment.
export PATH="$PATH:/home/xuzihan/sratoolkit/sratoolkit.2.10.5-ubuntu64/bin"
#restart
source ~/.bashrc

#Configure the software: First change the directory to 'bin' of sratoolkit, then use the command to open an interactive window.
cd /home/xuzihan/sratoolkit/sratoolkit.2.10.5-ubuntu64/bin
vdb-config --interactive

弹出的交互式窗口可以使用键盘和鼠标结合操作。

使用

sratoolkit的默认下载路径是~/ncbi/public。下载的原始数据格式为.sra格式。

此外，我使用了美国的服务器，所以使用起来没有网速问题。但是据说国内的朋友直接使用prefetch下载数据的速度是非常慢的，解决办法就是使用aspera协助下载（本处不提及，网上解决办法超多的）。

#download a single SRR file.
prefetch SRRxxxxxx

#download multiple files.
prefetch --option-file the_text_containing_one_series_no_per_row.txt

默认路径下载的sra文件放在了~/ncbi/public/sra中。

测序raw数据通常为.fastq，所以我们继续用sratoolkit中的fastq-dump对下载得到的.sra文件进行格式转换。由于下载的是单端测序数据，所以直接用如下命令就可以在同文件夹中获得相应fastq文件，作为一系列分析的起点。双端测序的格式转换命令是不同的！那再说吧...

fastq-dump /path/SRRxxxx.sra

开始分析

我把以上的命令集成了写成了一个shell脚本，通过输入GEO SRR_ID，GEO SRR_ID_phenotype，phenotype的字符数，和使用的接头序列路径，就可以进行一步完整分析了。

单端测序比双端要简单一些，那就先拿单端练练手吧。教程使用的原文件是2013年在ncbi上发表的卵母细胞（n=3）和合子（n=3）的转录组单端测序数据。

首先我们需要查询得到想从GEO下载的原始数据编号，然后使用nano文本编辑器编写成每行含有一个SRR ID的文件（id.txt），以及对应的分组表型文件：我这里自定义了一下，格式为SRR ID_表型分组信息，每行一个（group.txt）。

注：关于表型分组信息需要特别说明一下：由于脚本编写时文件的内部重命名依赖了分组表型的字符数，所以运行这个脚本需要分组名称的字符数相等。如本例中的分组：oocyte_1/2/3和zygote_1/2/3，都分别是8个字符。所以将8作为第三个参数输入。

mkdir ~/rna.seq.tutorial && cd ~/rna.seq.tutorial
nano id.txt
#SRR445718 SRR445719 SRR445720 SRR445721 SRR445722 SRR445723
nano group.txt
#SRR445718_oocyte_1 SRR445719_oocyte_2 SRR445720_oocyte_3 SRR445721_zygote_1 SRR445722_zygote_2 SRR445723_zygote_3
两个输入文档的展示

接下来编写脚本。

这个脚本不需要超级用户root。它依赖的系统目录构架如下：


系统中的目录构架

所以我们只需要配置好软件；把sratoolkit的下载路径设置为~/ncbi/public或是默认；把rna.seq.tutorial目录建立，并在其中写入脚本，id.txt和group.txt，就可以开始愉快运行脚本了。

脚本内容如下，还是建议通读一下的。

#!/bin/bash
#nano mrna.seq.pipeline.sh
#we can use it by entering: bash mrna_seq_pipeline.sh ~/rna.seq.tutorial/id.txt ~/rna.seq.tutorial/group.txt length_of_the_group_name absolute_directory_of_adapters_for_trimmomatic

#download sra files 
prefetch --option-file $1
cd ../ncbi/public/sra

#get the fastq raw data.
for files in `cat $1`
do 
echo ${files}' is being processed'
fastq-dump ${files}'.sra'
seriesngroup=`grep ${files} $2`
total_length=${#seriesngroup}
group_this_file=${seriesngroup:10:$3} #I postulate the SRR series are 10-digit series numbers.
mv ${files}'.fastq' ${group_this_file}'.fastq'
done
group_names=`ls *.fastq`
echo "fastq raw data of all samples have been all fetched."

#quality control by fastqc.
if [ -d './fastqc.output' ] 
then
    echo 'fastqc.output directory has existed. Good!'
else
    echo 'fastqc.output directory has not been established. Now we will make it!'
    mkdir ./fastqc.output
fi
fastqc *.fastq -o ./fastqc.output
echo "QC on raw reads has been done."

#use trimmomatic to filter the reads with low quality and trim low-quality bases.
if [ -d './trimmed_fastq' ] 
then
    echo 'trimmed_fastq directory has existed. Good!'
else
    echo 'trimmed_fastq directory has not been established. Now we will make it!'
    mkdir ./trimmed_fastq
fi
for files in ${group_names}
do 
echo ${files}' is being filtered and trimmed'
java -jar ~/trimmomatic/trimmomatic-0.39.jar SE -threads 8 -phred33 "./${files}" "./trimmed_fastq/clean_${files}" ILLUMINACLIP:$4:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:50
done
echo "Reads trimming has been done."

#perform quality control again on these clean data.
cd ./trimmed_fastq
if [ -d './clean.fastqc.output' ] 
then
    echo 'clean.fastqc.output directory has existed. Good!'
else
    echo 'clean.fastqc.output directory has not been established. Now we will make it!'
    mkdir ./clean.fastqc.output
fi
fastqc *.fastq -o ./clean.fastqc.output
echo "QC on trimmed reads has been done."

#use bowtie2 to align the reads to the reference.
if [ -d '~/hg.38.ref' ] 
then
    cd ~/hg.38.ref
    if [ -e 'Homo_sapiens.GRCh38.dna.primary_assembly.fa' ]
    then
    echo 'the reference genome has existed. Good!'
    else
    echo 'Now we will download the reference genome.'
        wget ftp://ftp.ensembl.org/pub/release-100/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
        gunzip Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
    fi
else
    echo 'Now we will make the directory and download the reference genome.'
    mkdir ~/hg.38.ref
    wget ftp://ftp.ensembl.org/pub/release-100/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
    gunzip Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
fi
if [ -d '~/hg.38.ref/hg38.index' ] 
then
    cd ~/hg.38.ref/hg38.index
    if [ -e 'hg.38.1.bt2' ] && [ -e 'hg.38.2.bt2' ] && [ -e 'hg.38.3.bt2' ] && [ -e 'hg.38.4.bt2' ] && [ -e 'hg.38.rev.1.bt2' ] && [ -e 'hg.38.rev.2.bt2' ]
    then
        echo 'The reference genome index has existed. Good!'
    else
        echo 'Bowtie2 is working on reference index establishment.'
        bowtie2-build ~/hg.38.ref/Homo_sapiens.GRCh38.dna.primary_assembly.fa hg.38
    fi
else
    echo 'Bowtie2 is working on reference index establishment.'
    mkdir ~/hg.38.ref/hg38.index && cd ~/hg.38.ref/hg38.index
    bowtie2-build ~/hg.38.ref/Homo_sapiens.GRCh38.dna.primary_assembly.fa hg.38
fi
cd ~/ncbi/public/sra/trimmed_fastq
for trimmed_files in `ls *.fastq`
do
    echo ${trimmed_files}' is being aligned with hg38 reference.'
    bowtie2 -p 8 -x ~/hg.38.ref/hg38.index/hg.38 -U ${trimmed_files} -S "${trimmed_files}.sam" --no-unal
done
echo "The bowtie2 alignment on samples has been done."

#use samtools to transform and order the aligned results.
if [ -d '../bam' ]
then
    echo 'bam directory has existed.'
else
    mkdir ../bam
fi
for sams in `ls *.sam`
do
    cd ../bam
    charlength=$((6+$3))
    bam_file=${sams:0:$charlength}
    samtools view -bS "../trimmed_fastq/$sams" > "${bam_file}.bam"
    samtools sort "${bam_file}.bam" -o "${bam_file}_ordered.bam"
    rm "${bam_file}.bam"
    sambamba view -h -t 8 -f bam -F "[XS] == null and not unmapped and not duplicate" "${bam_file}_ordered.bam" > "${bam_file}_ordered_filtered.bam"
    rm "${bam_file}_ordered.bam"
done
echo "alignment reordering and filtering has been done."

#get the annotation for the genome and transform .gff file to .gtf file for quantification.
if [ -e '~/hg.38.ref/Homo_sapiens.GRCh38.100.gtf' ]
then
    echo "GTF file has all set!"
else
    echo "Now download the GTF annotation file."
    wget ftp://ftp.ensembl.org/pub/release-100/gtf/homo_sapiens/Homo_sapiens.GRCh38.100.gtf.gz -P ~/hg.38.ref
    gunzip Homo_sapiens.GRCh38.100.gtf.gz
fi

#Use HTSeq to assemble the reads and complete the annotation.
if [ -d '../count_matrix' ]
then
    echo 'count_matrix directory has existed.'
else
    mkdir ../count_matrix
fi
cd ~/ncbi/public/sra/bam
for bams in `ls * | grep "filtered.bam"`
do
    cd ../count_matrix
    count_file_name=${bams:6:$3}
    htseq-count -f bam -s no "../bam/$bams" ~/hg.38.ref/Homo_sapiens.GRCh38.100.gtf > "${count_file_name}_count.txt"
done
echo "gene annotation and counts extraction have been done."

最终开始运行：共4个参数，第一个为id.txt，第二个group.txt，第三个为表型分组的字符数，第四个为测序用到的接头序列fasta。

bash ~/rna.seq.tutorial/mrna.seq.pipeline.sh ~/rna.seq.tutorial/id.txt ~/rna.seq.tutorial/group.txt 8 ~/trimmomatic/adapters/TruSeq3-SE.fa

运行过程中terminal的一些显式进程大致如下，我截了一些图：


脚本运行过程中的一些进程显示

80%的修剪后保留率，80%的比对率都提示数据的总体质量还不错。

我们看看数据的两次fastqc结果：


两次fastqc结果的比较（左-原始数据、右-trimmomatic修剪后数据）

左边是原始数据的QC report，右边是trimmomatic处理后的QC report。

（1）碱基质量分数：可见原始数据中，随着reads的延伸，碱基的测序质量快速降低，总体已经达到不能接受的程度。在过滤修剪后，每个碱基位置的总体质量分数分布均在绿色范围，说明过滤很好地提升了测序reads的质量。

（2）reads长度：可见原始数据中，测到的reads都是100bp长度。在修剪后，reads长度出现差异，分布在50-100bp间，说明部分reads中有碱基修剪；这也与我们设置可接受最小reads长度为50bp相符合。

（3）接头内容：可见原始数据中，reads的3’端检测到少量的接头内容，可能是由于mRNA片段本身过短所导致的接头被测序。通过修剪后，右边的report显示已经基本检测不到接头序列的存在，说明修剪成功。

最后我们可以得到如下文件，我在terminal中查看了相关文件夹中的文件：
Screen Shot 2020-05-11 at 11.10.30 AM.png

count_matrix中的6个count.txt就是每个样本的表达矩阵了！如下图所示，展示了其中一个样本的表达矩阵：


一个样本的表达矩阵概览

通过htseq获得的每个样本的表达矩阵有2列，第一列为基因的ENSEMBL ID，第二列为每个注释的raw counts。此外，最后五行描述了htseq计算的结果，从__not_aligned和__alignment_not_unique均为0即可知bowtie2和sambamba的过滤是有效的（在脚本中已经设置Bowtie2在比对时不输出未必对的reads，以及sambamba过滤多次比对的reads）。

最后我们将所有样本的表达矩阵合并成一个总矩阵后，就可以进行下游一系列分析了。

有兴趣的朋友可以学习，欢迎大家相互交流！觉得有帮助就点个赞吧！

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