SV工具包SURVIVOR的使用
in 生信工具 with 0 comment

SV工具包SURVIVOR的使用

in 生信工具 with 0 comment

SURVIVOR 是一套用于模拟/评估 SV、合并和比较样本内及样本间 SV 的工具集,其还可对SV进行重新格式化及进行统计总结,其他详细信息也可参考工具的wiki

安装

git clone https://github.com/fritzsedlazeck/SURVIVOR.git
cd SURVIVOR/Debug
make

快速使用

Program: SURVIVOR (Tools for Structural Variations in the VCF format)
Version: 1.0.7

Usage: SURVIVOR <command> [options]

Commands:
-- Simulation/ Evaluation
	simSV	    根据一个参考基因组模拟SV和SNP.
	scanreads	从 mapped reads 中获取错误信息用于模拟.
	simreads	模拟长 reads (Pacio or ONT).
	eval	    评估从模拟数据中鉴定的SV.

-- Comparison/filtering
	merge	比较或合并VCF以生成一致或多个样本的VCF.
	filter	对VCF进行特定SV大小或区域的过滤
	stats	统计VCF中的多个信息
	compMUMMer	使用MUMMer (Show-diff)找到的断点来注释VCF.

-- Conversion
	bincov	            将bins coverage vector 转为bed文件来过滤低MQ区域的SV
	vcftobed	        将VCF格式转为bed格式
	bedtovcf	        将bed格式转为VCF格式
	smaptovcf	        将smap文件转为VCF文件 (beta version)
	bedpetovcf	        将bedpe文件转为VCF文件(beta version)
	hapcuttovcf	        使用原始提供给Hapcut2的SNP文件将Hapcut2最终结果转为VCF文件 
	convertAssemblytics	将Assemblytics结果转换为VCF文件

常用功能

这里介绍下SURVIVOR工具包中的几个常用功能,其他功能及详细参数介绍参见Wiki的Methods & Parmater部分。

模拟和评估SV

关于SV模拟的详细参数介绍参见:Methods & Parmater
首先,需要生成参数配置文件:

SURVIVOR simSV parameter_file

注:可以更改参数配置文件中的具体参数,但不要更改每行的位置
然后,基于参考序列进行SV模拟,参考序列需为fasta格式:

SURVIVOR simSV ref.fa parameter_file 0.1 0 simulated

该步骤使用的read模拟工具为Mason,也可以选择其他模拟工具。
下一步,对模拟的SV进行评估:

SURVIVOR eval caller.vcf simulated.bed 10 eval_res

该步骤会使用之前的模拟数据集对caller.vcf进行评估,并允许断点上下游有10bp的差异,评估结果保存在eval_res中。
最后会输出一行评估的结果:

Overall: 20 11/0/0/0/9 0/0/0/0/0 0/0/0/0/0 1 0

这表明一共模拟了20个SV,然后的三组数据分别是真阳性结果(即模拟了并找到的),假阴性结果(模拟了但没找到),假阳性结果(未模拟但找到了)。每组中的五个数字分别代表五种SV类型(DEL/DUP/INV/TRA/INS)。倒数第二个表示灵敏度,最后一个数字表示FDR概率。总之,在这次模拟中,总共模拟了20个SV,从中再次发现11个缺失和9个插入,并且没有丢失或报告其他的任何SV。

多样本间SV的比较

为了降低一些SV鉴定工具的假阳性,同时保持较高的灵敏度,一般的建议是一个样本使用多个鉴定工具。例如,可以在短reads数据上运行Manta、GRIDSS和Lumpy。每个工具都会生成一个VCF文件作为结果。
这时SURVIVOR就得安排上了,它可比较这些VCF文件并生成一致性结果。注意,通常,SURVIVOR可以合并任何技术或其他工具生成的SV VCF文件
首先, 将VCF文件名输出到一个文件中:

ls *vcf > sample_files

然后, 生成合并的VCF结果:

SURVIVOR merge sample_files 1000 2 1 1 0 30 sample_merged.vcf

1000表示允许合并的SV间的距离最大为1000bp;
2表示仅输出2个工具均鉴定出的SV;
1表示仅输出2个工具鉴定出的同类型的SV;
1表示仅输出2个工具鉴定出的同方向的SV;
30表示仅考虑长度在30bp以上的SV

识别SV鉴定错误的区域

首先,从排好序的bam文件中提取比对质量较差的reads:

samtools view -H our.sort.bam > lowMQ.sam
samtools view our.sort.bam | awk '$5<5 {print $0}' >>  lowMQ.sam
samtools view -S -b -h lowMQ.sam > lowMQ.bam

现在就得到了MQ < 5 的reads的bam文件, 然后, 需要计算碱基覆盖度:

 samtools depth lowMQ.bam >  lowMQ.cov

下一步, 将覆盖度文件进行聚类生成bed文件, 从而用于SV的过滤:

SURVIVOR bincov lowMQ.cov 10 2 > lowMQ.bed

该步骤允许将距离最大10bp的区域聚在一起,并且仅在覆盖度大于2的情况下才考虑该区域。生成的bed文件可用于过滤SV时使用.

参考