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：http://08643.cn/p/01110d85e71f
