pbsv Structural Variant Caller

The standard pbsv workflow for calling and analyzing structural variants (SVs) in diploid genomes from PacBio SMRT reads consists of three main stages:

1. Alignment: Align PacBio reads (Subreads or CCS) to a reference genome using an aligner like pbmm2.
2. Discovery: Identify structural variation signatures from the aligned BAM files. This step reduces aligned reads to only those relevant to SVs, producing .svsig.gz files.
3. Calling: Call structural variants and assign genotypes jointly across all samples using the discovered signatures to produce a .vcf file.

Call structural variants jointly for all samples using their respective .svsig.gz files.

Important: If your input signatures were generated from CCS reads, you must include the --ccs flag.

pbsv call ref.fa ref.sample1.svsig.gz ref.sample2.svsig.gz ref.var.vcf

For large genomes or high coverage, process chromosomes separately to improve performance.

1. Generate separate .svsig.gz files for each chromosome using the --region flag in pbsv discover.
2. Call SVs using the -j flag to specify the number of threads.

# 1. Generate per-chromosome signatures
for i in $(samtools view -H hg38.movie1.bam | grep '^@SQ' | cut -f2 | cut -d':' -f2); do
    pbsv discover --region $i hg38.movie1.bam hg38.sample1.$i.svsig.gz
done

# 2. Call SVs using 8 threads
pbsv call -j 8 hg38.fa hg38.sample1.*.svsig.gz hg38.sample1.vcf

The pbsv suite can be installed using the Bioconda package manager.

For the recommended alignment step, you should also install pbmm2 via conda.

conda install pbsv
conda install pbmm2

Before using pbsv, you must align your reads to a reference genome. The sample name (stored in the SM tag of the read groups) is required for downstream joint calling.

Depending on your input type, use the following pbmm2 commands:

# Subreads BAM input
pbmm2 align ref.fa movie1.subreads.bam ref.movie1.bam --sort --median-filter --sample sample1

# CCS BAM input
pbmm2 align ref.fa movie1.ccs.bam ref.movie1.bam --sort --preset CCS --sample sample1

# CCS FASTQ input
pbmm2 align ref.fa movie1.Q20.fastq ref.movie1.bam --sort --preset CCS --sample sample1 --rg '@\tID:movie1'

You can find example datasets for Genome in a Bottle sample HG002 to test pbsv at the following locations:

• 10-fold CLR coverage: https://downloads.pacbcloud.com/public/dataset/HG002/Sequel-201810/
• 28-fold CCS coverage: https://bit.ly/2RW1b3I

Use pbsv discover to identify signatures of structural variation from aligned BAM files. This produces .svsig.gz files.

Recommendations:

• Provide a tandem repeat annotation .bed file via the --tandem-repeats flag to increase sensitivity and recall.
• Optionally index the resulting .svsig.gz files using tabix to allow random access during the calling stage via pbsv call -r.

Sample names are automatically transferred from the RG headers in the BAM to the .svsig.gz file.

pbsv discover ref.movie1.bam ref.sample1.svsig.gz

# Optional: index for random access
tabix -c '#' -s 3 -b 4 -e 4 ref.sample1.svsig.gz

Use the findTandemRepeats --merge command to identify and merge tandem repeat regions. This command requires a reference genome file in FASTA format and a TRF (Tandem Repeat Finder) output file in BED format.

findTandemRepeats --merge <REF>.fa <REF>.trf.bed

Thresholds used to determine if a variant is called and assigned a genotype.

Genotyping thresholds:
  -A,--call-min-reads-all-samples [3]
  -O,--call-min-reads-one-samples [3]
  -S,--call-min-read-per-one-sample [1]
  -B,--call-min-bnd-reads-all-samples [2]
  -P,--call-min-reads-per-strand-all-samples [20]

Note: For CCS input with `--ccs`, thresholds are relaxed to: -A 3 -B 2 -O 3 -S 0 -P 10.