How do I generate a variant list (i.e. VCF file) using Illumina reads from a human genome? This is a problem I have to solve frequently, and I'd be interested in knowing what other methods people use to solve the same problem.

About twice a year, I get asked to determine variants from Illumina reads, usually from either mouse or human. These are things that have good reference genomes (e.g. human genome from ensembl, mouse genome from ensembl), and are frequently tested out on various bioinformatics tools.

I have at my disposal three computers, in order of usage preference: a laptop with 4 processing threads and 8GB of memory, a desktop with 12 processing threads and 64GB of memory, and a server with 24 processing threads and 256GB of memory. Occasionally I get access to better computers when I need to do lots of processing, or to do things quickly, but these comprise my basic bioinformatics team.

I've had whole-genome projects that involved processing data from 4 individuals, up to projects with about 120 individuals, with each individual having about 10-40M paired-end reads (100-125bp each). Beyond this level, I suspect other researchers will typically have their own bioinformatics team and/or capabilities, which is why I don't get any of the larger-scale jobs that I hear about from salaried bioinformaticians.

What is the standard approach that you would use to process these reads into reference-mapped variant data (i.e. gzipped VCF) for use in downstream analysis tools (e.g. VEP)?
1. Adapter Trimming

One of the first things I do after encountering a set of reads is to remove the adapter sequences from the start and end of reads. Most basecalling software includes some amount of built-in adapter trimming, but it is almost always the case that some adapter sequence will remain. Removing adapters is helpful for mapping because it reduces the effort the mapper needs to go through to make a match (which may help some borderline reads to be mapped properly).

My preferred adapter trimmer is Trimmomatic. It has the ability to search for palindromic adapter read through, and includes an adaptive sliding window option. Trimmomatic is threaded, and seems to work reasonably fast when run through files one at a time in threaded mode. Here's an example script that I might run to trim reads from some samples. This puts trimmed FASTQ files into a trimmed subdirectory:

TRIMPATH=/data/all/programs/trimmomatic/Trimmomatic-0.36;
JARPATH="${TRIMPATH}/trimmomatic-0.36.jar";
ADAPTfile="${TRIMPATH}/adapters/TruSeq3-PE-2.fa";
CMDS=" ILLUMINACLIP:${ADAPTfile}:2:30:10:5:true LEADING:3 TRAILING:3  SLIDINGWINDOW:10:20 MINLEN:40 ";

mkdir -p trimmed;

for INFILE1 in *_1.fastq.gz; do
  base=$(basename "$INFILE1" _1.fastq.gz);
  echo ${base};
  INFILE2="${base}_2.fastq.gz";
  OUTFILE_P1="trimmed/${base}_P1.fastq.gz";
  OUTFILE_P2="trimmed/${base}_P2.fastq.gz";
  OUTFILE_U1="trimmed/${base}_U1.fastq.gz";
  OUTFILE_U2="trimmed/${base}_U2.fastq.gz";
  java -jar "${JARPATH}" PE -threads 20 -phred33 \
        "${INFILE1}"  "${INFILE2}" \
        "${OUTFILE_P1}" "${OUTFILE_U1}" "${OUTFILE_P2}"  "${OUTFILE_U2}" \
        "${CMDS}";
done;


Trimmomatic is somewhat sensitive to bad input data, and will loudly fail (i.e. stop running) if it sees quality strings of different length to sequence strings:

Exception in thread "Thread-1" java.lang.RuntimeException: Sequence and quality length don't match: 'TACATGGCCCTGAAATGACTTTCACCCAGGCAACCAGTGCCCCCTGTATAGACACATGCCTTGGGCGCTCCCCACCCTTCCTCGCGTGGCCACACCTCTGT' vs '-AAFFJFJJ7A-FJJJFJFJFJJFFA-FFAF-&lt;A-&lt;FFFJA-&lt;-A7-F&lt;&lt;FFJJAJJJJJJJJJ--&lt;--7-7AJ7&lt;AAJA--J7&lt;ACTGCTGTGGGGCACCCAGCCCCCCAGATAGCCTGGCAGAAGGATGGGGGCACAGACTTCCCAGCTGCACGGGAGAGAC'
        at org.usadellab.trimmomatic.fastq.FastqRecord.&lt;init&gt;(FastqRecord.java:25)
        at org.usadellab.trimmomatic.fastq.FastqParser.parseOne(FastqParser.java:89)
        at org.usadellab.trimmomatic.fastq.FastqParser.next(FastqParser.java:179)
        at org.usadellab.trimmomatic.threading.ParserWorker.run(ParserWorker.java:42)
        at java.lang.Thread.run(Thread.java:745)


2. Read Mapping

Read mapping finds the most likely location for a read within a target genome.
There are extremely fast approximate mappers available, but these don't yet work for variant calling, and an exact mapping approach is necessary. My current preferred mapper is HISAT2. I use this instead of BWA due to the double-read issue, and because of the local variant-aware mapping. HISAT2 is made by the same computing group as Bowtie2 and Tophat2 (JHUCCB), and is their recommended tool for replacing those other programs.


  HISAT2 is a successor to both HISAT and TopHat2. We recommend that HISAT and TopHat2 users switch to HISAT2.


The HISAT2 page includes a link to a genomic index file for the human genome, including SNP variants, which I use when mapping human reads. HISAT2 is also threaded, so I run it on the files one at a time and pipe through samtools to create a sorted BAM file:

mkdir mapped
for r1 in trimmed/*_P1.fastq.gz;
  do base=$(basename "${x}" _P1.fastq.gz);
  r2="trimmed/${base}_P2.fastq.gz";
  echo ${base};
  hisat2 -p 20 -t -x /data/all/genomes/hsap/hisat_grch38_snp/genome_snp -1 \
    "${r1}" -2 "${r2}" 2&gt;"mapped/hisat2_${y}_vs_grch38_stderr.txt" | \
    samtools sort &gt; "mapped/hisat2_${y}_vs_grch38.bam";
done


3. Variant Calling

I currently use samtools/bcftools to do this, but would be interested in other opinions as we recently had a grant reviewer response that samtools was a dinosaur program and better approaches were available. Samtools doesn't currently work with threads for variant calling, so I save the commands to a text file and then run them through GNU Parallel. This requires the genome files to be first downloaded and indexed via samtools faidx. This produces a set of gzipped VCF files in the variants directory:

mkdir -p variants;
(for x in mapped/*.bam;
   do echo samtools mpileup -v -f /data/all/genomes/hsap/hisat_grch38_snp/Homo_sapiens.GRCh38.dna.primary_assembly.fa "${x}" \| \
   bcftools call --ploidy GRCh38 -v -m -O z -o variants/$(basename "${x}" .bam).vcf.gz; done) &gt; call_jobs.txt
cat call_jobs.txt | parallel -u;


As a slow, but more accurate alternative, variants can be called for all samples at the same time:

samtools mpileup -v -f /data/all/genomes/hsap/hisat_grch38_snp/Homo_sapiens.GRCh38.dna.primary_assembly.fa mapped/*.bam | \
  bcftools call --ploidy GRCh38 -v -m -O z -o variants/hisat2_allCalled_vs_grch38.vcf.gz


If a faster processing is desired, then the variant calling can be parallelised across chromosomes using parallel and merged afterwards using bcftools norm. The implementation of this is left as an exercise to the reader.