SureTypeSCR: R package for rapid quality control and genotyping of SNP arrays from single cells

Metadata

The minimal set of input data for loading the Illumina SNP array data consists of a manifest file, cluster file, sample sheet and a set of GTC files, where each GTC file corresponds to one sample analyzed on the SNP array (Table 1). Both manifest file and cluster file are issued by Illumina per SNP array type. While manifest file describes SNP markers used on the array, cluster file contains information about genotype clusters per SNP marker gathered from population studies and used for scoring in the GenomeStudio software.¹³

Implementation

The core of the package is implemented in a Python library and SureTypeSCR communicates with this library using reticulate. SureTypeSCR further uses Illumina's Python library IlluminaBeadArray to load the GTC files and then utilizes functions from the tidyverse ecosystem (packages dplyr and magrittr) to implement functions for assessing data quality. The data classification process then assigns a quality score to each analyzed single cell genotype (Figure 1).

Figure 1. Workflow implemented within the SureTypeSCR package.

SureTypeSCR utilizes Illumina's IlluminaBeadArray library to load the metadata (Table 1) and raw genotype files. In case the data is in IDAT format, SureTypeSCR utilizes Illumina’s IAAP CLI software to convert it into GTC format. SureTypeSCR then implements various functions to check the quality of the data, perform intensity transformation, run dimension reduction algorithm and visualize the results. Subsequently, classification is performed using machine learning algorithm previously trained on large batch of single cell data with known ground truth.⁶ The algorithm is currently embodied in RF-GDA, which is part of the SureTypeSC Python library. An optional step is context dependent validation that can be implemented within SureTypeSCR in case parental or ploidy information are available.

Operation

R (>4.0) and Python (>3.6) are required for installing and running the SureTypeSCR package. The software is installed using devtools. To ensure maximal reproducibility across different platforms, a virtual Python environment is created and all necessary Python dependencies are installed in this environment using reticulate. Subsequently, SureTypeSCR is built, installed and linked to the Python virtual environment The package was tested on three major platforms (Linux/Win/Mac). Data processing times depend on the number of samples in the batch and is estimated at 20s per sample on a single CPU with 4 GB RAM.

Data initialization and QC

We start the analysis with initializing the package, data and metadata (see Table 1 and code below). The R data package containing the sperm data can be downloaded from GitHub using devtools. Function data(.) then initializes data frame metadata, which stores the family information and other metadata that can be used in the analysis and samplesheet containing path to the downloaded samplesheet with the data. Manifest and cluster file are part of the SuretypeSCR installation. Function scbasic(.) loads the data into an R data frame. We then filter out SNPs, termed intensity only SNPs, that are used to detect copy number variant but do not provide genotyping information (filter(.) and str_detect(.)).

library(devtools)
# install SureTypeSCR from github
devtools::install_github("Meiomap/SureTypeSCR")
# install data package with sperm data (compiled from GSE19247)
devtools::install_github("Meiomap/johnsonspermdata")

library(SureTypeSCR)
library(johnsonspermdata)
# load metadata and samplesheet location
data(metadata,samplesheet)
setwd(system.file("data",package="johnsonspermdata"))

manifest=system.file("files/HumanCytoSNP-12v2_H.bpm",package="SureTypeSCR")
cluster=system.file("files/HumanCytoSNP-12v2_H.egt",package="SureTypeSCR")

df = scbasic(samplesheet=samplesheet, bpm=manifest, egt=cluster)
#filtering out intensity-only SNPs
df %<>% filter ( Chr !=0 & !str_detect (Name ,"cnv"))

Calculating call rates per individual and genotype reveals a high degree of heterozygosity (Figure 2A, AB rates), suggestive of significant ADI, since sperm are haploid cells and there were no aneuploidies reported in these samples:⁸

df %>%
group_by(individual,gtype) %>%
callrate()

Figure 2. Visual outputs from the analysis of 23 sperm samples with SureTypeSCR.

(A) Call rates per individual and genotype calculated as proportions of called genotypes and all genotypes and no calls. No calls represent SNPs with itentisites below a baseline defined in the paper describing the original data.⁸ (B) PCA analysis accross all SNPs and genotypes from 23 samples that called in every sample. (C) MA plot on normalized intensities across 23 samples, the X axis corresponds to logarithmic average (a) and the Y axis is logarithmic difference (m). (D) MA plot across 23 samples after filtration with SureTypeSC. (E) Effect of used threshold on average heterozygosity (solid line) and average call rate (dashed line) across 23 sperm samples for both, SureTypeSC (grey line) and Illumina GenCall (yellow line). The ribbons represent the standard error of mean.

The principal component analysis (PCA) is performed using function plot_pca(.) that returns a ggplot object:

df %>%
plot_pca(features="gtype", metadata=metadata, by_chrom=FALSE)

As shown in the code example above, users can customize which features (columns of the data frame) to use for the PCA analysis with the features parameter. There is an option to customize and add metadata to the ggplot object (currently, family information is supported) and a choice whether the PCA should be run per chromosome (by_chrom parameter) or on the whole data frame. While the per chromosome analysis can reveal aneuploid chromosomes, the latter is useful for validating kinship of the samples. This is demonstrated in Figure 2B, where the 23 sperm samples are separated into two clusters corresponding to two families defined in the metadata.

Data transformation and classification

Transformation of the intensities into a logarithmic scale minimizes the variability between the SNPs and samples and allows the patterns of the genotyping clusters to be detected.⁶ In order to evaluate the single cell genotypes using our classification algorithm, we calculate the logarithmic difference and logarithmic average of the intensities (m and a, respectively, Figure 2C). The following code performs the data transformation by adding four additional columns to the original data frame, two for the raw intensities and two for the normalized intensities for the X and Y channels. The user can then control the plotting by adjusting the fraction of points to be visualized, whether a smoothing spline should be applied to the transformed data and whether to use normalized intensities for plotting (parameters n, smooth and normalized in plot_ma(.)).

df %>%
calculate_ma() %>%
plot_ma(n=0.1)

Note that, by default, plot_ma(.) visualizes the plots per sample and we use stat_bin_2d(.) in Figure 2C to illustrate the point density and error distribution across the whole dataset. The MA plot in Figure 2C reveals an erroneous heterozygous cluster where m is close to zero and a is low that we previously showed is due to ADI.⁶ We subsequently perform sample genotype classification with SureTypeSC using:

clf=system.file('files/rf.clf',package="SureTypeSCR")

df_model = df %>%
       calculate_ma() %>%
       group_by(individual) %>%
       nest() %>%
       mutate(model=map(data, function(df) suretype_model(df,individual, clf)))

The first layer of the classification algorithm (Random Forest) is loaded from the file. Then, the classification model is created per individual sample (group_by(.) and nest(.)) using Gaussian Discriminant Analysis to infer model parameters.⁶ The Gaussian Discriminant Analysis is conducted per individual sample rather than the combined dataset in order to avoid bias in the scoring function due to potential outliers in the data. The first two parameters of suretype_model(.) are formal and the last parameter defines the classifier (clf) to be used in the first layer (see the reference manual for a detailed description of all available parameters). After unnesting the df_model, the dataframe contains an additional column that contains the SureTypeSC classification score (rfgda_score). We can then apply a threshold (set_threshold(.)) and use MA plot again to observe how SureTypeSC has affected the quality of the data:

df_model_unnested = df_model %>%
unnest(c(data,model))

df_model_unnested %>%
set_threshold(clfcol='rfgda_score',threshold=0.5) %>%
plot_ma()

Figure 2D shows the results from the entire dataset (using stat_bin_2d(.)). Unlike Figure 2C, which contains the data prior to SureTypeSC, the heterozygous cluster (m close to 0 and low a) caused by ADI is effectively removed and the data are concentrated along m = 4 and m = −4 representing homozygous AA and homozygous BB genotypes, respectively.

Finally, we determined the call rate and % of heterozygous SNPs in the data as a function of the used threshold in both SureTypeSC and Illumina's GenCall (rfgda_score and score columns in the data frame, respectively):

callr_calc <- function(.data,algo,threshold)
{
  .data %>%
   set_threshold(clfcol=algo,threshold=threshold) %>%
   group_by(individual,gtype) %>%
   callrate() %>%
   pivot_wider(names_from=gtype, values_from=Callrate) %>%
   mutate(thr=threshold, alg=algo) %>%
   mutate(callrate=AA+BB+AB,thr=threshold,alg=algo)
}

thrs=seq(0.0,0.96,0.01)

suretype=(map(thrs,function(x) callr_calc(df_model_unnested,'rfgda_score',x))) %>% bind_rows()
gencall=(map(thrs,function(x) callr_calc(df_model_unnested,'score',x))) %>%
bind_rows()
performance=bind_rows(suretype,gencall)

Figure 2E confirms that SureTypeSC is more specific towards noise whilst retaining higher call rates as the threshold increases compared to GenCall. This is consistent with our validation study published previously.⁶