基于paired-end mapping methods 和 split-read analysis,检测结构变异。
写在前面
Delly is an integrated structural variant (SV) prediction method that can discover, genotype and visualize deletions, tandem duplications, inversions and translocations at single-nucleotide resolution in short-read massively parallel sequencing data. It uses paired-ends, split-reads and read-depth to sensitively and accurately delineate genomic rearrangements throughout the genome. Structural variants can be visualized using Delly-maze and Delly-suave.
参考资料
Manual: https://github.com/dellytools/delly
FAQ: https://groups.google.com/forum/#!topic/delly-users/fOL8J3Exod4
实战
版本:Delly v0.7.8 可以一次性call 所有类型的SVs;低版本的需要通过-t
指定SV类型
00 输入文件
- Bam files need to be sorted, indexed and ideally duplicate marked. If multiple libraries are present for a single sample these need to be merged in a single bam file with unique ReadGroup tags.
- .excl: a file include information of telomere and centromere regions and also all unplaced contigs, those regions will be removed from calling.
- 有群体时,建议将所有个体的
.bcf
文件合并后,再进行genotyping。
#install
git clone --recursive https://github.com/dellytools/delly.git
cd delly/
make all
#也可以直接下载编译后的二进制文件,直接使用,不需安装 https://github.com/dellytools/delly/releases/
#Running
DB="/root/cc/DB"
i="SXY"
samtools index $i.recal.bam
~/cc/biosoft/delly/delly_v0.7.8_linux_x86_64bit call -g $DB/hg19_chr.fa -x exc1.excl -o delly_out/$i.exc1.bcf $i.recal.bam
#bcftools view $i.target.bcf >$i.target.vcf
#Merge SV sites into a unified site list
delly merge -o sites.bcf s1.bcf s2.bcf ... sN.bcf
#Genotype this merged SV site list across all samples. This can be run in parallel for each sample.
delly call -g $DB/target-ref.fa -v sites.bcf -o s1.geno.bcf s1.bam
#Merge all genotyped samples to get a single VCF/BCF using bcftools merge
bcftools merge -m id -O b -o merged.bcf s1.geno.bcf s2.geno.bcf ... sN.geno.bcf
#Apply the germline SV filter which requires at least 20 unrelated samples
delly filter -f germline -o germline.bcf merged.bcf
bcftools view germline.bcf >germline.vcf
01 SV统计:长度分布、在各样本中的情况等
长度统计:DEL,DUP和INV,可以用INFO:END-POS获得;INS的长度值为 INFO:INSLEN,而[POS, INFO:END] 提供的大概的插入位置信息。
这里使用作者开发的一个小工具,统计size分布、每个样本的变异数目等。
git clone --recursive https://github.com/dellytools/svprops.git
cd svprops/
./svprops germline.vcf >germline.tab #size分布,各类型
./sampleprops germline.vcf >sample.tab #每个样本在多少SV中是RR,RA,AA(基因型统计)
#绘图
Rscript R/svprops.R germline.tab
Rscript R/sampleprops.R sample.tab
head -3germline.tab
#chr start chr2 end id size vac vaf singleton missingrate svtype ct precise ci inslen homlen ce refgq altgq gqsum rdratio medianrc refratio altratio maxaltratio PEsupport SRsupport supportsum
#chr22 17900897 chr22 17900914 DEL00000000 18 7 0.152174 NA 0 DEL 3to5 1 6 0 5 1.93828 75 10000 41614 0.563525 632 0 0.721311 0.787879 0 228 737
#chr22 18003439 chr22 18004197 DEL00000006 759 9 0.195652 NA 0 DEL 3to5 0 174 0 0 0 15 45 558 0.658926 107 0 0.25 0.333333 11 0 123
head -3 sample.tab
#sample missing homref het homalt
#sample1 0 76 37 14
#sample2 0 69 47 11
germline.tab各列信息如下
- vac: variant allele count (across all samples)
- vaf: variant allele frequency (across all samples)
- singleton: N/A if present in multiple samples or sample name if only present in one
- missingrate: how many samples have a missing genotype (useful after GQ calibration)
- ct: Delly's INFO:CT,Paired-end signature induced connection type
- precise: 1 if INFO:PRECISE,Precise structural variation
- ci: Delly's INFO:CIPOS,PE confidence interval around POS
- inslen: insertion length
- homlen: homology length
- ce: Delly's INFO:CE,Consensus sequence entropy
- refgq: median reference genotype quality for all SV non-carriers (GT==0/0)
- altgq: median alternative genotype quality for het. SV carriers
- gqsum: total GQ sum (useful to flag repetitive sites that are poorly genotyped)
- rdratio: read-depth ratio of SV carriers to non-carriers (useful for filtering CNVs)
- medianrc: median coverage
- refratio: median REF support for non-carriers
- altratio: median ALT support for carriers
- maxaltratio: max. ALT support for carriers
- PEsupport: total paired end support across all samples
- SRsupport: total split-read support across all samples
- supportsum: total depth across all samples
sample.tab各列信息如下
- missing:该个体有多少位点的genotype为missing
- homref:该个体有多少位点为ref的纯合位点
- het:该个体有多少位点为杂合位点
-
homalt:该个体有多少位点为alt的纯合位点
Note
样本数比较少时,利用
filter
可能使所有SV被过滤掉,作者建议根据INFO FIELD(PASS)来过滤SV。exc1.excl:排除区域(端粒等区域,DELLY github上有参考文件)
chr1 #排除整条染色体
chr22 0 10000 telomere
chr22 13000000 16000000 centromere
chr22 51294566 51304566 telomere
Error
bcf文件转为vcf时,发现里面一片空白,可能是bcftools版本过低,建议更新后再转换格式即可。
bcftools view file.bcf >file.vcf