mpbr: Summary statistics

Source: vignettes/mpbr.Rmd 
library(mpbr)
#> 
#> Attaching package: 'mpbr'
#> The following object is masked from 'package:stats':
#> 
#>     filter

Construct the SNPdata object

The functions in this package require a SNPdata object. Objects of this class are generated using the get_snpdata() function with the following arguments:

  1. vcf_file: the input VCF file (required).
  2. meta_file: a tab-delimited file with the samples metadata (required).
  3. output_dir: the path to the folder where the output files will be stored (optional).
  4. gaf: a tab-delimited file with the gene ontology annotation (optional). If not provided, the default file that was downloaded from https://plasmodb.org/plasmo/app/downloads/Current_Release/Pfalciparum3D7/gaf/ will be used.
  5. gff: a tab-delimited file with the gene annotation file (optional). If not provided, the default file that was downloaded from https://plasmodb.org/plasmo/app/downloads/Current_Release/Pfalciparum3D7/gff/ will be used.
  6. num_threads: the number of threads to use when importing the input data.

Type ?get_snpdata for more details.

snpdata <- get_snpdata(
  vcf_file = vcf,
  meta_file = meta_file,
  output_dir = output_dir,
  gof = gof,
  gff = gff
)

Use the print() function to visualize the structure of the SNPdata object created from the get_snpdata().

print(snpdata)

Calculate the SNPs minor allele frequency

the compute_maf() requires a SNPdata object, the name of the genotype table from which the MAF should be calculated (either “GT” for the raw genotypes, “Phased” for phased genotypes, “Phased_Imputed” for phased and imputed genotypes).
include_het argument determines whether to account for heterozygous sites during MAF calculation. The function returns a SNPdata object with the following 2 extra columns in the details table:

  • MAF: the column with the SNPs minor allele frequencies
  • MAF_allele: the column that determines which of the reference (0) or alternate (1) alleles is the minor allele on that site.

type ?compute_MAF for more details.

snpdata <- compute_maf(
  snpdata,
  include_het = FALSE,
  mat_name = "GT"
)
head(snpdata$details)

Filter the SNPdata object

Use the filter_snps_samples() function to remove the SNPs and samples that do not satisfy the user criteria. In the example below, we intend to discard SNPs and samples where the number of missing genotypes represent 20% or more the total number of observations. SNPs with MAF < 1% or a quality score < 10 will also be removed.

snpdata <- filter(snpdata,
                  min_qual            = 10,
                  max_missing_sites   = 0.2,
                  max_missing_samples = 0.2,
                  maf_cutoff          = 0.01)

All loci and samples that do not satisfy the defined conditions will be dropped from the input SNPdata object. Details about the filtered object can be viewed using:

print(snpdata)

The distribution of the MAF in the filtered SNPdata object looks like below:

hist(snpdata$details$MAF, 100)
grid()

Mixed genotypes phasing

Mixed genotypes are not accounted for in many population genetics methods as they introduce a extra layer of complexity. Genotype phasing is the most common way of dealing with mixed genotype. In this package, we use the phase_mixed_genotypes() to phase the mixed genotypes based on their allelic depth (AD) and the total allele count on the specific locus through the following procedure:

  • When there is no read supporting both alleles, the mixed genotype is replaced by the less common allele in the locus. If both alleles are equally represented in the population, it is replaced based on a Bernoulli simulation with a probability defined as referenceallelecountreferenceallelecount+alternateallelecount\frac{reference_allele_count}{reference_allele_count + alternate_allele_count}
  • When both allele of the mixed genotype are supported by more than 5 reads and the number of reads supporting either the reference or or alternate allele is 2 times the other allele, it will be replace by the allele with the less number of reads. Otherwise the mixed allele is replaced based on a Bernoulli simulation as describe above.
  • When there is no read supporting one of the allele and the read count of the other allele is > 5, the mixed allele is replaced with the one that is not supported. Other it is replaced based on a Bernoulli simulation as describe above.

The process above is repeated nsim times and a new genotype file is generated and stored in the output directory. After all the iterations, the correlation between the MAF before phasing and the MAF from each phased genotype data is calculated. The phased genotype with the highest correlation will be retained and stored in the output SNPdata object. We recommend nsim = 50 at least to maximize the change of getting the highest MAF correlation.

snpdata <- phase_mixed_genotypes(snpdata, nsim = 100)
head(snpdata$Phased[, 1:5])

Impute the missing data

Some population genetics methods require genotype data with no missing allele. In {mpbr}, missing genotypes can be imputed using the impute_missing_genotypes() use the impute_missing_genotypes function to impute the missing genotypes. The genotype argument determines which genotype table should be used when imputing the missing data (“GT” to impute from raw data, “Phased” to impute from the phased data).

For every missing allele, we compute the reference and alternate allele frequencies at the given locus and identify the minor allele frequency (MAF). The missing allele will be replaced based on a Bernoulli simulation with a probability that is set at the resulting MAF value.

This process is repeated nsim times and a new genotype file is generated and stored in the output directory. After all the iterations, the correlation between the MAF before imputation and the MAF from each imputed genotype data is calculated. The imputed genotype with the highest correlation will be retained and stored in the output SNPdata object. We recommend nsim = 50 at least to maximize the change of getting the highest MAF correlation.

snpdata <- impute_missing_genotypes(snpdata, genotype = "GT", nsim = 100)
head(snpdata$Imputed[, 1:5])

Select data for a given chromosome

The package also contains functions that can be used to subset the SNPdata object. To generate a SNPdata object with data from chromosome 7 only (“Pf3D7_07_v3”) from an existing object of the same class, use the select_chrom() function. Note that the meta table in the output object will not have the percentage_missing_sites, Fws, COI columns. This is due to the fact that the existing values in these columns are calculated from the entire dataset and are not specific the chosen chromosome.

## subset only chromosome 7 data
chrom7_snpdata <- select_chrom(snpdata, chrom = "Pf3D7_07_v3")
print(chrom7_snpdata)

drop a set of SNPs from the data

If you want to remove some loci from the current SNPdata object, use the drop_snps() function and provide it with a data frame with 2 columns named as: Chrom and Pos”, representing the genomic coordinates of the loci to be dropped. Another alternative will be to provide the chrom (chromosome name), start (the start position of the region to be discarded) and end (the ens position of the region to be discarded).

## drop a set of SNPs based on their genomic coordinates
snp_to_be_dropped <- chrom7_snpdata$details %>% dplyr::select(Chrom, Pos)
idx               <- sample(seq_len(nrow(snp_to_be_dropped)), 100,
                            replace = FALSE)
snp_to_be_dropped <- snp_to_be_dropped[idx, ]

# there is a problem in the way the index of the vcf file is set when
# dropping snps.
# solution: consider replacing index with an appropriated value
reduced_snpdata   <- drop_snps(chrom7_snpdata, snp_to_be_dropped)
print(reduced_snpdata)

# drop a region from a chromosome
reduced_snpdata   <- drop_snps(chrom7_snpdata,
                               chrom = "Pf3D7_07_v3",
                               start = 110031, end = 154906)
print(reduced_snpdata)

drop a set of samples from the data

When you want to remove some samples from the SNPdata object, use the drop_samples() function and provide it with a vector of sample IDs to drop.

idx                   <- sample(seq_len(nrow(reduced_snpdata$meta)), 10,
                                replace = FALSE)
samples_to_be_dropped <- reduced_snpdata$meta$sample[idx]
reduced_snpdata       <- drop_samples(reduced_snpdata, samples_to_be_dropped)
print(reduced_snpdata)