Skip to content
Detect complex indels from next-generation sequencing reads
Branch: master
Clone or download
Fetching latest commit…
Cannot retrieve the latest commit at this time.
Permalink
Type Name Latest commit message Commit time
Failed to load latest commit information.
LICENSE
README.md
indelseek.pl

README.md

INDELseek

Detect complex indels from next-generation sequencing reads

Download latest version in a ZIP package

INDELseek algorithm (Au et al. BMC Genomics 2017 18:16 Figure 1) INDELseek algorithm

Dependencies, as tested on 64-bit CentOS 5.5

Usage

INDELseek.pl analyzes SAM alignments in plain text to call complex indels in VCF format. Users are recommeneded to focus on a single BAM file and a small region at a time (e.g. a single PCR amplicon or single exon of protein-coding gene). Please refer to the following example to scale up the analysis to multiple BAM files and genomic regions in parallel.

samtools view bam region | ./indelseek.pl [options]

  • bam: indexed BAM alignment file
  • region: 1-based position coordinates, e.g. amplicon chr4:55589744-55589911. For details please refer to documentation of samtools view.
  • standard output (STDOUT): variant calls in VCF version 4.1

Options

  • --refseq FILE: indexed reference genome in FASTA format [ucsc.hg19.fasta]
  • --samtools FILE: path to samtools executable [samtools]
  • --max_distance INT: maximum distance of two variants to be considered together as a complex indel [5]
  • --quality_threshold INT: complex indels with mean read base quality score below INT are marked LowQual [20]
  • --skip_lowqual: skip LowQual complex indel calls in output
  • --min_depth INT: complex indels with depth (forward and reverse reads with the complex indel) below INT are marked LowDepth [50]
  • --skip_lowdepth: skip LowDepth complex indel calls in output
  • --min_af FLOAT: complex indels with allele frequendcy below FLOAT are marked LowAF [0]
  • --skip_lowaf: skip LowAF complex indel calls in output
  • --depth_bam FILE: indexed BAM alignment file is needed for allele frequency calculation
  • --phredoffset INT: read base quality score encoding offset. Example: 33 for Illumina 1.8+ or Sanger, 64 for Illumina 1.3+ / 1.5+ or Solexa. (https://en.wikipedia.org/wiki/FASTQ_format#Encoding) [33]

Examples

# Any complex indel in KIT exon 8 with minimum allele frequency 2%?
# (Actual output of sample 10 described in manuscript)
samtools view sample.bam chr4:55589744-55589911 | ./indelseek.pl --skip_lowqual --skip_lowdepth --skip_lowaf --min_af 0.02 --depth_bam sample.bam | tee sample.complexindel.vcf
##fileformat=VCFv4.1
##source=INDELseek
##reference=file://ucsc.hg19.fasta
##INFO=<ID=DP2,Number=2,Type=Integer,Description="# alt-foward and alt-reverse reads">
##INFO=<ID=QS2,Number=2,Type=Float,Description="Mean quality scores of alt-foward and alt-reverse bases">
##INFO=<ID=AF,Number=1,Type=Float,Description="Allele Frequency">
##INFO=<ID=RDP,Number=1,Type=Float,Description="Mean read depth of REF positions">
##FILTER=<ID=LowAF,Description="AF below 0.02">
##FILTER=<ID=LowDepth,Description="ALT depth below 50">
##FILTER=<ID=LowQual,Description="Mean quality scores below 20 or ALT contains N">
#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO
chr4	55589768	.	CTTACGACA	AACCTC	35.63	PASS	DP2=857,752;QS2=37.44,33.56;AF=0.119;RDP=13525.3
chr4	55589767	.	ACTTACGACA	GGATGGAACT	36.33	PASS	DP2=251,203;QS2=37.37,35.04;AF=0.033;RDP=13651.1
chr4	55589766	.	GACTTACGA	TTTCCG	35.76	PASS	DP2=208,184;QS2=37.46,33.83;AF=0.029;RDP=13724.2
chr4	55589769	.	TTACGACA	CTCCT	35.68	PASS	DP2=151,125;QS2=37.61,33.35;AF=0.021;RDP=13376.0

# Speed up complex indel detection in list of BAM files and genomic regions by GNU Parallel
parallel --tag "samtools view {1} {2} | indelseek.pl --skip_lowqual --skip_lowdepth --skip_lowaf --min_af 0.02 --depth_bam {1}" ::: *.bam :::: region.list | grep -v "#" > complexindel.tab

Citations

INDELseek algorithm

Example studies that used INDELseek

You can’t perform that action at this time.