GARCOM: A user-friendly R package for genetic mutation counts

Abstract

Next-generation sequencing (NGS) has enabled analysis of rare and uncommon variants in large study cohorts. A common strategy to overcome these low frequencies and/or small effect sizes relies on collapsing strategies, i.e. to bin variants within genes/regions. Several tools are now available for advanced statistical analyses; however, tools to perform basic tasks such as obtaining allelic counts within defined gene/region boundaries are unavailable or require complex coding. GARCOM (“Gene And Region Count Of Mutations”) library, an open-source freely available package in R language, returns a matrix with allelic counts within genes/regions per sample. GARCOM accepts input data in PLINK or VCF formats, with additional options to subset data for refined analyses.

Introduction

Genome-wide association studies (GWAS) have led to the identification of several genomic common variants associated with complex diseases, ¹ yet missing heritability remains extensive. Importantly, most of the disease-causing variants are rare in nature ² (minor allele frequency < 1%) where common variants serve as a proxy. Rapid decline in sequencing costs have enabled in-depth analysis of rare variants (RVs) through Whole-Genome sequencing (WGS) and Whole-Exome Sequencing (WES). Furthermore, large-scale reference panels have allowed for RV imputation. ^{3 - 5} Power to identify statistically significant RVs decreases as the minor allele frequency decreases: therefore, an ideal method to overcome this limitation is to group RV at the gene/region level, usually via a collapsing test.

Despite the availability of sophisticated tools for annotation, quality-control and association analyses, tools to perform basic tasks, for instance, obtaining allelic count within defined genetic boundaries (genes and/or regions) are lacking, to our knowledge. R libraries such as BEDMatrix and bigsnpr ⁶ provide allelic counts for each SNP per individual but algorithms to extract information within genetic boundaries in a collapsed fashion are unavailable.

Here we introduce a user-friendly R package, GARCOM (“ Gene And Region Count Of Mutations”) that computes allelic counts per individual within user-provided gene or genomic region boundaries.

Methods

GARCOM is written and developed in open-source R ⁷ statistical and programming language. GARCOM imports data.table, ⁸ vcfR, ⁹ bigstatr, bigsnpr and stats libraries for internal data transformation and processing. A stable version is released and publicly available on the CRAN repository.

install.packages(“GARCOM”)

Operation

GARCOM was developed on R (≥4.0) (RRID:SCR_017299) with other dependencies and minimum versions as: data.table (≥1.12.8), vcfR (≥1.12.0), bigsnpr (≥1.4.11). Full documentation of dependencies and installation is available at GARCOM GitHub repository. There is no minimum memory (RAM) requirement, but that may vary according to the nature and size of input genetics data. GARCOM was developed on Unix platform but can also be used on other platforms (e.g. Windows, Ubuntu).

Implementation

GARCOM operates through two main functions: “ gene_pos_counts” accepts PLINK ¹⁰ (RRID:SCR_001757) input data, whereas “ vcf_counts_SNP_genecoords” accepts VCF ¹¹ input format. After reading in the data, these functions perform operations to count variants within genes/genomic regions for each individual included in the file.

output <- gene_pos_counts(recoded_genetic_data, gene_boundaries, snp_locations)

output <- vcf_counts_SNP_genecoords(recoded_genetic_data, gene_boundaries, snp_locations)

where, “ output” is the object generated by GARCOM after a successful run of function; “ recoded_genetic_data” is the main input file in PLINK or VCF format; “ gene_boundaries”, and “ snp_locations” are additional input files for gene/regions and SNP information, respectively.

Typical workflow is shown in Figure 1. In brief, the “gene_pos_counts” function will process genetic input data (“ recoded_genetic_data”) generated from the PLINK software through the -- recode A option. Data are read in standard matrix format using the data.table R library. For VCF files, the “ vcf_counts_SNP_genecoords” function reads the VCF input file employing the extract.gt function from the vcfR library. The genotype values are read within the “GT” field.

Figure 1.

Workflow for standard GARCOM functions.

In addition to the -- recode A option, GARCOM requires genetic boundaries information and SNP information as shown in Table 2 and Table 3, respectively.

The output generated by GARCOM is a matrix, with M rows and N columns, where M represents the genes/genomic regions with at least one allele count and N represents the individuals. Genes/regions with zero allelic counts across all individuals are excluded from the final output. Missing values are counted as zero in final output. When no allelic counts are present in the user-defined genes, NULL value is returned. GARCOM allows missing values (NA) in input data.

In addition to the functions above described, GARCOM provides several options for user flexibility. For instance, GARCOM can be ran restricting analyses to 1) a list of genes, or 2) filter SNPs and extract individuals of interest. For instance, users can provide list of individuals using the “ keep_indiv” parameter; similarly, genes can be filtered in by using the “ filter_gene” parameter.

output <- gene_pos_counts(recoded_genetic_data, gene_boundaries, snp_locations, keep_indiv=mylist.txt)

output <- gene_pos_counts(recoded_genetic_data, gene_boundaries, snp_locations, filter_gene=mysetofgenes.txt)

Use cases

The input PLINK file has a matrix structure of N rows with M columns, where N rows represent individuals (one for each ID). The first six columns are family ID, individual ID, paternal ID, maternal ID, sex, and phenotype (standard output from PLINK ( Table 1)). The following columns consist of the variants included in the analyses.

Table 1.

Sample rows and columns for input genetics data recoded from PLINK software (--recode A).

FID	IID	PAT	MAT	SEX	PHENOTYPE	SNP1_A	SNP2_T	SNP3_G	SNP4_C	SNP5_C
FID1	IID_sample1	0	0	1	NA	1	1	0	NA	NA
FID2	IID_sample2	0	0	1	NA	0	1	0	NA	0
FID3	IID_sample3	0	0	1	1	0	0	1	0	0
FID4	IID_sample4	0	0	1	1	0	0	1	0	0
FID5	IID_sample5	0	0	1	1	0	0	1	0	0

Table 2.

Sample data for gene/region boundaries.

Data must contain “GENE”, “START” and “END” column names.

GENE	START	END
GENE1	100	180
GENE2	200	400

Table 3.

Sample data for SNP information where SNP and BP column names must be present in input data, where SNP is the single nucleotide polymorphism identifier and BP is base pair location.

SNP	BP
SNP1	100
SNP2	101
SNP3	201

Table 4.

Sample output, where GENE column identifies the gene/region names with corresponding allele counts per individual.

Individual_ID1, Individual_ID2 and Individual_ID3 are sample individual IDs, where values represent allelic counts within gene for individual.

GENE	Individual_ID1	Individual_ID2	Individual_ID3
GENE1	10	2	1
GENE2	2	1	0

The input VCF file follows the standard VCF formats (please refer to the vcfR library documentation).

Toy data (gene and SNP coordinates) are shared within the package as “ genecoord” and “ snppos”, respectively.

Simulation

We performed simulation on real data for chromosome 1 (“CHR1”, # of variants = 23,456) and CHR22 (# of variants = 4,814) on randomly sampled individuals (N = 100, 200, 500, 1000, 5000, 10,000) extracted from whole-exome sequencing dataset as described in the study by Tosto et al. ¹² Genetics data were recoded using PLINK -- recode A flag. For both chromosomes we found increased memory consumption and time ( Figure 2) as we increased the number of individuals included in the simulation. Memory consumption for CHR22 was significantly lower due to a smaller number of variants and genomic boundaries. Simulations were performed with 16GB memory (RAM) requested on computing cluster node.

Figure 2.

Comparison of memory (in MB) and CPU time (in seconds) for CHR1 and CHR22 on different sample sizes.

Graph A represents memory consumption; graph B shows processing time in seconds for various sample sizes on CHR1 and CHR22.

All simulations were conducted on R (v4.0), data.table (v1.13.6) with default 16 threads, GARCOM (v1.40), bigsnpr (v1.6.1).

Discussion

GARCOM is easy to use where basic knowledge of R programming language is desired. GARCOM is designed by harnessing existing libraries, such as data.table, that allow for efficient handling of large data. GARCOM data processing is independent of the reference genome build. GARCOM can be used on several platforms (e.g., Unix, Windows). GARCOM comes with certain limitations: genomic boundaries and variants' location need to be specified, as mentioned in the package documentation. GARCOM collapses variants based on base-pair (BP) location, hence processing multiple chromosomes at one time would add additional burden (time and memory) and can be hazardous because of identical BP across different chromosomes. In case of large-sized studies, e.g. UK Biobank (N ≥ 200K), processing data per chromosome is highly recommended due to above mentioned memory limitations. Lastly, GARCOM depends on public and freely available R packages.

Future

VCF format can accommodate locus annotation performed by software such as ANNOVAR. ¹³ To this end, GARCOM plans to accommodate annotation filters in addition to the existing ones. One challenge associated with annotated VCF is the resulting large file size; we will try to add this functionality, keeping RAM limitations and processing time in mind. We plan to add features to handle bgen format which stores large amount of genetics data, with appropriate R library ( http://www.well.ox.ac.uk/~gav/resources/rbgen_v1.1.5.tgz).

Data and software availability

Sample data associated with the package where applicable are provided within the library with proper documentation. No additional source data are required. We distribute the package under the MIT license. GARCOM can be downloaded from CRAN and GitHub from https://cran.r-project.org/web/packages/GARCOM/index.html and https://github.com/sariya/GARCOM respectively.

Reporting guidelines: Bugs and suggestions are welcome at the GitHub repository.

Author Contribution: SS, GT

Ethical Statement: Informed consent was obtained from all participants. For the whole-exome sequencing, the study protocol was approved by the Institutional Review Board (IRB) of Columbia university Medical Center (CUMC) (Approval number: AAAP0477). The study was conducted according to the principles expressed in the Declaration of Helsinki.