群体进化分析--基因流及群体历史分析——我愿称之为史上最全

基因流（也称基因迁移）
是指从一个物种的一个种群向另一个种群引入新的遗传物质，从而改变群体“基因库”的组成。通过基因交流向群体中引入新的等位基因，是遗传变异一个非常重要的来源，影响群体遗传多样性，产生新的性状组合

一、Treemix分析

Treemix软件使用全基因组的等位基因频率数据，推断多个群体的分化和混合的模式。该软件输入数据为多个群体的等位基因频率数据，可以生成这些群体的最大似然树，并且可以推测群体混合事件。

1.计算等位基因频率
# 1.1 vcf文件转换为tped格式，生成sample.tped和sample.tfam文件。
vcftools --vcf QC.sample-select-geno005-maf003.vcf --plink-tped --out sample
# 1.2 文件修改
sample.tfam文件修改第一列数据为breed ID（品种）。
# 1.3 提取文件sample.pop.cov，格式为：共三列，前两列与修改后的sample.tfam前两列一样，为群体ID和样本ID，第三列和第一列一致，tab分隔。
cat sample.tfam |awk '{print $1"\t"$2"\t"$1}' >sample.pop.cov
# 1.4 计算等位基因组的频率，生成plink.frq.strat和plink.nosex文件
plink --threads 12 --tfile sample --freq --allow-extra-chr --chr-set 29 --within sample.pop.cov 
# 1.5 压缩等位基因频率文件
gzip plink.frq.strat 
# 1.6 转换格式
#用treemix自带脚本进行格式转换，notes：输入输出都为压缩文件，plink2treemix.py使用python2并需要绝对路径（否则报错）。
python2 /home/sll/miniconda3/bin/plink2treemix.py plink.frq.strat.gz sample.treemix.in.gz
2.treemix推断基因流
touch treemix.sh
##======= 多次分析以评估最佳m值=======
## 比如m取1-10(常用1-5,1-10)，每个m值重复5次(至少两次)
for m in {1..5}
do
    {
    for i in {1..5}
    do
        {
           treemix -i sample.treemix.in.gz -o sample.${m}.${i} -root wBGU -m ${m} -k 500 -noss
        }
    done
    }
done
##===============================
nohup bash treemix.sh &
# -i 指定基因频率输入文件
# -o 指定输出文件前缀
# -root指定外类群(指定的是居群名称)，多个用逗号分隔；最好指定，否则后面plot_tree画树没找到更换外类群的参数会很麻烦
# -m 为the number of migration edges即基因渗入的次数
# -k 500 因为SNP之间不是独立位点，为了避免连锁不平衡，用k参数指定SNP数量有连锁，比如这里指定用1000个SNP组成的blocks评估协方差矩阵
# -noss 关闭样本量校正。TreeMix计算协方差会考虑每个种群的样本量，有些情况(如果有种群的样本只有1个)会过度校正，可以关闭。
3.最佳迁移边缘个数选择
# 将生成的llik、cov、modelcov文件放置在同一文件夹，使用R包OptM分析：

library(OptM)
linear = optM("./data") ##文件夹名
plot_optM(linear)

# 生成图中，当Δm值最小时的migration edges为最佳迁移边缘个数
4.基因渗入作图
## 使用m=最佳迁移边缘个数结果文件做图
source("plotting_funcs.R") #treemix scr文件夹中R脚本
plot_tree("TreeMix") #TreeMix为结果文件前缀

##=====绘制混合树====
prefix="sample.3"  ## treemix产生的结果文件前缀
library(RColorBrewer)
library(R.utils)
par(mfrow=c(2,3))
for(edge in 1:3){
  plot_tree(cex=0.8,paste0(prefix,".",edge)) ## ""中的东西取决于你结果文件前缀和后面的东西的连接符号
  title(paste(edge,"repetition"))
} # 放的是m012345，则0:5
treemix系统发育分析（ML）

群体间的系统发育树，每个分支代表一个群体。

treemix -i sample.treemix.in.gz -o sample.ML.tree -root wBGU -k 1000

结果展示：

ML-tree.png

二、f3、f4检验以及D检验

D 统计量（Patterson’s D statistic）是目前最广为使用的渗入检测方法之一。D 统计量需要四个群体，这里称之为 P1、P2、P3 和 O。他们的系统发育关系如图 1-15 所示，为(((P1,P2)P3,)O)。其中，P1 和 P2 为姐妹群，O 为外群。对于这四个群体中存在的双等位 SNP，以外群O 中的类型作为祖先型 A（即 ancestral allele），另一种则为衍生型 B（即 derived allele）。在已知的系统发育关系的情况下，基因组上绝大多数 SNP 在这四个群体中的排布模式应该为 BBAA，即 P1 和 P2 同为衍生型等位基因，而 P3 和 O 同为祖先型的等位基因。BBAA 的模式是符合系统发育关系的，但如果 P2 和 P3 之间发生了基因流，那么则会产生大量 ABBA 模式的位点，即 P2 和 P3 共享了同一种等位基因。反之，如果 P1 和 P3 之间发生了渗入，则会产生大量 BABA 模式的位点。
image.png

但实际上，由于 P1、P2 和 P3 的共同祖先在某些位点上可能是存在多态性的，这种多态性可能会随机分配给这三个群体，从而导致 ABBA 或者 BABA 的情况。这种现象通常被称为不完全谱系分选（incomplete lineage sorting，ILS）。因此单独去检测ABBA 或者 BABA 模式的位点数量无法判断它们是由于 ILS 还是渗入导致的。但由于ILS 是不受选择的，所以它产生的 ABBA 和 BABA 模式的位点数量应该是大致相当的。所以我们可以通过比较 ABBA 和 BABA 的数量是否有显著的差异来判断是单纯的 ILS还是发生了渗入。

image.png
其中，CABBA表示 ABBA 模式位点的数量。D 统计为符合 ABBA 模式的位点数量与 BABA 模式的位点数量之差，除以 ABBA 模式和 BABA 模式位点数量之和。故 D统计量为正时，ABBA 的数量更多，因此可能是 P2 和 P3 之间发生了渗入。反之，D统计量为负时，BABA 数量更多，可能是 P1 和 P3 之间发生了渗入。由于 D 统计量的计算是基于对 ABBA 和 BABA 模式位点的计数，因此它也被称之为 ABBA-BABA test。
需要注意的是，D statistic 中四个群体的排列顺序以及正负值代表的含义在不同的软件中可能并不相同，如 ADMIXTOOLS中D为正值时可能是 P1 和 P3 之间发生了渗入、ANGSD doAbbababa2和 Dsuite.

●不完全谱系分选(ILS,Incomplete Lineage sorting):位点树（基因树）和物种树不一致的现象。例如ABC物种分化前某一位点存在多态性，为0/1。随着C分化出去，而此多态性的位点在C中发生了固定，为1，而在AB祖先中是以多态存在。当AB发生分化时，此等位以随机的方式分别进入AB物种中。当B的位点和C相同时，即发生BC关系更近的现象。我们以为是BC之间存在基因流的现象，而实际上为不完全谱系分选。

1、f3统计量
image.png
                     f3(A,B;C)=(c-a)(c-b)；A, B为参考群体，C为测试群体
1）Treemix中的threepop

准备输入文件

先提取所用品种的样本
# 1.1 vcf文件转换为tped格式，生成sample.tped和sample.tfam文件。
vcftools --vcf QC.sample-select-geno005-maf003.vcf --plink-tped --out sample
# 1.2 文件修改
sample.tfam文件修改第一列数据为breed ID（品种）。
# 1.3 提取文件sample.pop.cov，格式为：共三列，前两列与修改后的sample.tfam前两列一样，为群体ID和样本ID，第三列和第一列一致，tab分隔。
cat sample.tfam |awk '{print $1"\t"$2"\t"$1}' >sample.pop.cov
# 1.4 计算等位基因组的频率，生成plink.frq.strat和plink.nosex文件
plink --threads 12 --tfile sample --freq --allow-extra-chr --chr-set 29 --within sample.pop.cov 
# 1.5 压缩等位基因频率文件
gzip plink.frq.strat 
# 1.6 转换格式
#用treemix自带脚本进行格式转换，notes：输入输出都为压缩文件，plink2treemix.py使用python2并需要绝对路径（否则报错）。
python2 /home/sll/miniconda3/bin/plink2treemix.py plink.frq.strat.gz sample.treemix.in.gz

threepop -i sample.treemix.in.gz -k 500 > 3pop
cat 3pop |grep -v Estimating |grep -v nsnp|tr ';' ' ' > 3pop.txt

如果F3算出来的值是一个负值的话（同时满足统计学检验的标准，Z-scores=F3/SE(F3)要小于-3），认为F3具有显著的统计学效应，即C中的一些基因来自于A和B的混合。值得注意的是，上面反过来不成立，F3大于0，并不认为C中不含有A和B的混合

如果令C物种为外来物种的话，那么C其实无论如何都不可能是AB的混合，而F3值越大，可以用来表征A,B之间的相似性越强(outgroup f3)
2） Admixtools中的pq3Pop(可做outgroup f3判断A和B种群的关系)：
  1. AdmixTools需要特征（eigenstrat）文件，要将vcf文件转换为eigenstrat。
bash convertVCFtoEigenstrat.sh QC.sample-select-geno005-maf003.vcf
convertVCFtoEigenstrat.sh
#!/bin/bash
# Script to convert vcf to eigenstrat format for ADMIXTOOLS
# Written by Joana Meier
# It takes a single argument: the vcf file (can be gzipped) and 
# optionally you can specify --renameScaff if you have scaffold names (not chr1, chr2...)
# Here, you can change the recombination rate which is currently set to 2 cM/Mb 
rec=2
# It requires vcftools and admixtools
# for some clusters, it is needed to load these modules:
# module load gcc/4.8.2 vcftools openblas/0.2.13_seq perl/5.18.4 admixtools
renameScaff="FALSE"
# If help is requested
if [[ $1 == "-h" ]]
then
 echo "Please provide the vcf file to parse, and optionally add --renameScaff if you have scaffolds instead of chromosomes"
 echo "Usage: convertVCFtoEigenstrat.sh <vcf file> --renameScaff (note, the second argument is optional)"
 exit 1
# If the second argument renameScaff is given, set it to True
elif [[ $2 == "--renameScaff" ]]
then
 renameScaff="TRUE"
# If no argument is given or the second one is not -removeChr, give information and quit
elif [ $# -ne 1 ]
then
 echo "Please provide the vcf file to parse, and optionally add --renameScaff if you have scaffolds instead of chromosomes"
 echo "Usage: ./convertVCFtoEigenstrat.sh <vcf file> --renameScaff (note, the second argument is optional)"
 exit 1
fi
# Set the first argument to the file name
file=$1
file=${file%.gz}
file=${file%.vcf}
# if the vcf file is gzipped:
if [ -s $file.vcf.gz ]
then
 # If renaming of scaffolds is requested, set all chromosome/scaffold names to 1
 if [ $renameScaff == "TRUE" ]
 then
  echo "setting scaffold names to 1 and positions to additive numbers"
  zcat $file".vcf.gz" | awk 'BEGIN {OFS = "\t";add=0;lastPos=0;scaff=""}{
    if($1!~/^#/){
       if($1!=scaff){add=lastPos;scaff=$1}
       $1="1"
       $2=$2+add
       lastPos=$2
     }
     print $0}' | gzip > $file.renamedScaff.vcf.gz
  # Get a .map and .ped file (remove multiallelic SNPs, monomorphic sites and indels)
  vcftools --gzvcf $file".renamedScaff.vcf.gz" \
         --plink --mac 1.0 --remove-indels --max-alleles 2 --out $file
 else
 # Get a .map and .ped file (remove multiallelic SNPs, monomorphic sites and indels)
 vcftools --gzvcf $file".vcf.gz" \
         --plink --mac 1.0 --remove-indels --max-alleles 2 --out $file
 fi
# if the file is not gzipped
else
 # If renaming of scaffolds is requested, set all chromosome/scaffold names to 1
 if [ $renameScaff == "TRUE" ]
 then
  echo "setting scaffold names to 1 and positions to additive numbers"
  awk 'BEGIN {OFS = "\t";add=0;lastPos=0;scaff=""}{
    if($1!~/^#/){
       if($1!=scaff){add=lastPos;scaff=$1}
       $1="1"
       $2=$2+add
       lastPos=$2
     }
     print $0}' $file.vcf | gzip > $file.renamedScaff.vcf.gz
  # Get a .map and .ped file (remove multiallelic SNPs, monomorphic sites and indels)
  vcftools --gzvcf $file".renamedScaff.vcf.gz" \
         --plink --mac 1.0 --remove-indels --max-alleles 2 --out $file
 else
 vcftools --vcf $file".vcf" --plink --mac 1.0 --remove-indels --max-alleles 2 --out $file
 fi
fi
# Change the .map file to match the requirements of ADMIXTOOLS by adding fake Morgan positions (assuming a recombination rate of 2 cM/Mbp)
awk -F"\t" -v rec=$rec 'BEGIN{scaff="";add=0}{
        split($2,newScaff,":")
        if(!match(newScaff[1],scaff)){
                scaff=newScaff[1]
                add=lastPos
        }
        pos=add+$4
    count+=0.00000001*rec*(pos-lastPos)
        print newScaff[1]"\t"$2"\t"count"\t"pos
        lastPos=pos
}' ${file}.map  | sed 's/^chr//' > better.map
mv better.map ${file}.map
# Change the .ped file to match the ADMIXTOOLS requirements
awk 'BEGIN{ind=1}{printf ind"\t"$2"\t0\t0\t0\t1\t"; 
 for(i=7;i<=NF;++i) printf $i"\t";ind++;printf "\n"}' ${file}.ped > tmp.ped
mv tmp.ped ${file}.ped
# create an inputfile for convertf
echo "genotypename:    ${file}.ped" > par.PED.EIGENSTRAT.${file}
echo "snpname:         ${file}.map" >> par.PED.EIGENSTRAT.${file}
echo "indivname:       ${file}.ped" >> par.PED.EIGENSTRAT.${file}
echo "outputformat:    EIGENSTRAT" >> par.PED.EIGENSTRAT.${file}
echo "genotypeoutname: ${file}.eigenstratgeno" >> par.PED.EIGENSTRAT.${file}
echo "snpoutname:      ${file}.snp" >> par.PED.EIGENSTRAT.${file}
echo "indivoutname:    ${file}.ind" >> par.PED.EIGENSTRAT.${file}
echo "familynames:     NO" >> par.PED.EIGENSTRAT.${file}
# Use CONVERTF to parse PED to eigenstrat
convertf -p par.PED.EIGENSTRAT.${file}
# change the snp file for ADMIXTOOLS:
awk 'BEGIN{i=0}{i=i+1; print $1"\t"$2"\t"$3"\t"i"\t"$5"\t"$6}' $file.snp > $file.snp.tmp
mv $file.snp.tmp $file.snp
  1. 修改.ind文件的第三列为品种id
  2. 提供A、B、C群体的pop文件，3列(其中C为目标种群)（做outgroup f3时，C为外群，判断A和B的关系，脚本一致）
  3. 修改脚本文件par.PED.EIGENSTRAT.QC.sample-select-out-geno005-maf003，里的东西为：
genotypename:   QC.sample-select-out-geno005-maf003.eigenstratgeno
snpname:        QC.sample-select-out-geno005-maf003.snp
indivname:      QC.sample-select-out-geno005-maf003.ind
popfilename:    pop.txt
inbreed: YES   ##(做outgroup f3应删除这一行，目标群体存在近交，则加上这行)
  1. qp3Pop分析
qp3Pop -p par.PED.EIGENSTRAT.QC.sample-select-out-geno005-maf003 > 3pop_qp3pop

结果展示：

image.png
2、f4统计量
image.png
                    f4(A,B;C,D)=(a-b)(c-d); 判断A,B与C,D之间是否存在基因流
                    其中A、B和C、D直接互为对立姊妹群
Treemix中的fourpop
fourpop -i sample.treemix.in.gz -k 500 > 4pop
cat 4pop |grep -v Estimatin|grep -v nsnp |tr ',' ' '|tr ';' ' ' > 4pop.txt

若B,C之间发生基因流，B,C等位基因频率趋同，介于A,D之间，f4<0，|Z| > 3
若A,D之间发生基因流，A,D等位基因频率趋同，介于B,C之间，f4<0, |Z| > 3
A,C或B,D之间发生基因流，f4>0, |Z| > 3,有统计学意义

值得注意的是，将A当作外来物种，那么我们就可以直接通过F4的正负来判断基因流了。若F4>0，则B与D之间有基因流；若F4<0，则B与C之间有基因流

这个是有关f3、f4的详细价绍，可以将看看，讲的很好【概念】关于基因流中的F3和F4统计量 - 知乎 (zhihu.com)

3、D统计量

也叫ABBA-BABA检验，用于检测是否存在由“Admixture”造成的显著不符合拓扑树结构事件的方法。使用Z检验作为显著性检验算法（X,Y; W,Outgroup）

image.png
4 populations (W, X, Y, Z)
num = nBABA-nABBA= (w − x)(y − z )
den = nBABA+nABBA= (w + x − 2wx)(y + z − 2yz )
D = num/ den （这个是官网上给的公式）
也就是                       D=(nBABA-nABBA)/(nBABA+nABBA)

有的文献里面指出D=(nABBA-nBABA)/(nABBA+nABAB)，可能是因为不同软件的缘故，Dsuite软件中是这个和公式，但是不管他的公式如何，只要记住，Z为外群的情况下，D > 0，W与Y之间存在基因流； D < 0，X与Y之间存在基因流（官网上的说明，原话如下，Z-score = D/std. err (标准差)和D值是同正负的）

If the Z-score is +ve, then the gene flow occured either between W and Y or X and Z 
If the Z-score is -ve, then the gene flow occured either between W and Z or X and Y.
AdmixTools中的qpDstat(这里的123步的这些文件和前面的f3 qp3Pop用的一致)

1、修改.ind文件的第三列为品种id
2、提供A、B、C、D群体的pop文件，4列(注：D为outgroup，判断ABC的关系)
3、修改脚本文件par.PED.EIGENSTRAT.QC.sample-select-out-geno005-maf003，里的东西为：

genotypename:   QC.sample-select-out-geno005-maf003.eigenstratgeno
snpname:        QC.sample-select-out-geno005-maf003.snp
indivname:      QC.sample-select-out-geno005-maf003.ind
popfilename:    pop.txt
f4mode: YES ##此选项为进行f4检验

4、qpDstat分析

qpDstat -p par.PED.EIGENSTRAT.QC.sample-select-out-geno005-maf003 > 4pop.txt
结果展示：
image.png
这个是整理后的结果，输出结果没有列名，我这个列名是按照他官网上说明来的，不过这里他好像用的公式是
D=(nBABA-nABBA)/(nABBA+nBABA)，道理是一样的，看上面那个图就能理解了

D统计量=0，说明总体ABBA和BABA的数量相同，不存在明显的渗入；如果不等于0，则或许存在渗入。
Z>3和<-3分别是正负显著
也就是用于判断X,Y与W之间是否发生基因流

三、渗入片段检测（单倍型）

基因组上的哪些片段是通过基因交流获得的以及它们的生物学功能

●适应性基因渗透包含两个过程：基因渗透和适应。
●基因渗透是指遗传成分从一个群体通过有性生殖的方式进入到另一个群体。产生杂交的前提是两个群体之间没有生殖隔离或只有不完全生殖隔离。只有这样才能通过与可育Ｆ１代回交达到遗传成份垂直传递的目的。
●通过基因渗透得到的区域，如果在该群体内没有受到选择，那么会因为随机遗传漂变的影响，逐渐在基因组中消失。只有当基因渗透得到的区域受到自然或人工选择，它才能在基因组中保留一定的频率。

软件对比：

与HAPMIX, LAM-LD, RFmix需要设置多个参数，HAPMIX需要提供遗传图谱、充足率、突变率和平均代时等较难获得的生物学参数。Loter除了单倍型数据以外不需要其他任何生物学参数即可完成祖先血统推断。与HAPMIX, LAM-LD, RFmix软件相比，在考虑人工模拟和人工混群的情况下，随着基因渗入次数的增加，Loter软件收到的影响较小。

1、Loter软件：

提供足够的祖先物种的基因组数据来作为来源群体（source populations），同时也需要发生渗入样本的基因组数据。对于每一个渗入个体的每一个变异位点，loter均会进行祖先来源判断，从所提供的来源群体中选择一个可能性最大的群体作为其来源。

需要提供一个几乎没有发生渗入或者渗入比例极低的群体作为对照。

loter_cli -r ref_phase.vcf -a admix_phase.vcf -f vcf -o interval.txt

 -r REF [REF ...], --ref REF [REF ...]  files storing input reference haplotypes.可以是多个，空格分开输入
 -a ADM, --adm ADM  file storing input admixed haplotypes.
 -o OUTPUT, --output OUTPUT
                        file to store the infered ancestries of admixed haplotypes, saved as a numpy array by default, or as a text CSV file if the file extension is '.txt'.
 -f 输入文件的格式；有npy(默认), txt（csv）, vcf可选，所有的输入文件格式一致
 -n CPU数量
 -pc Run the phase correction module after the inference algorithm (only available for data from diploid organisms)

结果里面均为0101。本人推断0为祖先型，1为衍生型，可用Excel的条件格式去看他的区域。

2、RFMix软件：

1）、过滤（去除多等位基因位点）
2）、beagle填充

java -jar -Xmn12G -Xms24G -Xmx48G  /home/software/beagle.25Nov19.28d.jar gt=tibetan-36.filter-nchr.recode.vcf  out=tibetan-36.filter-nchr.recode.vcf.beagle ne=36

3)、提取样本并分染色体（多个参考群体在一个vcf文件中）

bcftools view -S id.txt  common_89_cattle_851_ASIA.vcf  -Ov > sample-select.vcf

for i in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28;
do /home/software/vcftools/src/cpp/vcftools --gzvcf tibetan-36.filter-nchr.recode.vcf.beagle.vcf.gz --chr ${i} --recode --recode-INFO-all --out ${i}.test;
done

4)、RFMix推断渗入区域

rfmix -f BC.1.test.recode.vcf -r ANG.1.test.recode.vcf MSU.1.test.recode.vcf -m map.txt -g genetic.map -o outer --chromosome=1

-f 目标群体
-r 参考群体
-m map文件包含两列；1：参考群体样本ID，2：参考群体品种ID
-g genetic.map 文件包含三列；1：染色体号，2：物理距离(bp)，3：遗传距离（cM）, 1cM=1Mb

结果会产生四个文件。

四、LD衰减-- PopLDdecay

LD 的衰减是受重组率和重组代数影响；研究 LD 的衰减可以揭示群体重组的历史，与等位基因的频率相结合，LD 衰减也可以用于检测正向选择
PopLDdecay（软件由西安交通大学杨铁林老师团队开发 （PopLDdecay: a fast and effective tool for linkage disequilibrium decay analysis based on variant call format files），他们实验室的宗旨是，如果觉得某个软件不好用，那就自己写一个，确实厉害

1、计算LDdeacy
整个vcf文件中群体的LDdeacy
PopLDdecay -InVCF QC.sample-select-out-geno005-maf003.vcf -MaxDist 1000 -OutType 3 -OutStat overlap.all.stat
计算vcf文件中的某个群体的LDdeacy(一般用这个)
示例：
PopLDdecay -InVCF QC.sample-select-out-geno005-maf003.vcf -OutStat  GroupA_sample.stat -MaxDist 1000 -SubPop GroupA_sample.list

实操：
制作一列的ID的每个群体txt文件

(保证当前目录没有其他txt文件的情况下，运行以下脚本）
ls *.txt | cut -d '.' -f 1 | sort -u | while read id; do (nohup PopLDdecay -InVCF QC.sample-select-out-geno005-maf003.vcf -OutStat ${id}_sample.stat -MaxDist 1000 -SubPop ${id}.txt &); done

# -InVCF：输入vcf文件。
# -MaxDist：最长Decay距离。
# -OutType：输出文件格式。
# -OutStat：输出文件前缀。
2、可视化LDdeacy
单个群体的可视化：
Plot_OnePop.pl -inFile overlap.all.stat.gz -bin1 500 -bin2 7000 -break 5000 -output overlap
多个群体的可视化（一张图）：
Plot_MultiPop.pl  -inList  Pop.ResultPath.txt -bin1 500 -bin2 7000 -break 5000 -output Fig

# -bin1 -bin2 -break 可以看结果来加大，若结果不太平滑的话
# - List Format :[Pop.ResultPath  PopID ]--制作一个第一列为上面生成群体对应的路径，第二列为群体ID的txt文件

结果展示：


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

推荐阅读更多精彩内容