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-<A-<FFFJA-<-A7-F<<FFJJAJJJJJJJJJ--<--7-7AJ7<AAJA--J7<ACTGCTGTGGGGCACCCAGCCCCCCAGATAGCCTGGCAGAAGGATGGGGGCACAGACTTCCCAGCTGCACGGGAGAGAC'
at org.usadellab.trimmomatic.fastq.FastqRecord.<init>(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>"mapped/hisat2_${y}_vs_grch38_stderr.txt" | \
samtools sort > "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) > 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.
Aside: I see that GATK's best practise guidelines call per-sample with HaplotypeCaller, then process those VCF files as aggregated data. @Geraldine-vdauwera says this two-step joint calling workflow was designed to be scalable for any size cohort, so you can apply it pretty much the same way whether you have two samples or twenty thousand. Once you go higher than that there are some additional tweaks to keep scaling up.
