Hprc banner tamu.png

Difference between revisions of "Bioinformatics:Sequence Variants"

Jump to: navigation, search
(Sequence Variants (SNPs, indels))
(Call Variants)
Line 4: Line 4:
== Call Variants ==
== Call Variants ==

Latest revision as of 09:40, 13 October 2021

Sequence Variants (SNPs, indels)

Back to Bioinformatics Main Menu

Call Variants


GCATemplates available: grace

 module spider FreeBayes

FreeBayes is a Bayesian genetic variant detector designed to find small polymorphisms, specifically SNPs (single-nucleotide polymorphisms), indels (insertions and deletions), MNPs (multi-nucleotide polymorphisms), and complex events (composite insertion and substitution events) smaller than the length of a short-read sequencing alignment.


GCATemplates available: grace

 module spider Platypus

Platypus homepage github page

Platypus is a tool designed for efficient and accurate variant-detection in high-throughput sequencing data. By using local realignment of reads and local assembly it achieves both high sensitivity and high specificity. Platypus can detect SNPs, MNPs, short indels, replacements and (using the assembly option) deletions up to several kb.


GCATemplates available: no

GATK 4.0

Use the following command to run GATK version 4.0+ instead of using the .jar file


If whole genome processing for a large genome is taking a long time, you can use intervals to parallelize the analysis and analyse sequential regions of a chromosome or more specifically just coding regions.

Below is a step-by-step breakdown of the Best Practices workflow, with a detailed explanation of why -L should or shouldn’t be used with each tool.

Tool -L Why or why not?
BaseRecalibrator YES This excludes off-target sequences and sequences that may be poorly mapped, which have a higher error rate. Including them could lead to a skewed model and bad recalibration.
PrintReads NO Output is a bam file; using -L would lead to lost data.
UnifiedGenotyper/Haplotype Caller YES We’re only interested in making calls in exome regions; the rest is a waste of time & includes lots of false positives.
Next steps NO No need since subsequent steps operate on the callset, which was restricted to the exome at the calling step.

GATK 3.0

Recommended workflows for variant discovery analysis with GATK

 module spider GATK

To execute GATK run:

 java -jar $EBROOTGATK/GenomeAnalysisTK.jar

The GATK is the industry standard for identifying SNPs and indels in germline DNA and RNAseq data

The GATK tools are primarily designed to process exomes and whole genomes generated with Illumina sequencing technology, but they can be adapted to handle a variety of other technologies and experimental designs.


GCATemplates available: yes (look at bwa templates)

 module spider SAMtools

SAM Tools provide various utilities for manipulating alignments in the SAM format, including sorting, merging, indexing and generating alignments in a per-position format.

Example of sorting and creating a bam file without creating an intermediate sam file (reserving the entire compute node 20 cores).

Notice that using the $TMPDIR for temporary files is useful since files in $TMPDIR do now count against your quota.

bwa mem -M -t 20 -R "@RG\tID:my_readgroup\tLB:my_library\tSM:my_sample\tPL:ILLUMINA" \
 my_ref_genome.fasta paired_end_R1.fastq.gz paired_end_R2.fastq.gz | samtools view -h -Sb - | \
 samtools sort -m 2G --threads 10 - -T $TMPDIR/sorted -o my_output_bam


GCATemplates available: terra

IMSindel An accurate intermediate-size indel detection tool incorporating de novo assembly and gapped global-local alignment with split read analysis.

module load IMSindel/1.0.1-foss-2019b-Ruby-2.7.1

The output directory option --outd is required and the directory must exist prior to running imsindel.

The memory usage is low and you can complete the analysis in a timely manner by using half the cores on a compute node since the early step uses multiple cores but the majority of the remaining processing uses one core.


GCATemplates available: no

 module load Parsnp/1.2-Linux64

Parsnp is a command-line-tool for efficient microbial core genome alignment and SNP detection. Parsnp was designed to work in tandem with Gingr, a flexible platform for visualizing genome alignments and phylogenetic trees; both Parsnp and Gingr form part of the Harvest suite.


GCATemplates available: no

 module load kSNP/3.021-Linux_package

kSNP identifies the pan-genome SNPs in a set of genome sequences, and estimates phylogenetic trees based upon those SNPs. SNP discovery is based on k-mer analysis, and requires no multiple sequence alignment or the selection of a reference genome, so kSNP can take 100's of microbial genomes as input.


GCATemplates available: ada

Delly homepage

Delly is an integrated structural variant (SV) prediction method that can discover, genotype and visualize deletions, tandem duplications, inversions and translocations at single-nucleotide resolution in short-read massively parallel sequencing data. It uses paired-ends, split-reads and read-depth to sensitively and accurately delineate genomic rearrangements throughout the genome. Structural variants can be visualized using Delly-maze and Delly-suave.

module spider Delly

Delly primarily parallelizes on the sample level. Hence, OMP_NUM_THREADS should be always smaller or equal to the number of input samples.

If you have 5 samples, for example, you would use the following in your job script along with the appropriate #SBATCH values:


Annotate Variants


GCATemplates available: terra

snpEff homepage

 module spider snpEff

Genetic variant annotation and effect prediction toolbox.

To see a list of available snpEff genome databases

 java -jar $EBROOTSNPEFF/snpEff.jar databases

Download pre-built database

This must be done on a login node and not in a job script.

1) First load the snpEff version of your choice

module load snpEff/4.3t-foss-2018b-Python-3.6.6-Java-1.8.0

2) Copy the main snpEff.config file to your current directory

cp $EBROOTSNPEFF/snpEff.config ./

3) See if the genome is available (save output to a file then search the file):

java -Xmx4g -jar $EBROOTSNPEFF/snpEff.jar databases > avail_dbs.txt

4) Download database using your snpEff.config file which will download it to the current directory in a directory called data:

java -Xmx4g -jar $EBROOTSNPEFF/snpEff.jar download -c snpEff.config Sporisorium_reilianum

5) Use the -config option when running snpEff commands in your job script:

java -Xmx4g -jar $EBROOTSNPEFF/snpEff.jar -config snpEff.config Sporisorium_reilianum my_file.vcf > annotated.vcf

Create your own database

You can create your own custom database with a sequences.fasta and genome.gff3 file (or .gff).

This must be done on a login node and not in a job script.

The following example creates a database for AaegL5.0

1) First load the snpEff version of your choice

module load snpEff/4.3t-foss-2018b-Python-3.6.6-Java-1.8.0

2) Make a directory for your new config file and database files

mkdir -p $SCRATCH/snpeff/4.3t

3) Copy the main snpEff.config file to your new directory

cp $EBROOTSNPEFF/snpEff.config $SCRATCH/snpeff/4.3t/

4) Edit your snpEff.config file and change the data.dir to the directory in your $SCRATCH (use your NETID in the full path)

data.dir = /scratch/user/NETID/snpeff/4.3t

5) Add a unique line to your snpEff.config file with your custom genome name using spaces:

AaegL5.0.genome : AaegL5.0

6) Make a directory with the same name that you added to the snpEff.config file and copy your unzipped reference genome and gff3 files to the new directory $SCRATCH/snpeff/4.3t/

mkdir $SCRATCH/snpeff/4.3t/AaegL5.0
cp my_reference_genome.fa $SCRATCH/snpeff/4.3t/AaegL5.0/
cp my_reference_genome.gff $SCRATCH/snpeff/4.3t/AaegL5.0/

7) Create symlinks in the new directory to your two files with names sequence.fa and genes.gff

cd $SCRATCH/snpeff/4.3t/AaegL5.0/
ln -s my_reference_genome.fa sequences.fa
ln -s my_reference_genome.gff genes.gff

8) Make the database using the -config option to specify the new config file.

java -jar $EBROOTSNPEFF/snpEff.jar build -config $SCRATCH/snpeff/4.3t/snpEff.config -gff3 AaegL5.0

This will create a snpEffectPredictor.bin file in the build.dir

9) You can delete the reference and gff files from the directory after you create the database

cd $SCRATCH/snpeff/4.3t/AaegL5.0/
unlink sequences.fa
unlink genes.gff
rm my_reference_genome.fa
rm my_reference_genome.gff

10) Use the -config option when running snpEff commands in your job script:

java -jar $EBROOTSNPEFF/snpEff.jar -config $SCRATCH/snpeff/4.3t/snpEff.config AaegL5.0  my_file.vcf > annotated.vcf

Variant Effect on Protein


GCATemplates available: ada

PROVEAN homepage

 module spider PROVEAN

PROVEAN (Protein Variation Effect Analyzer) is a software tool which predicts whether an amino acid substitution or indel has an impact on the biological function of a protein.

PROVEAN is useful for filtering sequence variants to identify nonsynonymous or indel variants that are predicted to be functionally important.

Since it uses BLAST and psiBLAST against the nr database, it is slow and may take several minutes per query.

If you are working with human variants, you can download a file of precomputed scores for every amino acid change at every position for all human genes at: http://provean.jcvi.org/downloads.php