Long Read SV calling and JBrowse viewing - mestato/EPP622_2026 GitHub Wiki

We need some new software

conda activate svs
conda install bcftools survivor

Set up.

mkdir 03_cutesv
cd 03_cutesv

slurm script. using recommended hifi params from cutesv documentation, also not specifying memory, cutesv doesn't need much ram.

#!/bin/bash
#SBATCH --job-name cutesv
#SBATCH --nodes=1
#SBATCH --cpus-per-task=8
#SBATCH -A ISAAC-UTK0505
#SBATCH --partition=campus
#SBATCH --qos=campus
#SBATCH --time=12:00:00
#SBATCH --mail-type=BEGIN,END,FAIL
#SBATCH [email protected]

# Activate the conda environment
eval "$(conda shell.bash hook)"
conda activate svs

# Run cutesv on two sequence sets
for BAM in /lustre/isaac24/proj/UTK0505/youruserid/module3/02_minimap2/*bam
do
 # The */ matches any text ending with a forward slash.
  SAMPLE="${BAM##*/}"
  # second, remove .bam from the file name
  SAMPLE="${SAMPLE%.bam}"

  # We need a working directory, separate for each sample, for cutesv
  mkdir -p ${SAMPLE}_workdir

  cuteSV \
  $BAM \
  /lustre/isaac24/proj/UTK0505/module3_data/AppSpring_1.0.0_Hap1.Chr.fasta \
  $SAMPLE.vcf \
  ${SAMPLE}_workdir \
  --max_cluster_bias_INS 1000 \
  --diff_ratio_merging_INS 0.9 \
  --max_cluster_bias_DEL 1000 \
  --diff_ratio_merging_DEL 0.5 \
  --genotype \
  --threads 8

done

Lets get stats from SURVIVOR

SURVIVOR stats AppSpring2.vcf -1 -1 -1 AppSpring2.vcf.stats.txt
SURVIVOR stats CherBrave3.vcf -1 -1 -1 CherBrave3.vcf.stats.txt

Lets also check out our 0/0, 0/1, and 1/1 rates.

grep -c '0/0' *vcf 
grep -c '0/1' *vcf
grep -c '1/1' *vcf

What is suspicious about those 0/0 calls? How many are IMPRECISE? grep -c 'IMPRECISE' *vcf

Lets look at a well behaved call from Cher Brave that is 1/1.

grep '1/1' CherBrave3.vcf | head -n 1000 | tail -n 10

Why might I not want the first ten?

Lets check out this one. Line from the vcf:

AppSpring_1.0.0_Hap1_Chr09	11910377	cuteSV.DEL.1948	TATCTATGCAAGGACAAGGCATAAAATCTAT	T	610.7	PASS	PRECISE;SVTYPE=DEL;SVLEN=-30;END=11910407;CIPOS=-515,515;CILEN=-5,5;RE=126;AF=0.8182;STRAND=+-	GT:DR:DV:PL:GQ	1/1:28:126:611,84,0:83
AppSpring_1.0.0_Hap1_Chr09	11913481	cuteSV.INS.1708	A	ACACGAACTGAGCTAAAACCTCATACGAATTGGCACGCTAACTCAATGTCAACTG	648.9	PASS	PRECISE;SVTYPE=INS;SVLEN=54;END=11913481;CIPOS=-16,16;CILEN=-1,1;RE=111;AF=0.8346	GT:DR:DV:PL:GQ	1/1:22:111:649,103,0:102

Lets decode the genotype stats

GT = genotype
DR = high quality reference reads
DV = high quality variant reads
PL = Phred-scaled genotype likelihoods rounded to the closest integer
GQ = Genotype quality	

I'm going to downsample the bam to select the alignments right around this variant - 100kb up and down stream. This makes the bam much smaller.

cd ../02_minimap2
samtools view  -b -h -o CherBrave3.Chr9_12Mb.bam CherBrave3.bam AppSpring_1.0.0_Hap1_Chr09:11810377-12010377
samtools index CherBrave3.Chr9_12Mb.bam
ls -lah Cher*bam

Our bam went from 26Gb to 47Mb.

On a laptop terminal, transfer files to your laptop. Bam, bam.bai, reference genome fasta, and ref genome fai.

scp isaac:/lustre/isaac24/proj/UTK0505/mstaton1/module3/02_minimap2/CherBrave3.Chr9_12Mb.bam\* .
scp isaac:/lustre/isaac24/proj/UTK0505/module3_data/AppSpring_Hap1_Chr9.fasta\* .
scp isaac:/lustre/isaac24/proj/UTK0505/mstaton1/module3/03_cutesv/CherBrave3.vcf .

Load into JBrowse and lets take a look.