In edsgard/scphaser: Phase Variants Within Genes Using Allele Counts

Create data-structure

Create an "acset" data-structure. An acset is a list that at a minimum contains four data-structures:

featdata: A data-frame which maps variant names to feature names and which contains reference and alternative alleles. It must contain four columns named "feat", "var", "ref" and "alt", and the rownames must be set to the "var" column.
phenodata: A data-frame which annotates the samples. It must contain a column named "sample". If not provided it is created by the "new_acset" function.
refcount: A matrix with allelic counts for the reference allele. The rownames have to match the values in the "var" column in the featdata, and the colnames the values in the phenodata "sample" column.
altcount: A matrix with allelic counts for the alternative allele, with the same row- and col-names as refcount.

As part of this package a dataset with allele counts from human are provided, called "marinov". The dataset is a list containing the four required data-structures to create an "acset".

##Extract data-structures from the marinov list
library('scphaser')
invisible(marinov)
featdata = marinov[['featdata']]
refcount = marinov[['refcount']]
altcount = marinov[['altcount']]
phenodata = marinov[['phenodata']]

Create an "acset" data-structure

acset = new_acset(featdata, refcount, altcount, phenodata)

Print the elements of the acset data-structure and the dimensions of these, illustrating the number of variants and cells.

lapply(acset, dim)

Print the number of genes

length(unique(acset$featdata$feat))

Filter on at least n variants per feature

Features with less than "nminvar" variants are removed

nminvar = 2
acset = filter_feat_nminvar(acset, nminvar)

Print the element dimensions after filtering out genes with less than 2 variants

lapply(acset, dim)

Call genotypes

Transcribed genotypes are called as 2 or 0 if there at least "min_acount" reads and the fold-change >= 3 or <= 1/3, where fold-change = alternative allele count / reference allele count. For entries that do not meet the criteria, such as when bi-allelic expression close to a 50/50 expression ratio between the alleles, the genotype is set to 1.

min_acount = 3
fc = 3
acset = call_gt(acset, min_acount, fc)

Print the elements after calling genotypes. Note that a new element called "gt" was created, containing the called genotypes.

lapply(acset, dim)

Randomize original counts. The randomized dataset will be used below.

acset_rnd = racset(acset, type = 'gt')

Filter variants

Filter variants on having at least "nmincells" cells with monoallelic calls in at least "nminvar" variants within a feature. After variant filtering, filter features on having at least "nminvar" variants.

nmincells = 5
nminvar = 2
acset = filter_acset(acset, nmincells, nminvar)

Number of variants after filtering

nrow(acset$featdata)

Number of features after filtering

length(unique(acset$featdata$feat))

Phase

There are three phasing arguments that can be provided to the phasing function, each with two possible values. Below we call it using genotypes, exhaustive clustering and without weighing and also illustrate calling using the other possible values (first commented lines).

##acset = phase(acset, input = 'ac', weigh = FALSE, method = 'exhaust')
##acset = phase(acset, input = 'ac', weigh = FALSE, method = 'pam')
##acset = phase(acset, input = 'ac', weigh = TRUE, method = 'exhaust')
##acset = phase(acset, input = 'ac', weigh = TRUE, method = 'pam')

##acset = phase(acset, input = 'gt', weigh = TRUE, method = 'exhaust')
##acset = phase(acset, input = 'gt', weigh = TRUE, method = 'pam')
##acset = phase(acset, input = 'gt', weigh = FALSE, method = 'pam')
acset = phase(acset, input = 'gt', weigh = FALSE, method = 'exhaust', verbosity = 0)

As an overview of the elements in the acset datastructure we print the dimension of each element

lapply(acset, dim)

The haplotype output is contained in an element that was added by the phasing function and is named "phasedfeat"

head(acset[['phasedfeat']])

The variants that have been swapped compared to the input reference and alternative alleles are contained in "varflip". The number of swapped variants is the length of this vector.

length(acset[['varflip']])

Assess genotype concordance

To assess the success of the phasing one can calculate the degree of variability remaining if all cells with haplotype 2 are set to haplotype 1. As a rough measure of this we here calculate the variability as the number of cells that differ from the inferred haplotype for every gene with two variants. The differing number of cells per gene we here denote inconcordance.

##get gt concordance before and after phasing
acset = set_gt_conc(acset)

Inconcordance before phasing

table(acset$gt_conc$notconc$feat2ncell)

Inconcordance after phasing

table(acset$gt_phased_conc$notconc$feat2ncell)

Phasing of a randomized genotype matrix

Filtering

nmincells = 3
nminvar = 2
acset_rnd = filter_acset(acset_rnd, nmincells, nminvar)

Dimensions after filtering

lapply(acset_rnd, dim)

Number of genes after filtering

length(unique(acset_rnd$featdata$feat))

Phasing

acset_rnd = phase(acset_rnd, input = 'gt', weigh = FALSE, method = 'exhaust', verbosity = 0)

Get genotype matrix concordance before and after phasing

acset_rnd = set_gt_conc(acset_rnd)

Inconcordance before phasing

table(acset_rnd$gt_conc$notconc$feat2ncell)

Inconcordance after phasing

table(acset_rnd$gt_phased_conc$notconc$feat2ncell)

Plot

The removal of variability from the cell distribution after phasing indicates successful phasing. The right subfigure show that phasing cannot be done if the genotype matrix is scrambled, since the transcribed genotype of each cell then do not anylonger correspond to one of two underlying haplotype states.