Thanks to @tfwillems for helpful comments on this post!

Introduction

Our lab is very interested in genotyping short tandem repeats (STRs) from whole genome sequencing datasets. For us, the idea of PCR free libraries is extremely exciting, since it has the potential to greatly reduce one of our major sources of errors: stutter noise. Noise at STRs is thought to result from polymerase slippage during amplification, and results in many reads being off from the true allele by a unit (or more than one unit) of the repeat motif. Below are some example alignments where some reads are a result of stutter noise (ok ok, I also put these here as an excuse to try out/show off the new "export figure" feature of my pybamview tool, which you should check out if you often find yourself looking at sequence alignments for indels).

The examples above are from male chrX calls, so there is only one true allele present. In the first example, one read is off by one repeat unit from the true allele (it is missing an AGAT). In the second, there's an example of a stutter read off by 2 repeat units. We see examples such as these quite commonly. It is even worse at diploid loci where there are potentially 2 alleles, plus stutter, making it difficult to determine the true genotype. In our program lobSTR, we attempt to model this noise when scoring genotypes. But in low coverage data, such as the 1000 Genomes Project samples, where we often have only a single read for each call, the model's not going to help us much. If you have one read, and that read is from stutter, you're almost certainly going to make a wrong call at that locus.

The PCR-free method has the potential to dramatically reduce the PCR stutter problem at STRs. Below I compare two whole genome sequencing datasets of similar quality, one PCR free and one not, and show evidence that the PCR free protocol cleans up a lot of the noise at STRs, reducing the stutter rate by several fold.

Datasets and STR genotyping

• HGDP00542 (with PCR): This dataset was obtained from the CEPH website and can be downloaded here. This is a Papuan individual sequenced to 25.9x. The dataset was published with the Meyer et al. Denisova paper in Science 2012.
• High coverage PCR free dataset (This dataset is not yet published). This is another HGDP sample sequenced to 37.6x coverage using the PCR-free protocol.

Both datasets consist of 100bp paired end reads for male samples. Besides the difference in the PCR free protocol vs. not, these datasets are expected to have similar quality.

I used our STR genotyping tool, lobSTR to genotype STRs in these samples. For alignment, I used v2.0.4. For allelotyping, I used v2.0.5. For anyone that wants all the gory details the specific commands I used are at the bottom of this post.

PCR-free dataset had 4-fold less "stutter" reads

For each sample, I looked at loci on chrX to estimate how many reads are the result of stutter noise. I took all loci with at least 5x coverage and a unique modal allele (the allele supported by the most reads) which I called the true allele. Since the samples are male, chrX calls should be hemizygous, and any read not matching the true genotype is from stutter (or an alignment error, which definitely happens, but I'm ignoring that here). This is the same method lobSTR uses to train its stutter model.

The first metric I looked at was simply the percentage of reads determined to be from stutter noise. Below are the results broken up by motif size. The stutter rate clearly decreases with motif length (homopolymers are by far the worst, and the bane of anyone's existence who works with STR calling, whereas longer motifs have very few reads from stutter). On average, the PCR free dataset had a stutter rate more than 4-fold lower than the dataset using the non-PCR-free protocol! Pretty encouraging.

While the rate of stutter reads was greatly reduced in the PCR-free dataset, the distribution of stutter error sizes was about the same. Below is the distribution of error sizes for each motif length (I only showed periods 2-5 to save space, but the plots for periods 1 and 6 are similar). You can see what we have observed many times: stutter noise tends to come in multiples of the repeat unit length and has a strong bias toward decreasing the number of repeats from the true allele. This is true for both datasets, perhaps with the PCR-free dataset having a slightly lower bias toward decreasing allele length.

Another metric I looked at was the quality score of STR calls. For each genotype, lobSTR reports a quality score. This score is calculated as:

A \(Q\) value close to 1 means a quite confident call, whereas a \(Q\) close to 0 means the call is not much better than the other possible genotype calls.

I compared the average \(Q\) for different coverage levels in the two datasets. The y-axis gives \(-1\log10(1-Q)\) (higher=better). The scores are significantly higher at all coverage levels in the PCR-free dataset:

Conclusion

In conclusion, the Illumina PCR-free protocol seems to have great promise for improving STR calls from high throughput sequencing datasets. While I only compared two datasets here, we see similar trends in other datasets: PCR free samples have much lower stutter noise at STRs. Of particular note is the improvement at homopolymers, which have until now been nearly impossible to call with any reasonable accuracy. This holds great promise for genotyping STRs as more and more sequencing datasets are generating using this protocol.

Appendix - lobSTR commands

 
# lobSTR alignment
lobSTR \
  -f ${BAM_SORTED_BY_NAME} \
  --bampair \
  --index-prefix resource_bundles/hg19/lobstr_v2.0.3_hg19_ref/lobSTR_ \
  --rg-sample ${SAMPLE} \
  --rg-lib ${SAMPLE} \
  --out ${SAMPLE} \
  -p 20 \
  --fft-window-size 16 \
  --fft-window-step 4 \
  --bwaq 15 \
  -v
samtools sort ${SAMPLE}.aligned.bam ${SAMPLE}.sorted
samtools index ${SAMPLE}.sorted.bam

# Train noise model
samtools view -b ${SAMPLE}.sorted.bam chrX > train_${SAMPLE}.sorted.bam
samtools index train_${SAMPLE}.sorted.bam
allelotype \
  --command train \
  --noise_model ${SAMPLE} \
  --index-prefix resource_bundles/hg19/lobstr_v2.0.3_hg19_ref/lobSTR_ \
  --strinfo resource_bundles/hg19/lobstr_v2.0.3_hg19_strinfo.tab \
  --haploid chrX \
  --bam train_${SAMPLE}.sorted.bam \
  --min-border 5 \
  --min-bp-before-indel 7 \
  --maximal-end-match 15 -v

# Run allelotype
allelotype \
  --command classify \
  --noise_model ${SAMPLE} \
  --index-prefix resource_bundles/hg19/lobstr_v2.0.3_hg19_ref/lobSTR_ \
  --strinfo resource_bundles/hg19/lobstr_v2.0.3_hg19_strinfo.tab \
  --bam train_${SAMPLE}.sorted.bam \
  --min-border 5 \
  --min-bp-before-indel 7 \
  --maximal-end-match 15 -v \
  --out ${SAMPLE}
vcf-sort ${SAMPLE}.vcf | bgzip -c > ${SAMPLE}.vcf.gz
tabix -p vcf ${SAMPLE}.vcf.gz