组装生物信息群体遗传学

生信分析|Minimap2+sniffles calling S

2021-08-24  本文已影响0人  行走的脱发机double

1、minimap2

Minimap2 is a fast sequence mapping and alignment program that can find overlaps between long noisy reads, or map long reads or their assemblies to a reference genome optionally with detailed alignment (i.e. CIGAR). At present, it works efficiently with query sequences from a few kilobases to ~100 megabases in length at a error rate ~15%. Minimap2 outputs in the PAF or the SAM format.
pacbio reads mapping reference genome

minimap2 -a ref.fa pacbio_reads.fasta --MD | samtools view -@ 10 -bS | samtools sort -@ 10 -o pacbio_reads_mapping_ref.sorted.bam

minimap2比对必选参数

-a  Generate CIGAR and output alignments in the SAM format. Minimap2 outputs in PAF by default.
--MD    Output the MD tag (see the SAM spec).(must for sniffles)
-t  Number of threads [3].

2、sniffles

Sniffles is a structural variation caller using third generation sequencing (PacBio or Oxford Nanopore). It detects all types of SVs (10bp+) using evidence from split-read alignments, high-mismatch regions, and coverage analysis. Please note the current version of Sniffles requires sorted output from BWA-MEM (use -M and -x parameter), Minimap2 (sam file with Cigar & MD string) or NGMLR.

sniffles -m pacbio_reads_mapping_ref.sorted.bam -s 4 -v Q1reads_mapping_Q1ref.vcf

sniffles call SVs必选参数

-m  bam file
-v out vcf file
-s <int>,  --min_support <int> Minimum number of reads that support a SV. [10]

SV calling for a population

obtain a multi sample vcf file from long reads using Sniffles and SURVIVOR.
Download SURVIVOR:

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

Merge all the vcf files across all samples

lack the information if a SV identified in one sample but not in the other is really absent.

#1.Please sort the files prior to merging.
#2.put all file paths over all vcf files together e.g.:
ls *sample.vcf > vcf_files_raw_calls.txt
#3.call SURVIVOR to merge these into one vcf file:
./SURVIVOR merge vcf_files_raw_calls.txt 1000 1 1 -1 -1 -1 merged_SURVIVOR_1kbpdist_typesave.vcf

Force call all the SVs across all the samples

obtain all calls even if they might be 0/0 in all samples.

./sniffles -m my_sample.bam -v my_sample_gt.vcf --Ivcf merged_SURVIVOR_1kbpdist_typesave.vcf
ls *sample_gt.vcf > vcf_files_gt_calls.txt
./SURVIVOR merge vcf_files_gt_calls.txt 1000 -1 1 -1 -1 -1 merged_gt_SURVIVOR_1kbpdist_typesave.vcf
#The -1 for the minimum SV caller is necessary

SURVIVOR merge参数

vcf_files_gt_calls.txt    File with VCF names and paths
1000    max distance between breakpoints (0-1 percent of length, 1- number of bp) 
-1    Minimum number of supporting caller
1    Take the type into account (1==yes, else no)
-1   Take the strands of SVs into account (1==yes, else no)
-1    Disabled.
-1    Minimum size of SVs to be taken into account.
merged_gt_SURVIVOR_1kbpdist_typesave.vcf    Output VCF filename

1000表示允许合并的SV间的距离最大为1000bp;a maximum allowed distance of 1kb, as measured pairwise between breakpoints (begin1 vs begin2, end1 vs end2).
-1表示强制输出所有SV;
1表示仅输出2个工具鉴定出的同类型的SV;
-1表示输出2个工具鉴定出的所有方向的SV;
-1长度表示不过滤长度较小的SV

参考:

minimap2帮助文档:https://lh3.github.io/minimap2/minimap2.html
https://github.com/fritzsedlazeck/Sniffles
https://github.com/lh3/minimap2
SURVIVOR:https://www.jianshu.com/p/01110d85e71f

上一篇下一篇

猜你喜欢

热点阅读