最全bcftools使用说明--只看本文就够了

Bcftools说明文档

1 简介

BCFtools 是一款多种实用工具的集合,它可以用于处理VCF文件和二进制的BCF文件。它可以接受VCF格式、压缩的VCF格式以及BCF格式,并能自动检测输入的格式类型。在有索引文件存在的条件下,BCFtools 可以应用于所有场景,在没有索引文件存在时,BCFtools只能应用于部分场景。通常,在同时读取多个VCFs时,必须对它们建立索引和压缩。BCFtools 被设计用于处理文件流,它将标准输入文件 "-" 作为标准输入,并将输出作为标准输出。因此可以将多个命令与Unix的管道符进行组合使用。本文档对 BCFtools 大部分的功能命令进行详细的介绍,详细英文教程请参考http://samtools.github.io/bcftools/bcftools.html

2 版本及注意事项

本说明文档最后更新日期为 2019-12-07 与 bcftools v1.10 版本相匹配。需要注意的是bcftools 软件一般倾向于接受 BCF 格式。为了更有效的处理文件,一般首先对 VCF 文件进行 BCF 格式的转换、构建索引操作。

BCF 文件是VCF 文件的二进制压缩格式,采用 bgzip 命令进行压缩。bgzip(Blocked GNU-Zip 格式)是 GZIP 的一种变体,是生信特有的一种压缩格式。BCF、BAM 和 TABIX 格式都采用了 BGZIP 压缩格式。软件的下载请参考 https://sourceforge.net/projects/samtools/files/tabix/

以下为压缩和建立索引操作:
bgzip view.vcf #压缩
bgzip -d view.vcf.gz #解压
bcftools index view.vcf.gz #建立索引

3 BCF1格式

当 Samtools 软件版本小于 0.1.19 时,输出的BCF1格式并不能与此版本的 bcftools 相兼容。为了读取老版本的BCF1文件,可以用老版本的bcftools的view命令先将其转换成VCF文件,再用bcftools软件进行处理。

samtools-0.1.19/bcftools/bcftools view file.bcf1 | bcftools view 

4 变异检测

bcftools 的变异检测来源于samtools mpileup 命令。对于高版本的 bcftools,用户可以选择旧的 samtools 检测变异模型(-c/--consensus-caller参数)或者新的多元检测模型(-m/--multiallelic-caller)。对于大多数任务,多元检测模型较为推荐。

5 bcftools功能列表

当直接运行 bcftools 不输入任何参数时,程序会自动列出参数列表,以下为各个功能的简要介绍:

命令 功能
annotate 编辑VCF文件,添加或去除注释信息
call SNP/INDEL 变异检测,旧版本参数为“view”
cnv CNV 变异检测
concat 将相同样品组的VCF/BCF文件进行连接
consensus 通过 VCF变异创建一致性序列
convert 将 VCF/BCF文件转换为其他格式或其他格式转换成 VCF/BCF格式
csq 单倍体变异检测
filter 用固定的阈值过滤VCF/BCF文件
index 建立VCF/BCF文件索引
isec 对 VCF/BCF 文件取交集、补集
merge 从非重复的样本集中合并 VCF/BCF 文件
mpileup 用 pileup方法检测变异
query 将 VCF/VCF文件转换成用户自定义的格式
reheader 改变 VCF/BCF的表头,改变样品名称
sort 对 VCF/BCF文件进行排序
stats 对 VCF/BCF文件进行统计
view 对 VCF/BCF文件进行取子集、过滤及转换操作

6 命令及参数

6.1 通用命令行参数

以下是 bcftools 软件的通用命令行参数。具体的功能对应的参数见各功能的参数列表。bcftools 的所有功能参数几乎都会用到这些通用命令参数,为了节省篇幅,在解析 bctools 单个功能具体参数时,通用命令行参数不会再重复解析。

command Description
FILE 输入文件可以是VCF、BCF未压缩格式或BGZF压缩格式,"-" 符号作为标准输入
-c,--collapse snps/indels/botah/all/some/none/id 控制如何处理重复位置记录,当有多个输入文件时用此参数定义记录的兼容性,这里的兼容是指工具定义的一致性位点;对于一致性的位点一般只输出第一条记录
none :只有记录的REF和ALT等位基因相同时才认为兼容
some :只有记录的ALT等位基因相一致才认为兼容
all:当位置一致时,所有的记录都是兼容的,不论ALT等位基因是否一致
snps:所有的SNP记录都是兼容的,不论ALT等位基因是否匹配
indels:所有的INDEL记录都是兼容的,无论REF或者ALT等位基因是否匹配
both:"-c indels -c snps" 的简写
id:只有ID列相同才认为是兼容的。该命令只应用于bcftools merge命令
-f,--apply-filters LIST 该参数用于筛选VCF的FILTER列的信息,如果想筛选FILTER列为PASS的记录,可用-f PASS参数。
--no-version 不在输出文件的表头显示 bcftools 版本信息和命令行操作记录
-o,--output FILE 将结果输出到FIEL文件中,不进行标准输出
-O,--output-type b/u/z/v 规定输出文件的格式,b 表示压缩的BCF 格式,u 表示未压缩的BCF 格式,z 表示压缩的VCF 格式,v表示未压缩的VCF 格式
-Ou 参数可以在bcftools管道符命令中使用,添加该参数后可以减少VCF、BCF格式间的转换,加速处理过程
-r,--regions chr/chr:pos/chr:from-to/chr:from-[,...] 逗号分割的区域文件,与 - R 参数类似。不能与 - R 参数联合使用。
-R,--regions-file FILE 区域文件可以在命令行输入,也可以VCF、BED、制表符分隔的文件(默认格式)进行输入。制表符分隔的文件各列分别为:CHORM,POS,POS_TO。软件会默认以1-based对待制表符分隔的文件。BED文件各列分别为:CHORM,POS,POS_TO。软件会以0-based对待BED文件,但BED文件必须以".bed"或者".bed.gz"作为结尾。
以下几点需要注意:染色体编号必须准确输入,“chr20”和“20”是不一致的;使用-R参数需要首先对VCF/BCF文件建立索引;-R参数与-r参数不能同时使用。
-s,--samples [^]LIST 接逗号分隔的样品名列表,用于包含和排除(添加“^”)输入的样品。除了bcftools view命令外,使用-s 参数后,样品结果的输入格式也会按照样品列表名称的顺序进行排列。
-S,--sample-file FILE 接样品名文件,用于包含和排除(添加“^”)输入的样品。每个样品名按照单独一行排列。该参数与-s,--sample参数功能一致,最后输出文件结果也按照输入样品名的顺序排列。
bcftools call 命令可以接受两列文件,第二列可以是倍性(0,1,2)或者性别(“F”,"M")文件。
-t,--targets [^]chr/chr:pos/chr:from-to/chr:from-[,...] 与-r,--regions参数类似,但是该参数访问下一个位置是通过遍历整个VCF/BCF文件实现的,不是根据 tbi/csi 索引文件。
可以同时使用-r 和 -t 参数:-r 参数可以根据索引跳到指定区域,-t 参数可以跳过不在目标中的位置。与 -r 不同,-t 参数可以添加“”前缀来获得补集。例如,“X,Y,MT” 表示跳过X, Y, MT 染色体。-r 和 -t 参数还有一个区别是, -r 参数会检查 indels 的起始和终止位置,-t 参数只会检测起始位置。
需要注意的是-t 和 -T 参数不能同时使用。
-T,--targets-file [^]FILE 与-t,--targets参数类似,但是区域文件是来自于一个文件。
--threads INT 线程数。目前该参数只应用于压缩文件的输出流,只有--output-type为 b 或者 z 时才起作用。默认值为0
(1)-f 参数示例
bcftools view -f PASS,Q20 file.vcf 
(2)-r 参数示例
bcftools view G1_PC_T11_T_uvc1.vcf.gz -f PASS,Q20 --no-version -r 2:16085972-16085990
(3)-R 参数示例
bcftools view -f PASS,Q20 G1_PC_T11_T_uvc1.vcf.gz -R test.bed
(4)-s 参数示例
bcftools view -Ou -s sample1,sample2 file.vcf  | bcftools query -f '%INFO/AC\t%INFO/AN\n'  

6.2 bcftools annotate [OPTIONS] FILE

该命令用于添加或去除注释信息

command Description
-a,--annotations file 接bgzip压缩格式及tabix索引的注释文件。该文件可以是VCF、BED、制表符分隔文件,该文件包含强制列 CHROM 、POS (或者 FROM 和 TO ), 可选列 REF、ALT 及任意数量的注释列。BED文件必须以".bed"或者".bed.gz"作为结尾,否者会默认文件以制表符分隔,按1-based 处理。
当REF 和 ALT 同时存在时,VCF中只有完全匹配的记录才会被注释。当注释文件中有多个ALT等位基因(以逗号形式进行分隔)形式时,至少一个等位基因形式与VCF中的记录保持一致。
-c,--columns list 接逗号分隔的字段列表,目的是从 --annotations file 注释文件中筛选特定的列进行注释替换。
一般情况下,注释文件为 VCF 格式,注释文件的 ID、QUAL、FILTER、INFO、FORMAT 列都可以被编辑。INFO 列的 TAG 可以写作 “INFO/TAG” 或者 TAG;FORMAT 列的 TAG 可以写作 “FORMAT/TAG” 或者 “FMT/TAG”。
如果想完全使用注释文件 INFO 列的全部信息,可以直接使用 -c INFO 参数,如果想使用特定参数信息,可以使用 -c INFO/TAG 参数,如想使用除 TAG 列的其他信息,可以使用 ^INFO/TAG。
默认情况下,在设定注释筛选信息后,原始 VCF 中的信息会被注释文件中的信息所替换,如果不想原始信息被替换,可以使用 +TAG 追加参数。
-e, --exclude EXPRESSIION 接表达式命令,当时表达式返回为真值时,排除该位点。表达式命令和使用见表达式部分
--force 添加该参数可以忽略运行中的错误,但是有可能会导致生成BCF文件的异常,使用该参数需要最后彻底的核查
-I, --set-id [+] FORMAT 动态的指定ID列。在默认情况下会替换现有的ID列。如果添加“+”命令,只会替换缺失的ID列,如下所示:
<pre>bcftools annotate --set-id +'%CHROM_%POS_%REF_%FIRST_ALT' file.vcf</pre>
-i, --include EXPRESSION 接表达式命令,当时表达式返回为真值时,输出该位点。表达式命令和使用见表达式部分
-k, --keep-sites 保留没有通过 -i 和 -e 表达式过滤条件的记录,而不是直接进行抛弃
-l, --merge-logic tag:'first'/append/unique/sum/avg/min/max[,…] 该选项定义了在一个区域存在多个注释值时,如何进行处理。first (表示使用第一个值); append (表示进行追加);unique (表示去除重复的注释); sum (进行加和,这个只适用于数值型);avg (取平均值);min (取最小值);max(取最大值)
-m, --mark-sites TAG 添加该参数会在结果中增加一个新的 TAG 标记,用于注释记录在 -a 参数的输入文件是否存在,存在标记为“+”,不存在标记为“-”
-x, --remove list 要删除的注释列表。使用 "FILTER" 删除所有 filters 注释信息或使用 “FILTER/SomeFilter” 删除FILTER列特定的注释信息。与此相似,“INFO” 可以删除所有 "INFO" 列标签,“FORMAT” 可以删除所有除 GT 外的所有标签信息。删除整列信息后,会用 “.” 进行占位

bcftools annotate 示例:

# 从src.bcf注释文件中,筛选ID、QUAL、INFO/TAG 列的注释信息对 dst.bcf 文件进行注释,INFO/TAG 参数信息不是进行替换,而是进行追加
bcftools annotate -a src.bcf -c ID,QUAL,+TAG dst.bcf

# 注释 INFO 列及除 FORMAT/GT 以外的 FORMAT 列的信息
bcftools annotate -a src.bcf -c INFO,^FORMAT/GT dst.bcf

# 删除 VCF文件中 ID 、INFO/DP、FORMAT/DP 列的信息
bcftools annotate -x ID,INFO/DP,FORMAT/DP file.vcf.gz

# 删除 VCF文件中 INFO 列及除 FORMAT/GT、FORMAT/PL 以外的所有 FORMAT 列的信息
bcftools annotate -x INFO,^FORMAT/GT,FORMAT/PL file.vcf<br>

6.3 bcftools call [OPTIONS] <in.vcf.gz>

该命令用于 SNP/INDEL 的检测,输入文件为 VCF/BCF 格式。该命令常与 bcftools mpileup 命令连用。取代了之前的 bcftools view caller 命令。

command Description
-A, --keep-alts 输出该位置的所有等位基因类型,不进行质量过滤
-f, --format-fields list 输出 FORMAT 类型,目前只有GP,GQ 类型备选,具体解释参照VCFv4.2文档
-F, --prior-freqs <AN,AC> 添加千人基因组计划的等位基因频率
-G, --group-samples FILE/- 默认所有的样品来自同一种群。该参数可以将样品设置为不同的种群,以种群为单位进行处理
-V, --skip-variants snps/indels 忽略snps 或 indels位点
-v, --variants-only 只输出有变异的位点
-c, --consensus-caller samtools/bcftools 旧版变异检测参数,与 -m 参数冲突
-m, --multiallelic-caller bcftools 优化后的变异检测参数,与 -c 参数冲突。参数背后的模型进行了优化,较为推荐

bcftools call 示例:

# 利用bcftools mpileup 生成初始VCF文件,对 VCF 文件进行变异检测,输出所有变异位点 
bcftools mpileup -Ou -f reference.fa alignments.bam | bcftools call -mv -Ob -o calls.bcf

6.4 bcftools cnv [OPTIONS] FILE

该命令用于检测拷贝数变异。与其他用于检测单样品或者混合样品(正常细胞和肿瘤细胞混合)拷贝数变化的软件不同,该软件特定用于检测两个细胞系之间拷贝数的差异。该命令的输入文件为 VCF 文件,该VCF 文件需要 BAF(B-allele frequency) 和 LRR(Log R Ratio intensity)注释信息。利用 HMM(hidden Markov models) 模型,软件可以给出以下拷贝数状态:CN 2 (normal), 1 (single-copy loss), 0 (complete loss), 3 (single-copy gain)

# BAF 和 LRR 计算方法<br>LRR ~ fold change i.e. log2 depth sample / depth other samples
BAF = (ALT depth)/(Total depth)<br>
#输入 VCF 文件,格式示例
##fileformat=VCFv4.2
##FORMAT=(ID=GT,Number=1,Type=String,Description="Genotype")
##FORMAT=(ID=GC,Number=1,Type=Float,Description="GenCall score")
##FORMAT=(ID=IA,Number=1,Type=Float,Description="Intensity of the A Allele")
##FORMAT=(ID=IB,Number=1,Type=Float,Description="Intensity of the B Allele")
##FORMAT=(ID=BAF,Number=1,Type=Float,Description="B Allele Frequency")
##FORMAT=(ID=LRR,Number=1,Type=Float,Description="Log R Ratio")
#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  test
10      135656  rs10904561      T       G       .       .       .       GT:LRR:BAF:IA:IB:GC     0/1:0.0415:0.5176:0.543:0.661:0.9477
10      135708  rs7917054       G       A       .       .       .       GT:LRR:BAF:IA:IB:GC     1/0:0.2441:0.501:0.742:0.659:0.9604
10      178434  rs7089889       T       G       .       .       .       GT:LRR:BAF:IA:IB:GC     0/1:0.1440:0.4833:0.836:0.827:0.8884
10      188805  rs12146291      A       G       .       .       .       GT:LRR:BAF:IA:IB:GC     0/1:0.2144:0.4715:0.776:0.743:0.9610
10      203471  rs10903451      A       G       .       .       .       GT:LRR:BAF:IA:IB:GC     0/0:0.4201:0.0000:0.619:0.022:0.9337
10      252693  rs2379078       T       C       .       .       .       GT:LRR:BAF:IA:IB:GC     0/0:-0.2673:0.0068:0.806:0.081:0.9018
10      273344  rs2448378       A       C       .       .       .       GT:LRR:BAF:IA:IB:GC     0/0:-0.7407:0.0000:0.497:0.012:0.9411
10      279248  rs2496278       C       A       .       .       .       GT:LRR:BAF:IA:IB:GC     1/0:-0.4228:0.6268:0.285:0.958:0.6772
command Description
-c, --control-sample string 对照样品名称
-s, --query-sample 查询样品名称
-f, --AF-file file 等位基因频率文件,含有以下列 CHR,POS,REF,ALT,AF
-o, --output-dir path 结果输出目录

bcftools cnv 示例:

# 利用bcftools cnv 生成拷贝数变化过程文件,对拷贝数变化进行绘图 
bcftools cnv -o cnv/ outdir/test.vcf.gz
python cnv/plot.test.py

6.5 bcftools concat [OPTIONS] FILE1 FILE2 […]

该命令用于连接或合并多个 VCF/BCF 文件。所有的输入文件必须含有相同的样品列,并且样品列的排列顺序要一致。该命令可以用于将多个以染色体为单位的 VCF 文件合并成一个 VCF 文件,或者将不同突变类型的SNP VCF、INDEL VCF进行合并。输入的 VCF文件的必须按照染色体和位置进行排序。如果想输出排序后的结果,一般需要添加 -a, --allow-overlaps 参数,否则该命令只是是简单的进行合并。如果使用 --naive 参数,输入文件合并时会不进行解压缩操作,速度会非常快,该参数不能与-a, --allow-overlaps 参数联用。

command Description
-a, --allow-overlaps 在输出结果中,突变的次序按照染色体、位置信息排序,不按照输入文件的突变排序
-d, --rm-dups snps/indels/both/all/none 去除指定类型的重复记录,当记录有重复时只输入一种记录信息
-D, --remove-duplicates 与 -d none 参数效果相同,去除所有重复记录
-f, --file-list FILE 从FILE文件中读取输入文件信息,每个文件占据一行
-n, --naive 合并 VCF/BCF 文件并不进行解压缩操作。使用该参数时需要输入文件格式类型相同,VCF头文件相同。在合并过程中,程序会根据表头检查所有合并文件的兼容性,如果发现异常,程序会抛出错误。
--naive-force 与 --naive 参数作用一致,但是不会进行兼容性检测,使用该参数有一定的危险性

bcftools concat 示例:

# 合并SNP 、INDEL类型信息,并去除重复记录 
bcftools concat -a snps.vcf.gz indels.vcf.gz -D -Ob -o concat.vcf.gz

6.6 bcftools consensus [OPTIONS] FILE

该命令可以将 VCF 文件中的突变应用于参考基因组的 fasta 序列,创建有对应突变的参考基因组序列文件。在默认情况下,该程序可以将 ALT 变异替换基因组序列碱基,获得最终序列信息。

command Description
-c, --chain FILE 写出链文件的长度信息
-f, --fasta-ref FILE 参考基因组序列,fasta 格式
-H, --haplotype 1/2/R/A/LR/LA/SR/SA/1pIu/2pIu 根据 FORMAT/GT 信息,筛选特定等位基因生成一致序列:1:GT中的第一个等位基因,忽略相位
2:GT中的第二个等位基因,忽略相位
R:杂合子突变中的REF型等位基因
R:杂合子突变中的ALT型等位基因
LR, LA:碱基序列较长的等位基因
SR, SA:碱基序列较短的等位基因
1pIu,2pIu:等位基因相位已知时输出第一、第二个等位基因,相位未知时输出简并碱基
-I, --iupac-codes 输出简并碱基类型序列
-m, --mask FILE BED 格式文件或以 TAB 分隔的文件。在此区内的碱基会被替换为N
-M, --missing CHAR 如果 genotypes类型缺失,不会进行忽略,而是以 CHAR 字符进行替换
-s, --sample NAME 应用突变碱基替换的样品名称

bcftools consensus 示例:

# 对 “NA001” 样品的突变进行碱基替换,突变碱基替换为简并碱基形式,获得最终参考基因组序列文件
bcftools consensus -I -s NA001 -f in.fa in.vcf.gz > out.fa

# 对参考基因的特定区域进行突变碱基替换,获得最终参考基因组序列文件
samtools faidx ref.fa 8:11870-11890 | bcftools consensus in.vcf.gz -o out.fa

6.7 bcftools convert [OPTIONS] FILE

该命令可以实现 VCF/BCF 格式与其他格式的相互转换。

command Description
-G, --gensample2vcf prefix or gen-file,sample-file 将 IMPUTE2 转化为 VCF格式。第二列必须为"CHROM:POS_REF_ALT"格式
-g, --gensample prefix or gen-file,sample-file 将 VCF 格式转换为 gen/sample 格式。gen/sample 格式如下:
<pre> .gen
----
1:111485207_G_A 1:111485207_G_A 111485207 G A 0 1 0 0 1 0
1:111494194_C_T 1:111494194_C_T 111494194 C T 0 1 0 0 0 1
.sample
-------
ID_1 ID_2 missing
0 0 0
sample1 sample1 0
sample2 sample2 0</pre>
--gvcf2vcf 将 gVCF 格式转换为 VCF格式。将gVCF中的REF blocks 扩展成位点
-f, --fasta-ref file 参考基因组序列,fasta 格式
--tsv2vcf file 将Tab键分隔的变异格式(例如 23andMe 格式)转换为VCF格式。23andMe 格式如下:
<pre>rs6139074 20 63244 AA
rs1418258 20 63799 CC
rs6086616 20 68749 TT
rs6039403 20 69094 AG
-c, --columns list 逗号分隔的类型列表,对 --tsv2vcf 参数的输入文件各列信息进行解释,目前支持ID,CHROM,POS,AA类型,AA的第一个字符代表参考基因组,第二个字符代表二倍体的另一个等位基因

bcftools convert 示例:

# 将 VCF 格式转换为 gen/sample格式
bcftools convert -g - ch1.vcf.gz -Ob -o sample.vcf

# 将 gVCF 格式转换为 VCF格式
bcftools convert --gvcf2vcf  --fasta-ref ref.fa gvcf.vcf -o test.vcf

# 将 23andme 结果转换成 VCF格式
bcftools convert -c ID,CHROM,POS,AA -s SampleName -f 23andme-ref.fa --tsv2vcf 23andme.txt -Oz -o out.vcf.gz

6.8 bcftools csq [OPTIONS] FILE

该命令利用单倍体预测算法(Haplotype aware consequence predictor),可以精准的注释复杂类型变异,获得正确的变异注释结果。bcftools csq 命令需要输入参考基因组文件、GFF3格式注释文件、需要注释的VCF文件,最终注释信息会以INFO/BCSQ形式进行展现。bcftools csq 算法会以单倍体为单位综合考虑复杂变异所处密码子的位置,获取正确的注释信息,如果只是单独考虑各个突变位点,注释结果会产生较大的错误,例如:终止突变被解读为普通氨基酸变化,普通的氨基酸变化当做移码突变。如下图所示:


csq.png
command Description
-f, --fasta-ref FILE 参考基因组序列,fasta 格式
-g, --gff-annot FILE GFF3 格式的注释文件,下载可参考 ftp://ftp.ensembl.org/pub/current_gff3/homo_sapiens

bcftools csq 示例:

# 对 VCF 文件进行注释
bcftools csq -f ref.fa -g Homo_sapiens.GRCh38.103.chromosome.1.gff3.gz ch1.vcf.gz -Ob -o out.bcf

# BCSQ信息以 “|” 进行分隔,分别为:变异类型、基因名称、转录本名称、基因类型、正负链信息、氨基酸变化、碱基变化
BCSQ=missense|PLCH2|ENST00000278878|protein_coding|+|338V>338L|2489790G>C

6.9 bcftools filter [OPTIONS] FILE

该命令用于固定阈值过滤

command Description
-e, --exclude EXPRESSIION 接表达式命令,当时表达式返回为真值时,排除该位点。表达式命令和使用见表达式部分
-g, --SnpGap INT 过滤 indel 附近的 SNPs 位点。下面的例子展示了 --SnpGap 3 参数删除INDEL附近3bp内的SNPs位点:
<pre>4号位置发生了一个del,1号、7号位置的SNPs被过滤掉;0号、8号位置的SNPs没有被过滤掉: 0123456789
ref .G.GT..G..
del .A.G-..A..
4号位置发生了一个ins,1号、6号位置的SNPs被过滤掉;0号、7号位置的SNPs没有被过滤掉:
0123-456789
ref .G.G-..G..
ins .A.GT..A..
</pre>
-G, --IndelGap INT 过滤 indels 聚集区的第 INT 个变异,只保留一个突变。下面的例子展示了 --INdelGap 2 参数删除 INDELs 聚集区的突变位点:
<pre>第二个indel被过滤掉:
012345678901
ref .GT.GT..GT..
del .G-.G-..G-..
相似的,第二个indel被过滤掉:
01 23 456 78
ref .A-.A-..A-..
ins .AT.AT..AT..
</pre>
-i, --include EXPRESSION 接表达式命令,当时表达式返回为真值时,输出该位点。表达式命令和使用见表达式部分
-m, --mode [+x] FILTER annotations 存在时行为模式设置。默认FILTER列缺失时,将其设置为PASS。“+” 表示不替换,追加新的 FILTER 字符信息。“x” 表示重设过滤位点
-s, --soft-filter STRING/+ 满足筛选条件时将其注释为特定字符

bcftools filter 示例:

# 将QUAL值小于10的FILTER列注释为 lowQual
bcftools filter view.vcf -s lowQual -e 'QUAL<10' -Ob -o out.bcf

6.10 bcftools index [OPTIONS] in.bcf | in.vcf.gz

该命令用于对 bgzip 压缩的 VCF/BCF 文件建立索引,方便获取随机位置信息。该命令默认会产生 CSI (位置排序信息) 格式。CSI 格式的索引文件所支持的染色体长度阈值是2^31。使用 bcftools index 的 -t 参数可以产生 TBI 格式的索引文件,TBI 格式索引文件支持染色体长度的阈值是2^29。当启用索引文件时,bcftools 软件会先搜索 CSI 格式,再搜索 TBI 格式。

command Description
-c, --csi 生成 CSI 格式的索引文件,该参数为默认参数
-f, --force 当已存在索引文件时,强制进行覆盖
-m, --min-shift INT 设置索引文件的最小区间大小,默认为2^14
-o, --output FILE 设置输出文件的名称。如果不进行设定,默认在输入文件的末尾添加 .csi 或 .tbi 后缀
-t, --tbi 生成 TBI 格式的索引文件
-n, --nrecords 根据CSI 或 TBI 索引文件打印记录的个数
-s, --stats 打印索引文件构建过程中的统计信息

bcftools index 示例:

# 对 VCF 文件构建索引,生成 .csi 格式的索引文件
bcftools index view.vcf

6.11 bcftools isec [OPTIONS] A.vcf.gz B.vcf.gz […]

该命令用于在多个 VCF 文件之间创建交集、差集。根据参数的不同,程序可以从一个或多个文件中输出这些记录。

command Description
-C, --complement 输出只在第一个文件中存在,不在其他文件中存在的变异位点
-e, --exclude EXPRESSIION 接表达式命令,当时表达式返回为真值时,排除该位点。
如果 -e 或者 -i 参数只出现一次,表示针对所有输入文件都采用相同的过滤规则。否则每个输入文件都需要 -e 或 -i 参数提供一套过滤规则。如果输入文件不需要进行过滤,表达式用“-”参数进行替代。使用示例如下:
表达式命令和使用见表达式部分。
-i, --include EXPRESSION 接表达式命令,当时表达式返回为真值时,输出该位点。使用规则与 -e 参数相同
表达式命令和使用见表达式部分
-n, --nfiles [+-=]INT/~BITMAP -n=2 : 表示在 2 个输入文件中出现;-n-2 :表示在2个或小于2个文件中出现;-n+2 :表示在 2 个或者大于 2 个文件中出现;-n~1100 :表示在第1、第 2 个输入文件中出现,不在第3、第4个文件中出现。1 表示存在,0 表示缺失
-p, --prefix DIR 添加目录参数时,在目录中输出多个文件的交集和补集。多与 -w 参数连用
-w, --write LIST 后跟多个 VCF 文件列表。-w1 : 只输出第1个文件的交集和补集

bcftools isec 示例:

# 创建两个文件的交集和补集将其保存在 dir文件夹<br>bcftools isec -p dir A.vcf.gz B.vcf.gz<br>
# 对A文件(INFO/MAF>=0.01)和B文件(INFO/dbSNP)进行过滤,但不对C文件进行过滤;在 dir 文件夹生成 3 个文件的交集位点,这些位点至少在 2 个输入文件中出现
bcftools isec -e'MAF<0.01' -i'dbSNP=1' -e- A.vcf.gz B.vcf.gz C.vcf.gz -n +2 -p dir<br>
# 在 dir 目录输出A文件和B文件的交集,只按照A文件的格式输出
bcftools isec -p dir -n=2 -w1 A.vcf.gz B.vcf.gz<br>
# 在 dir 目录输出A文件和B文件的差集
bcftools isec -p dir -n-1 -c all A.vcf.gz B.vcf.gz<br>
# 打印出只在A、B文件中出现,不在C、D文件中出现的位点
bcftools isec -n~1100 -c all A.vcf.gz B.vcf.gz C.vcf.gz D.vcf.gz

6.12 bcftools merge [OPTIONS] A.vcf.gz B.vcf.gz […]

该命令用于将多个没有重复样品的VCF/BCF文件合并成一个多样品的VCF/BCF文件。如果合并的文件中有重复性的样品名称,程序会自动退出,除非添加 --force-samples 参数,添加参数后会在重复样品名称前添加标记进行输出。例如 A.vcf.gz 含有的样品为S1、S2、S3, B.vcf.gz含有的样品为S3、S4,输出入结果中会包含4个样品名称,分别是S1、S2、S3、2:S3 和 S4。

command Description
--force-samples 当合并的文件中有重复名称时,强制进行合并。添加参数后会在重复样品名称前添加索引标记进行输出,如2:S3
--print-header 只输出合并后的表头文件
--use-header FILE 自定义表头文件,进行输出
-0 --missing-to-ref 假设缺失位点的基因型为 0/0
-i, --info-rules TAG:METHOD[,…] 规定 INFO 列信息合并规则。合并规则有 sum,avg,min,max,join 。例如,可以用 -i tDP:sum 将所有输入文件的 tDP 进行合并,作为最终输出。
如果没有指定合并规则,则默认从第一个输入文件中获取该值。
合并的 QUAL 值目前设置为输入文件中的最大值,用户无法控制。
-l, --file-list FILE 从FILE文件中指定输入合并文件,每行设定一个输入文件
-m, --merge snps/indels/both/all/none/id 该选项可以控制哪些类型可以创建多等位基因记录,默认参数为 both。
<pre>-m none #不输出多等位信息
-m snps #SNP位点输出多等位信息
-m indels #INDEL位点输出多等位信息
-m both #SNP和INDEL位点都输出多等位信息
-m all #SNP位点可与INDEL位点进行合并输出多等位信息
-m id #根据ID列进行等位基因的合并</pre>

bcftools merge 示例:

# 合并 A、B 两个 VCF 文件的变异,不在一个位置输出多个等位基因
bcftools merge -m none  A.vcf.gz B.vcf.gz -Ou -o AB.vcf.gz 

6.13 bcftools mpileup [OPTIONS] -f ref.fa in.bam [in2.bam […]]

该命令可以对输入的 BAM/CRAM 文件进行初步的变异检测,生成 VCF/BCF 文件。该命令生成的初始 VCF/BCF 文件记录了所有位点的 reads 支持情况,包括比对错误和测序错误,要获得真实的 SNV/INDEL变异,需与bcftools call 配合使用。bcftools mpileup 命令来源于 samtools mpileup 命令(-v 参数产生 VCF文件;-g 参数产生 BCF 文件),但是不能产生samtools mpileup 命令默认的 pileup 格式文件。bcftools 软件添加 mpileup 命令是为了避免 samtools 软件和 bcftools 软件版本出现兼容问题。

command Description
-6, --illumina1.3+ 假定质量值以 illumina1.3+ 编码
-A, --count-orphans 变异检测过程中不忽略异常序列
-b, --bam-list FILE 批量输入 BAM 文件,每个 BAM 文件占据一行
-B, --no-BAQ 不计算 BAQ(比对质量值),使用该参数可以大大减少 SNPs 的假阳
-C, --adjust-MQ INT 当 reads 有较多错配时,使用该参数可以降低比对质量值
-d, --max-depth INT 各个变异位点,单个输入文件 reads 覆盖深度的最大值,默认值为 250。设置参数主要目的是防止输入文件较多,内存过大
-E, --redo-BAQ 动态计算 BAQ 值,忽略现有的 BQ 值
-f, --fasta-ref FILE FASTA 格式的参考基因组文件,需提前建立索引。基因组文件可以用 bgzip 命令进行压缩。
--no-reference 不输入参考基因组文件,使用该参数时,可以不输入 --fasta-ref 参数
-G, --read-groups FILE 输出结果列的样品名称列表,每个样品名称占据一行。当添加“^”前缀时,则不输出列表中的样品。该文件还可以为样品重新命名。输入文件示例如下:
<pre>RG_ID_1 #样品名称
RG_ID_2 SAMPLE_A #样品名称重新命名
</pre>
-q, -min-MQ INT 用于变异检测时的最小比对质量值,默认为0
-Q, --min-BQ INT 用于变异检测时的最小碱基质量值,默认为13
--ignore-RG 忽略 RG tags。每个BAM文件,当做一个样品
-x, --ignore-overlaps 不对重复区域进行变异检测

** bcftools mpileup 示例:**

# 检测 SNV/INDEL 并将QUAL值大于20,深度过高(DP >100)位点的 FILTER 列标记为 lowQUAL
bcftools mpileup -Ou -f ref.fa aln.bam | \
bcftools call -Ou -mv | \
bcftools filter -s LowQual -e '%QUAL<20 || DP>100' > var.flt.vcf

6.14 bcftools query [OPTIONS] file.vcf.gz [file.vcf.gz […]]

该命令可从 VCF/BCF 文件中提取特定字段信息,并按照自定义格式进行输出

command Description
-e, --exclude EXPRESSIION 接表达式命令,当时表达式返回为真值时,排除该位点。表达式命令和使用见表达式部分
-f, --format FORMAT 设定的数据结果模式类型,具体可参见示例部分
-H, --print-header 打印出表头
-i, --include EXPRESSION 接表达式命令,当时表达式返回为真值时,输出该位点。表达式命令和使用见表达式部分
-l, --list-samples 列出样品名称并退出
-u, --allow-undef-tags 如果 FORMAT 列有未定义的 tags 出现,不抛出错误,用 “.” 进行替换
-v, --vcf-list FILE 处理多个 VCF 文件时,可以将其放入一个文件进行处理

bcftools query 示例:

# -f参数,设定的格式示例:
%CHROM    代表 CHROM 列,与POS, ID, REF, ALT, QUAL, FILTER 等列用法相同
%END      代表 REF 等位基因的终止位置
%END0     代表 REF 等位基因的终止位置,0-based格式,比 %END 等位基因位置小1
%FORMAT   代表输出整个 FORMAT 列及相应样本的信息,常与 -s 或者 -S 参数联用
%GT       代表输出 GT 信息
%INFO     代表输出整个 INFO 列信息
%INFO/TAG 代表输出 INFO 列的TAG参数信息,例如 %INFO/tFA
%TAG{INT} 代表输出TAG 参数的第几个字段,例如 %ALT{0} 代表 ALT 的第1字段 
%POS0     代表输出 POS 的位置,0-based格式,比 %POS 等位基因位置小1
%SAMPLE   代表样品名称
%TYPE     代表突变类型,程序会自动检测突变类型并输出,类型有REF, SNP, MNP, INDEL, BND, OTHER
[]        []内放置 FORMAT 列的信息,例如 [%DP] 代表输出 FORMAT/DP 信息
\n        代表换行符
\t        代表制表分隔符

# 输出 VCF 文件中 CHORM,POS,REF,ALT 列信息
bcftools query -f '%CHROM  %POS  %REF  %ALT{0}\n' file.vcf.gz

# 输出 VCF 文件中 CHORM,POS,REF,ALT,FOMAT/GT 列信息,并以 制表符进行分隔
bcftools query -f '%CHROM\t%POS\t%REF\t%ALT[\t%SAMPLE=%GT]\n' file.vcf.gz

# 输出 VCF 文件中 FORMAT/GT 和 FORMAT/GT 列信息
bcftools query -f 'GQ:[ %GQ] \t GT:[ %GT]\n' file.vcf

# 筛选 VCF 文件中 GT 类型为 alt 的突变,输出 CHORM,POS,SAMPLE,FOMAT/GT 列信息
bcftools view -i'GT="alt"' file.bcf -Ou | bcftools query -f'[%CHROM:%POS %SAMPLE %GT\n]'

6.15 bcftools reheader [OPTIONS] file.vcf.gz

该命令用于更改 VCF/BCF 文件的表头和样品名称

command Description
-h, --header FILE 输入的 VCF 表头文件
-o, --output FILE 输出文件名称
-s, --samples FILE 新样品名称文件,每个样品名称占一行

bcftools reheader 示例:

# 获取 VCF 文件的表头文件
bcftools view -h file1.vcf >header.txt

# 替换 VCF 文件的表头文件
bcftools reheader -h header.txt file2.vcf -o new.header.vcf

6.16 bcftools sort [OPTIONS] file.bcf

该命令用于对 VCF/BCF 文件进行排序

command Description
-m, --max-mem <float>[kMG] 最大内存阈值
-o, --output FILE 输出文件名称
-T, --temp-dir <dir> 临时文件夹,用于存放临时文件,生成排序文件后会自动进行删除

bcftools sort 示例:

# 对 VCF 文件进行排序
bcftools sort -m 1G file.vcf.gz -o file.sort.vcf.gz -Ou -T ./temp

6.17 bcftools stats [OPTIONS] A.vcf.gz [B.vcf.gz]

该命令会对 VCF/BCF 文件的变异进行统计,生成文本格式的统计文件。当输入两个 VCF 文件时,程序会分别统计两个文件差集的变异结果、两个文件交集的变异结果。生成的统计文件可以用 plot-vcfstats 脚本进行绘图,脚本下载见 https://github.com/samtools/bcftools/blob/develop/misc/plot-vcfstats

command Description
-1, --1st-allele-only 当一个位点有多个等位基因时,只统计第一个等位基因
-d, --depth INT,INT,INT 设置深度阈值的分布:最小值,最大值,bin大小 默认为[0,500,1]
-I, --split-by-ID 按照 ID 列统计变异分析

bcftools stats 示例:

# 对 VCF 文件进行排序
bcftools stats file.vcf >stat.txt

6.18 bcftools view [OPTIONS] file.vcf.gz [REGION […]]

该命令可以根据位点和筛选规则,对 VCF/BCF 文件进行查看、过滤、取子集等操作。也可以用于 VCF 和 BCF 之间的相互转换。旧版本之前的命令为 bcftools subset 。

command Description
-G, --drop-genotypes 与 -s 参数连用,不输出 genotype 和 FORMAT 列信息
-h, --header-only 只输出 VCF 的表头文件
-H, --no-header 不输出 VCF 的表头文件
-l, --compression-level [0-9] 输出文件的压缩等级。0 代表不进行压缩,1 代表速度最快的压缩,9代表压缩效果最好、
-a, --trim-alt-alleles ALT 列有多个等位基因时,进行去除
-c, --min-ac INT[:nref/:alt1/:minor/:major/:'nonmajor'] 达到最小等位基因位置阈值(INFO/AC)的位点会被打印输出
-C, --max-ac INT[:nref/:alt1/:minor/:'major'/:'nonmajor'] 达到最大等位基因位置阈值(INFO/AC)的位点会被打印输出
-e, --exclude EXPRESSIION 接表达式命令,当时表达式返回为真值时,排除该位点。表达式命令和使用见表达式部分
-g, --genotype [^][hom/het/miss] 用于包含和排除(添加“^”)纯合(hom)、杂合(het)、缺失(miss)基因型
-i, --include EXPRESSION 接表达式命令,当时表达式返回为真值时,输出该位点。表达式命令和使用见表达式部分
-k, --known 只打印输出 ID 列不为 “.” 的位点
-m, --min-alleles INT 最小等位基因个数阈值
-M, --max-alleles INT 最大等位基因个数阈值
-n, --novel 只打印输出 ID 列为 “.” 的位点
-q, --min-af FLOAT[:nref/:alt1/:minor/:major/:nonmajor] 达到最小等位基因频率阈值(INFO/AC)的位点会被打印输出
-Q, --max-af FLOAT[:nref/:alt1/:minor/:major/:nonmajor] 达到最大等位基因频率阈值(INFO/AC)的位点会被打印输出
-v, --types snps/indels/mnps/other 输入逗号分隔的变异类型列表,筛选输出特定类型的变异
-V, --exclude-types snps/indels/mnps/ref/bnd/other 输入逗号分隔的变异类型列表,筛选不输出特定类型的变异

bcftools view 示例:

# 筛选满足深度、频率要求的snps,indels型突变
bcftools view -v snps,indels -i 'INFO/tFA>0.2 && FORMAT/DP < 200 ' file.vcf.gz -Ou -o outfile.vcf.gz

# 筛选 2 号染色体上,满足DP、GQ 要求的,并且 FILTER 列为 PASS 的位点
bcftools view -r 2  -i 'FMT/DP>10 & FMT/GQ >20 '  -f PASS  file.vcf.gz -Ou -o outfile.vcf.gz

7 表达式(EXPRESSIONS)

bcftools 大部分命令都接受过滤表达式,有效的过滤表达式需要包含以下信息:

(1) 数值常量,字符串常量,样品名称(目前只支持在过滤 ID 列时使用)<br>1,1.0,1e-4<br>"string"<br>@file_name<br>
(2) 数学运算符<br> +,*,-,/<br>
(3) 比较运算符
== (same as =), >, >=, <=, <, !=<br>
(4) 正则运算符,例如 “~” ,“!~” 正则表达式是大小写敏感的,除非添加 “/i” 参数
INFO/HAYSTACK ~ "needle"
INFO/HAYSTACK ~ "NEEDless/i"<br>
(5) 圆括号
(, )<br>
(6) 逻辑运算符
&&,  &, ||,  |<br>

**&& 与 & 运算符及 || 与 | 运算符的区别**
以&& 与 & 进行介绍:
<1> 当VCF文件中含有多个样品时, & 表示每个样品的参数都需要满足共同的条件:
$ bcftools query -i'FMT/DP>10 & FMT/GQ>20' -f'%POS[\t%SAMPLE:DP=%DP GQ=%GQ]\n' file.bcf
49979   SampleA:DP=10 GQ=50     SampleB:DP=20 GQ=40

<2> 当VCF文件中含有多个样品时, && 表示不必在每个样本中都满足共同的条件,只需在所有样本中共同满足条件就可以:
$ bcftools query -i'FMT/DP>10 && FMT/GQ>20' -f'%POS[\t%SAMPLE:DP=%DP GQ=%GQ]\n' file.bcf
31771   SampleA:DP=10 GQ=50     SampleB:DP=40 GQ=20
49979   SampleA:DP=10 GQ=50     SampleB:DP=20 GQ=40


(7) INFO 标签, FORMAT 标签,以及列名
INFO/DP or DP
FORMAT/DV, FMT/DV, or DV
FILTER, QUAL, ID, CHROM, POS, REF, ALT[0]<br>
(8) 1(或者 0 ),用来测试标签的存在(或者缺失)
FlagA=1 && FlagB=0<br>
(9) “.” 用来测试缺失值
DP=".", DP!=".", ALT="."<br>
(10) 使用这些表达式 (".|.", "./.", "."),无论倍性如何,可以匹配缺失的基因型
GT=".|.", GT="./.", GT="."<br>
(11) 样本的基因型:参考型、突变型、纯合型、杂合型、单倍体、二倍体
GT="ref"
GT="alt"
GT="mis"
GT="hom"
GT="het"
GT="hap"
GT="RR"
GT="AA"
GT="RA" or GT="AR"
GT="Aa" or GT="aA"
GT="R"
GT="A"<br>
(12) 根据REF,ALT列筛选筛选符合要求的变异类型(indel,snp,mnp,ref,bnd,other,overlap),其中正则表达式"\~"表示列出至少一种等位基因, "="表示列出所有等位基因
TYPE="snp"
TYPE~"snp"
TYPE!="snp"
TYPE!~"snp"<br>
(13) 表达式中的数组是以 0-based 为起始的,“*” 表示任何元素,“-” 表示范围。对于查询FORMAT向量中的元素,“:” 可以用于筛选向量中的一个元素和样品。具体示例如下:
INFO/AF[0] > 0.3             .. first AF value bigger than 0.3
FORMAT/AD[0:0] > 30          .. first AD value of the first sample bigger than 30
FORMAT/AD[0:1]               .. first sample, second AD value
FORMAT/AD[1:0]               .. second sample, first AD value
DP4[*] == 0                  .. any DP4 value
FORMAT/DP[0]   > 30          .. DP of the first sample bigger than 30
FORMAT/DP[1-3] > 10          .. samples 2-4
FORMAT/DP[1-]  < 7           .. all samples but the first
FORMAT/DP[0,2-4] > 20        .. samples 1, 3-5
FORMAT/AD[0:1]               .. first sample, second AD field
FORMAT/AD[0:*], AD[0:] or AD[0] .. first sample, any AD field
FORMAT/AD[*:1] or AD[:1]        .. any sample, second AD field
(DP4[0]+DP4[1])/(DP4[2]+DP4[3]) > 0.3
CSQ[*] ~ "missense_variant.*deleterious"<br>
(14) 在存在多个样品时,可以提供一份多个样品名称的文件(每个行一个样品)可以更有效的进行统计
GT[@samples.txt]="het" & binom(AD)<0.01<br>
(15) 函数标签名称,可以作用于 FORMAT 标签和 INFO 标签,例如求最大值、最小值、平均数 ...
MAX, MIN, AVG, MEAN, MEDIAN, STDEV, SUM, STRLEN, ABS, COUNT<br>
(16) 双尾二项分布函数
binom(FMT/AD)                .. GT can be used to determine the correct index
binom(AD[0],AD[1])           .. or the fields can be given explicitly
phred(binom())               .. the same as binom but phred-scaled<br>
(17) 变量名称(如果不存在可以动态进行计算),包括以下变量名:N_ALT(突变等位基因数量);N_SAMPLES(样本数量);<br>AC(突变等位基因数量);MAC(最小等位基因数量);AF(突变等位基因频率 AF=AC/AN);MAF(最小突变等位基因频率 MAF = MAC/AN);<br>AN(等位基因数量);N_MISSING(基因型缺失样本数量);F_MISSING(基因型缺失样本比例);ILEN(插入缺失长度,缺失为负值,插入为正值)
N_ALT, N_SAMPLES, AC, MAC, AF, MAF, AN, N_MISSING, F_MISSING, ILEN<br>
(18) N_PASS(符合表达式条件的样品个数),F_PASS(符合表达式条件的样品比例)
N_PASS(GQ>90 & GT!="mis") > 90
F_PASS(GQ>90 & GT!="mis") > 0.9<br>
(19) 自定义的perl脚本过滤
perl:path/to/script.pl; perl.severity(INFO/CSQ) > 3<br>
表达式注意事项:
<1> 字符串的比较和正则表达式是不区分大小写的
<2> 变量名称和函数名称是不区分大小,但是标签名称区分。例如“qual”可以用来代替“QUAL”,“strlen()”可以代替“STRLEN()”,但“dp”不能代替“DP”
<3> 当查询多个值时,程序将测试所有元素,并在结果上使用"或"逻辑。例如,当查询“TAG=1,2,3,4”,会进行如下判断:<br>
-i 'TAG[*]=1'   .. true, the record will be printed
-i 'TAG[*]!=1'  .. true
-e 'TAG[*]=1'   .. false, the record will be discarded
-e 'TAG[*]!=1'  .. false
-i 'TAG[0]=1'   .. true
-i 'TAG[0]!=1'  .. false
-e 'TAG[0]=1'   .. false
-e 'TAG[0]!=1'  .. true<br>
表达式示例:
MIN(DV)>5
MIN(DV/DP)>0.3
MIN(DP)>10 & MIN(DV)>3
FMT/DP>10  & FMT/GQ>10 .. both conditions must be satisfied within one sample
FMT/DP>10 && FMT/GQ>10 .. the conditions can be satisfied in different samples
QUAL>10 |  FMT/GQ>10   .. true for sites with QUAL>10 or a sample with GQ>10, but selects only samples with GQ>10
QUAL>10 || FMT/GQ>10   .. true for sites with QUAL>10 or a sample with GQ>10, plus selects all samples at such sites
TYPE="snp" && QUAL>=10 && (DP4[2]+DP4[3] > 2)
COUNT(GT="hom")=0      .. no homozygous genotypes at the site
AVG(GQ)>50             .. average (arithmetic mean) of genotype qualities bigger than 50
ID=@file       .. selects lines with ID present in the file
ID!=@~/file    .. skip lines with ID present in the ~/file
MAF[0]<0.05    .. select rare variants at 5% cutoff
POS>=100   .. restrict your range query, e.g. 20:100-200 to strictly sites with POS in that range.

表达式在Shell中的应用:
需要注意的是表达式必须经常用引号括起来,因为一些特殊字符在shell中具有特殊的含义。以下是一个正确示例
bcftools view -i '%ID!="." & MAF[0]<0.01'
©著作权归作者所有,转载或内容合作请联系作者
禁止转载,如需转载请通过简信或评论联系作者。
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 151,829评论 1 331
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 64,603评论 1 273
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 101,846评论 0 226
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 42,600评论 0 191
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 50,780评论 3 272
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 39,695评论 1 192
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 31,136评论 2 293
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 29,862评论 0 182
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 33,453评论 0 229
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 29,942评论 2 233
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 31,347评论 1 242
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 27,790评论 2 236
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 32,293评论 3 221
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 25,839评论 0 8
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 26,448评论 0 181
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 34,564评论 2 249
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 34,623评论 2 249

推荐阅读更多精彩内容