Call inversions with scoreInvHap
Carlos Ruiz, Juan R. Gonzalez
2019-10-29
Package
scoreInvHap 1.8.0
Contents
1 Introduction
scoreInvHap infers haplotype inversion’s status of a set of samples using SNP data. This method computes a similarity score between the sample SNPs in your cohort and the reference haplotypes. Samples are classified into the haplotype having the highest score. There are other approaches to perform this task such us inveRsion, invClust or PFIDO that are based on multivariate clustering procedures. However, these approaches may have several limitation that are addressed by using scoreInvHap. These limitations include:
- analizing limited or moderate number of samples. Clustering methods require to have large number of samples to perform clustering to properly discriminate haplotype inversion status. In particular, calling of a single sample cannot be performed without having a reference data.
- analyzing data from different SNP arrays. PFIDO and invClust may provide different answers depending on which array is used (i.e., Illumina 1M, Affy 6.0, …) since dfferent coverage may lead to perform the calling using different number of SNPs of a given inverted region.
- combining calling from different SNP arrrays. invClust and PFIDO clasify each individual as homozygous normal, heterozygous or homozygous inverted depending on cluster frequency. In those cases where inversion with haplotype frequency is close to 50%, both methods may have problems to determine which is the reference haplotype. Performing calling in different cohort of samples may lead to wrong results when combining data from different studies.
scoreInvHap overcomes these difficulties by using a set of reference genotypes that has been used to properly characterize the inversion of interest. Genomic information such as linkage desequillibrium (R2) and hetgerozygous sequence has been determined for each of the region of interest (ROI) and are used to create the score that discriminate the inversion status.
2 Inversion characterization
The package can be loaded by typing:
library(scoreInvHap)Calling procedure requires three objects that characterize the inversion of interest:
- the reference genotypes,
- the SNPs R2, and
- the heterozygous reference.
scoreInvHap already contains these required objects of four inversions. Two of them are well known and characterized genomic inversions: 8p23 and 17q21.31. We have also included two inversions that are described in invFEST database: 7p11.2 and Xq13.2. These objects have been created using VCF phase 3 data of 1000 Genomes Project. The code used to generate them can be found in the /inst/scripts folder of the package. This code can be modified to create the required files of any other inversion of interest. Required information of other inversions will be incorporated in the future when checked in our group or described in the literature.
2.1 Reference genotypes
As previously stated, the method uses the frequency of the SNP genotypes of each SNP located into the inversion region. This information is provided for he different haplotype populations. This information is encoded in an object called Refs that can be loaded by typing:
data("Refs")
names(Refs)## [1] "inv8p23.1" "inv17q21.31" "inv7p11.2" "invXq13.2"Each element of this list is a reference for one of the four available inversions. For instance, the reference of inversion inv8p23.1 can be obtained by:
ref <- Refs$inv8p23.1
class(ref)## [1] "list"ref[1:2]## $rs141039449
##
## haplos CC CG GG
## NI/NI 0.9310345 0.06321839 0.005747126
## NI/I 0.8595745 0.12765957 0.012765957
## I/I 0.6914894 0.25531915 0.053191489
##
## $rs138092889
##
## haplos AA AG GG
## NI/NI 0.005747126 0.02298851 0.9712644
## NI/I 0.012765957 0.05531915 0.9319149
## I/I 0.042553191 0.19148936 0.7659574This object is a list of matrices containing the frequency of each genotype (columns) in each inversion genotype (rows). Each component is named with the SNP id contained in the ROI. Notice that the alleles of the heteryzogous genotypes are alphabetically ordered.
2.2 SNPs R2 with inversion
The second object required is a vector containing the R2 between the inversion status and the SNPs genotypes. The SNPs with higher R2 will have more influence when computing the similarity score. We have formatted this object as a numeric vector, named with the SNPs ids. This information is provided in the object called SNPsR2. This is a list that contains the R2 of the four available inversions. For instance, this information for the inversion inv8p23.1 can be get by:
data("SNPsR2")
names(SNPsR2)## [1] "inv8p23.1" "inv17q21.31" "inv7p11.2" "invXq13.2"R2s <- SNPsR2$inv8p23.1
head(R2s)## rs141039449 rs138092889 rs138217047 rs143682252 rs550257716
## 0.0378562621 0.0289266440 0.0264322431 0.0264322431 0.0005282751
## rs2955571
## 0.00045703362.3 Heterozygous reference
The last required information is the heterozygous genotypes of the SNPs included in the references. This information is used to ensure that input SNPs have the same coded alleles than those used in the reference. This information can be retrieve from the object called hetRefs that can be inspect by typing:
data("hetRefs")
names(hetRefs)## [1] "inv8p23.1" "inv17q21.31" "inv7p11.2" "invXq13.2"hRefs <- hetRefs$inv8p23.1
head(hRefs)## rs141039449 rs138092889 rs138217047 rs143682252 rs550257716 rs2955571
## "CG" "AG" "AC" "AG" "CT" "AG"In that case, hRefs is a character vector contaning the heterozygous genotypes of the SNPs used as references in the ROI. It should be noticed that, in the heterozygous genotype, the alleles MUST BE ordered ALPHABETICALLY.
3 Running scoreInvHap: calling inversions
scoreInvHap deals with data either in a SNPMatrix or as Bioconductor VCF class. In the case of SNPMatrix, a list with two elements is required:
- genotypes: a SNPMatrix with individuals in rows and SNPs in columns
- map: a data.frame with the SNPs annotation. It must contain the columns allele.1 and allele.2 with the alleles of the SNPs.
We can load our data from a ped file or from plink binary format (.bed, .bim) to a SNPMatrix using snpStats:
library(snpStats)
## From a bed
snps <- read.plink("example.bed")
## From a pedfile
snps <- read.pedfile("example.ped", snps = "example.map")In both cases, snps is a list containing the elements genotypes and map. This object can be passed to scoreInvHap functions.
We can load a vcf file into a VCF object using the VariantAnnotation. We have included a small vcf in scoreInvHap package to illustrate how to deal with this data. This file contains a subset of SNPs of 30 European individuals belonging to the 1000 Genomes project. All these SNPs are located at the region 7p11.2. This vcf file contains imputed data. We can load the vcf with the following code:
library(VariantAnnotation)
vcf_file <- system.file("extdata", "example.vcf", package = "scoreInvHap")
vcf <- readVcf(vcf_file, "hg19")
vcf## class: CollapsedVCF
## dim: 380 30
## rowRanges(vcf):
## GRanges with 5 metadata columns: paramRangeID, REF, ALT, QUAL, FILTER
## info(vcf):
## DFrame with 4 columns: AF, MAF, R2, ER2
## info(header(vcf)):
## Number Type Description
## AF 1 Float Estimated Alternate Allele Frequency
## MAF 1 Float Estimated Minor Allele Frequency
## R2 1 Float Estimated Imputation Accuracy
## ER2 1 Float Empirical (Leave-One-Out) R-square (available only fo...
## geno(vcf):
## SimpleList of length 3: GT, DS, GP
## geno(header(vcf)):
## Number Type Description
## GT 1 String Genotype
## DS 1 Float Estimated Alternate Allele Dosage : [P(0/1)+2*P(1/1)]
## GP 3 Float Estimated Posterior Probabilities for Genotypes 0/0, ...We can observe that the object vcf contains 380 SNPs and 30 samples. Now we are ready to classify the samples with regard to the inversion inv7p11.2 by using the function scoreInvHap. This function requires four pieces of information:
- sample genotypes and their annotation (argument
SNPlist), - the list of matrices with the frequency of each genotypes in each inversion population (argument
Refs), - the numeric vector with the SNPs R2 (argument
SNPsR2), and - the vector with the heterozygous genotypes of the SNPs (argument
hetRefs).
The inv7p11.2 inversion status of the 30 samples from 1000 genomes is obtained by:
res <- scoreInvHap(SNPlist = vcf, SNPsR2 = SNPsR2$inv7p11.2,
hetRefs = hetRefs$inv7p11.2, Refs = Refs$inv7p11.2)
res## scoreInvHapRes
## Samples: 30
## Genotypes' table:
## IaIa IaIb IbIb NaIa NaIb NaNa NaNb NbIa NbIb NbNb
## 2 0 4 5 5 0 4 3 4 3
## - Inversion genotypes' table:
## NN NI II
## 7 17 6
## - Inversion frequency: 48.33%The results of scoreInvHap are encapsulated in a object of class scoreInvHapes. This object contains the classification of the samples and the simmilarity scores. We can obtain this data with the following getters:
# Get classification
head(classification(res))## HG00096 HG00097 HG00099 HG00100 HG00101 HG00102
## IbIb IbIb NaIb NbIb IaIa NbNb
## Levels: IaIa IaIb IbIb NaIa NaIb NaNa NaNb NbIa NbIb NbNb# Get scores
head(scores(res))## IaIa IaIb IbIb NaIa NaIb NaNa
## HG00096 0.7636176 0.7659578 0.9945596 0.6346131 0.7729563 0.7725742
## HG00097 0.7636176 0.7659578 0.9945596 0.6346131 0.7729563 0.7725742
## HG00099 0.6438462 0.7267291 0.7685422 0.7768338 0.9758600 0.7703418
## HG00100 0.3431981 0.4263283 0.4892990 0.2604141 0.4764038 0.3691768
## HG00101 0.9841178 0.7603894 0.7584714 0.7188997 0.6280646 0.7182172
## HG00102 0.2518848 0.1732271 0.3167020 0.1441663 0.1994766 0.3004813
## NaNb NbIa NbIb NbNb
## HG00096 0.4156252 0.3784514 0.5670885 0.3911247
## HG00097 0.4156252 0.3784514 0.5670885 0.3911247
## HG00099 0.5317241 0.3025523 0.5345369 0.2689791
## HG00100 0.6083141 0.6376241 0.6834476 0.5605652
## HG00101 0.3558992 0.4545927 0.4165003 0.3266407
## HG00102 0.5812806 0.5527802 0.5480628 0.72239223.1 Quality control
We can retrieve other values that are useful to evaluate the quality of the classification. For each sample, we can obtain its highest similarity score and the difference between the highest similarity score and the second highest:
# Get max score
head(maxscores(res))## HG00096 HG00097 HG00099 HG00100 HG00101 HG00102
## 0.9945596 0.9945596 0.9758600 0.6834476 0.9841178 0.7223922# Get difference score
head(diffscores(res))## HG00096 HG00097 HG00099 HG00100 HG00101 HG00102
## 0.2216032 0.2216032 0.1990263 0.0458235 0.2237283 0.1411116A classification is good when the highest score is close to 1 and the other scores are small. This means that the samples SNPs are almost the same than in one of the reference haplotypes and that they are different to the other references. We use the difference between the highest score and the second highest score as a measure of how different is the highest score from the rest. We can have a visual evaluation of these quality parameters with the function plotScores:
plotScores(res, pch = 16, main = "QC based on scores")
The horizontal line of the plot is a suggestion of the minimum difference between the highest and the second score that we accept. By default, this value is set to 0.1 but it can be changed with the parameter minDiff. This default value equals to considering that the sample SNPs are at least 10% more similar to the top reference than to the other references. plotScores relies on the base plot function, so we can pass additional parameters to customize the plot.
The other quality control estimate are based on the number of SNPs used in the computation. scoreInvHap allows having some missing measurements in the input data. However, this measurements are excluded from the computation of the scores. To reflect this issue, we have two measurements: the number of SNPs used in the computation and the proportion of non-missing measurements, or call rate:
# Get Number of scores used
head(numSNPs(res))## HG00096 HG00097 HG00099 HG00100 HG00101 HG00102
## 307 307 307 307 307 307# Get call rate
head(propSNPs(res))## HG00096 HG00097 HG00099 HG00100 HG00101 HG00102
## 1 1 1 1 1 1The number of SNPs must always be taken into account to evaluate the performance of the computation. It is highly recommended to use, at least, 15 SNPs in the computation. We have also included the function plotCallRate to plot the call rate of the samples:
plotCallRate(res, main = "Call Rate QC")
The vertical line is the minimum recommended call rate to consider a sample. By default, it is set to 0.9 but can be changed with the parameter callRate. Again, plotCallRate relies on the base plot function, so we can customize the plot.
The function classification have two parameters that selects samples based on these two QC parameters. The argument minDiff sets the minimum difference between the highest and the second highest score to consider a sample. The argument callRate sets the minimum call rate of a sample to pass the QC. By default, both arguments are set to 0 so all the sample are included:
## No filtering
length(classification(res))## [1] 30## QC filtering
length(classification(res, minDiff = 0.1, callRate = 0.9))## [1] 27Finally, the function classification has the argument inversion that, when it is set to TRUE, the haplotype based classification is transformed to an inversion based classification. This is useful for inversions inv7p11.2 and invXq13.2 that have more than one haplotype per inversion status:
## No filtering
table(classification(res))##
## IaIa IaIb IbIb NaIa NaIb NaNa NaNb NbIa NbIb NbNb
## 2 0 4 5 5 0 4 3 4 3## QC filtering
table(classification(res, inversion = TRUE))##
## II NI NN
## 6 17 73.2 Imputed data
When SNPs data are imputed, we obtain three different types of results: the best-guess, the dosage and the posterior probabilities. By default, scoreInvHap use the best-guess when computing the simmilarity scores. However, we can also use posterior probabilities setting the argument imputed to TRUE:
res_imp <- scoreInvHap(SNPlist = vcf, SNPsR2 = SNPsR2$inv7p11.2, hetRefs = hetRefs$inv7p11.2, Refs = Refs$inv7p11.2, imputed = TRUE)
res_imp## scoreInvHapRes
## Samples: 30
## Genotypes' table:
## IaIa IaIb IbIb NaIa NaIb NaNa NaNb NbIa NbIb NbNb
## 2 0 4 5 5 0 4 3 4 3
## - Inversion genotypes' table:
## NN NI II
## 7 17 6
## - Inversion frequency: 48.33%In this case, the samples were identically classified in both cases:
table(PostProbs = classification(res_imp),
BestGuess = classification(res))## BestGuess
## PostProbs IaIa IaIb IbIb NaIa NaIb NaNa NaNb NbIa NbIb NbNb
## IaIa 2 0 0 0 0 0 0 0 0 0
## IaIb 0 0 0 0 0 0 0 0 0 0
## IbIb 0 0 4 0 0 0 0 0 0 0
## NaIa 0 0 0 5 0 0 0 0 0 0
## NaIb 0 0 0 0 5 0 0 0 0 0
## NaNa 0 0 0 0 0 0 0 0 0 0
## NaNb 0 0 0 0 0 0 4 0 0 0
## NbIa 0 0 0 0 0 0 0 3 0 0
## NbIb 0 0 0 0 0 0 0 0 4 0
## NbNb 0 0 0 0 0 0 0 0 0 34 Other features
There are two additional parameters of scoreInvHap than can reduce computing time: R2 and BPPARAM. R2 indicates the minimum R2 that a SNP should have with the inversion to be included in the score. The less number of SNPs the less time it takes. By default, all SNPs are included in the computation. On the other hand, BPPARAM requires an instance of BiocParallelParam, which allows to parallelize the computation of the score. You can find more information about this class in its help page (?bpparam) and in the BiocParallel vignette.
sessionInfo()## R version 3.6.1 (2019-07-05)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.3 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.10-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.10-bioc/R/lib/libRlapack.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=C
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats4 parallel stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] snpStats_1.36.0 Matrix_1.2-17
## [3] survival_2.44-1.1 VariantAnnotation_1.32.0
## [5] Rsamtools_2.2.0 Biostrings_2.54.0
## [7] XVector_0.26.0 SummarizedExperiment_1.16.0
## [9] DelayedArray_0.12.0 BiocParallel_1.20.0
## [11] matrixStats_0.55.0 Biobase_2.46.0
## [13] GenomicRanges_1.38.0 GenomeInfoDb_1.22.0
## [15] IRanges_2.20.0 S4Vectors_0.24.0
## [17] BiocGenerics_0.32.0 scoreInvHap_1.8.0
## [19] BiocStyle_2.14.0
##
## loaded via a namespace (and not attached):
## [1] httr_1.4.1 splines_3.6.1
## [3] bit64_0.9-7 assertthat_0.2.1
## [5] askpass_1.1 BiocManager_1.30.9
## [7] BiocFileCache_1.10.0 blob_1.2.0
## [9] BSgenome_1.54.0 GenomeInfoDbData_1.2.2
## [11] yaml_2.2.0 progress_1.2.2
## [13] pillar_1.4.2 RSQLite_2.1.2
## [15] backports_1.1.5 lattice_0.20-38
## [17] glue_1.3.1 digest_0.6.22
## [19] htmltools_0.4.0 XML_3.98-1.20
## [21] pkgconfig_2.0.3 biomaRt_2.42.0
## [23] bookdown_0.14 zlibbioc_1.32.0
## [25] purrr_0.3.3 tibble_2.1.3
## [27] openssl_1.4.1 GenomicFeatures_1.38.0
## [29] magrittr_1.5 crayon_1.3.4
## [31] memoise_1.1.0 evaluate_0.14
## [33] tools_3.6.1 prettyunits_1.0.2
## [35] hms_0.5.1 stringr_1.4.0
## [37] AnnotationDbi_1.48.0 compiler_3.6.1
## [39] rlang_0.4.1 grid_3.6.1
## [41] RCurl_1.95-4.12 rappdirs_0.3.1
## [43] bitops_1.0-6 rmarkdown_1.16
## [45] DBI_1.0.0 curl_4.2
## [47] R6_2.4.0 GenomicAlignments_1.22.0
## [49] knitr_1.25 dplyr_0.8.3
## [51] rtracklayer_1.46.0 bit_1.1-14
## [53] zeallot_0.1.0 stringi_1.4.3
## [55] Rcpp_1.0.2 vctrs_0.2.0
## [57] dbplyr_1.4.2 tidyselect_0.2.5
## [59] xfun_0.10