Hands On Session: Exploring genotyping error with R package ‘whoa’

Introduction

Now we all have the chance to use Eric’s R-package whoa to estimate heterozygote miscall rates from GBS data stored in VCF files.

This is a small, lightweight package that lets users investigate the distribution of genotypes in genotype-by-sequencing (GBS) data where they expect (by and large) Hardy-Weinberg equilibrium, in order to assess rates of genotyping errors and the dependence of those rates on read depth.

The name comes from the bolded letters in this sentence:

Where’s my Heterozygotes at? Observations on genotyping Accuracy.

It also fits well with Eric’s reaction when he started investigating heterozygote miscall rates (rates at which true heterozygotes are incorrectly called as homozygotes) in some RAD-seq data sets—His eyes bugged out and he said, “Whoa!”

Installing whoa

Unfortunately, CRAN was shuttered this last week, so we were not able to get this package onto CRAN. So, there are a few options:

Option 1: From github

If you have all the tools to compile R packages that use C++ code on your computer then this option is fine.

devtools::install_github(repo = "eriqande/whoa")

Option 2: Precompiled versions for Mac and Windows

If you are working on a Linux machine than you probably have all the tools for compiling packages with C++. However, not all Mac and Windows users will have that capacity, so I have some precompiled binary packages for those users.

First, you will have to ensure that all the dependencies (necessary other packages) are also installed. If you think you don’t have the packages, you can install them easily with:

install.packages(c("tidyverse", "Rcpp", "vcfR", "viridis"))

Then, depending on if you are a Mac or a Windows user there are two further options.

Installing the binary for Mac OS X

# download the package 
download.file("https://www.dropbox.com/s/cnc6rfqe8afp8fd/whoa_0.0.1.tgz?dl=1", 
              destfile = "whoa_0.0.1.tgz")

# then install the package
install.packages("whoa_0.0.1.tgz", repos = NULL)

Installing the binary for Windows

# download the package 
download.file("https://www.dropbox.com/s/f48h75oft1rlzr8/whoa_0.0.1.zip?dl=1", 
              destfile = "whoa_0.0.1.zip")

# then install the package
install.packages("whoa_0.0.1.zip", repos = NULL)

A first run through

The package comes with a small bit of data from lobster to play with. The rest of this document shows a quick run through a few of the functions to do an analysis of a data set.

Packages

# load up the package:
library(whoa)

Lobster data

Read about the lobster data here. Execute this if you want:

help("lobster_buz_2000")

The main thing to know is that it is a vcfR object. You can make such an object yourself by reading in a VCF file using vcfR::read.vcfR().

Make a quick genotype frequency scatter plot

# first get compute expected and observed genotype frequencies
gfreqs <- exp_and_obs_geno_freqs(lobster_buz_2000)
# then plot those.  Set max_plot_loci so that all 2000
# loci will be plotted
geno_freqs_scatter(gfreqs, max_plot_loci = 2000)

Now infer an overall heterozygote miscall rate.

If we want to estimate the het miscall rate (over all read depth bins) we just set the minimum bin size to a very large value so it make just one bin:

overall <- infer_m(lobster_buz_2000, minBin = 1e15)

Preparing data structures for MCMC
Running MCMC
Tidying output

Now look at that:

overall$m_posteriors

Wow! (Or should we say “WHOA!”) A het miscall rate of around 25%.

Now infer a miscall rate for read depth bins

See the total_n above is about 65,000. That means 65,000 genotypes. (2000 loci typed at 36 individuals, with some missing data).
We will bin those up so that there are at least 2000 genotypes in each bin and then estimate the het miscall rate for each read depth bin.

binned <- infer_m(lobster_buz_2000, minBin = 2000)

Preparing data structures for MCMC
Running MCMC
Tidying output

And then we can plot the posterior mean and CIs for each read depth bin.

posteriors_plot(binned$m_posteriors)

Again, WHOA! The het miscall rate at low read depths is super high!

Another analysis from a VCF file

Now, let’s look at doing this while beginning with a VCF file.

For those that aren’t familiar with VCF files (that stands for Variant Call Format), you will want to get to know the format well. It is one of the main standards for storing information about SNPs (and other variants) that were obtained from sequencing. It allows the storage of extra information like variant quality scores and read depths. You can read more about VCF formats on Wikipedia.

Download a VCF file

We can get the red-drum data like this:

download.file("https://www.dropbox.com/s/twyylsui15q65yj/red_drum_Final_Filtered_SNPs.vcf.gz?dl=1",
              destfile = "red_drum_Final_Filtered_SNPs.vcf.gz")

This is a gzipped file, so it is a little hard to look at. Here is what the first 20 lines of the file look like:

gzcat red_drum_Final_Filtered_SNPs.vcf.gz | head -n 20

##fileformat=VCFv4.1
##source=freeBayes v0.9.20
##reference=reference.fasta
##phasing=none
##filter="AB > 0.25 | AB < 0.01 genotypes filtered with: QR > 0 | QA > 0 "
##filter="QUAL / DP > 0.2"
##filter="MQM / MQMR > 0.25 & MQM / MQMR < 1.75"
##filter="SAF / SAR > 100 & SRF / SRR > 100 | SAR / SAF > 100 & SRR / SRF > 50"
##filter="PAIRED > 0.05 & PAIREDR > 0.05 & PAIREDR / PAIRED < 1.75 & PAIREDR / PAIRED > 0.25 | PAIRED < 0.05 & PAIREDR < 0.05"
##INFO=<ID=NS,Number=1,Type=Integer,Description="Number of samples with data">
##INFO=<ID=DP,Number=1,Type=Integer,Description="Total read depth at the locus">
##INFO=<ID=DPB,Number=1,Type=Float,Description="Total read depth per bp at the locus; bases in reads overlapping / bases in haplotype">
##INFO=<ID=AC,Number=A,Type=Integer,Description="Total number of alternate alleles in called genotypes">
##INFO=<ID=AN,Number=1,Type=Integer,Description="Total number of alleles in called genotypes">
##INFO=<ID=AF,Number=A,Type=Float,Description="Estimated allele frequency in the range (0,1]">
##INFO=<ID=RO,Number=1,Type=Integer,Description="Reference allele observation count, with partial observations recorded fractionally">
##INFO=<ID=AO,Number=A,Type=Integer,Description="Alternate allele observations, with partial observations recorded fractionally">
##INFO=<ID=PRO,Number=1,Type=Float,Description="Reference allele observation count, with partial observations recorded fractionally">
##INFO=<ID=PAO,Number=A,Type=Float,Description="Alternate allele observations, with partial observations recorded fractionally">
##INFO=<ID=QR,Number=1,Type=Integer,Description="Reference allele quality sum in phred">

gzcat: error writing to output: Broken pipe
gzcat: red_drum_Final_Filtered_SNPs.vcf.gz: uncompress failed

Read in the VCF file

We use read.vcfR() from the vcfR package for this

drum_vcf <- vcfR::read.vcfR("red_drum_Final_Filtered_SNPs.vcf.gz")

Scanning file to determine attributes.
File attributes:
  meta lines: 59
  header_line: 60
  variant count: 7382
  column count: 214

Meta line 59 read in.
All meta lines processed.
gt matrix initialized.
Character matrix gt created.
  Character matrix gt rows: 7382
  Character matrix gt cols: 214
  skip: 0
  nrows: 7382
  row_num: 0

Processed variant 1000
Processed variant 2000
Processed variant 3000
Processed variant 4000
Processed variant 5000
Processed variant 6000
Processed variant 7000
Processed variant: 7382
All variants processed

Analysis

From this point on, its just the same as we did with the lobster data:

# first get compute expected and observed genotype frequencies
drumfreqs <- exp_and_obs_geno_freqs(drum_vcf)
# then plot those.  Set max_plot_loci so that all 2000
# loci will be plotted
geno_freqs_scatter(drumfreqs, max_plot_loci = 5000)

That looks a lot better…

Now, how about an overall het miscall rate?

drum_overall <- infer_m(drum_vcf, minBin = 1e15)

Preparing data structures for MCMC
Running MCMC
Tidying output

drum_overall$m_posteriors

That takes a bit longer (more individuals), but the results are good—about a 5% het miscall rate.

But note that the miscall rate is clearly higher at lower read depths.

First, check how many individuals and how many loci we have here:

dim(drum_vcf@gt)

[1] 7382  206

So, 7382 loci and 205 individuals. That means close to 1.5 million genotypes. So, if we want to break that up into read depth bins, we could put 50,000 in each bin and still have a large number of bins:

drum_binned <- infer_m(drum_vcf, minBin = 50000)

Preparing data structures for MCMC
Running MCMC
Tidying output

posteriors_plot(drum_binned$m_posteriors)

And we see a clear trend there.

Now, please use your own data set

If you have RAD data in vcf format lying around, please run it through whoa!

Note that if you have multiple, genetically distinct populations in your VCF, you can select individuals from just one population by indexing them out of the vcfR object. For example, if we wanted only a subset of samples from the drum VCF file we could do like this:

sams <- c("AR_001", 
          "AR_003",
          "AR_004",
          "AR_005",
          "AR_008",
          "AR_010",
          "AR_012",
          "AR_014",
          "AR_015",
          "AR_016")
drum_subset <- drum_vcf[, c("FORMAT", sams)]
# check the dimensions of the genotypes in that object:
dim(drum_subset@gt)

[1] 7382   11

Notice how you have to include the column “FORMAT” when you index the object. This is critical—that is the column that tells you how all the auxillary information that comes with the genotypes (like read depths and quality scores) is formatted.