diff --git a/.github/workflows/c-cpp.yml b/.github/workflows/c-cpp.yml new file mode 100644 index 0000000..0a8d354 --- /dev/null +++ b/.github/workflows/c-cpp.yml @@ -0,0 +1,28 @@ +name: C/C++ CI + +on: + push: + branches: [ "main" ] + pull_request: + branches: [ "main" ] + workflow_dispatch: + +jobs: + build: + runs-on: ${{ matrix.config.os }} + name: ${{ matrix.config.os }} + strategy: + matrix: + os: [ ubuntu-latest, ubuntu-latest-arm, macos-15 ] + + steps: + - uses: actions/checkout@v4 +# - name: configure +# run: ./configure + - name: make + run: make + - name: make check + run: make check +# - name: make distcheck +# run: make distcheck + diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml new file mode 100644 index 0000000..c5185ae --- /dev/null +++ b/.github/workflows/main.yml @@ -0,0 +1,27 @@ +name: Main C/C++ CI + +on: + push: + branches: [ "main" ] + pull_request: + branches: [ "main" ] + workflow_dispatch: + +jobs: + build: + + runs-on: ubuntu-latest + + steps: + - name: clone repo + uses: actions/checkout@v4 + - name: make + run: make + - name: make check + run: make check + +# - name: make distcheck +# run: make distcheck +# - name: configure +# run: ./configure + diff --git a/.gitignore b/.gitignore index b881939..d3f64c9 100644 --- a/.gitignore +++ b/.gitignore @@ -1,6 +1,10 @@ build/ +tmp_* *.swp *.swo *~ *.DS_Store *.vscode/ +scratch/ +data/ +tags diff --git a/CLAUDE.md b/CLAUDE.md new file mode 100644 index 0000000..84133aa --- /dev/null +++ b/CLAUDE.md @@ -0,0 +1,40 @@ +# Purpose +- Compute the genetic similarity matrix from genetic data stored in the +vcf, vcf.gz, and bcf files. +- Use htslib C library to query / read genentic data. +- Write data to a binary file with a header and payload. + * header data must include all information that is + required to reproduce the calculation. + * payload consists of the computed values +- Easy to use command line interface + + +# External libraries +- htslib, code available on github at: https://github.com/samtools/htslib +- argparse, code available on github at: https://github.com/robert-vogel/argparse + + +# Code style +- variable, function, class, etc. names: + * should be descriptive (self-documenting) and short. + * classes and structs use Upper camel case + * functions and variables use snake case + - global constants use all upper case +- keep orthogonal services of the code in distinct modules +- module header files + * written to the 'include' directory. + * API documentation should be on line(s) proceeding the entities they + describe. +- implementation files in the 'src' directory, make sure to describe the +purpose / design of complex code blocks in the comments +- Unit tests should be written in the 'tests' directory and use Google +test and mock frameworks + + +# File and directory structure +- 'scratch' directory is not tracked by version control and is for testing +ideas by implementing in small programs +- 'build' directory is the target path for program and unit test builds +and should not be tracked by version control. +- 'include' directory for header files +- 'src' directory for 'main.cpp' file and module implementation files diff --git a/CMakeLists.txt b/CMakeLists.txt deleted file mode 100644 index 47f55ec..0000000 --- a/CMakeLists.txt +++ /dev/null @@ -1,109 +0,0 @@ - -# By: Robert Vogel -# Affiliation: Palmer Lab at UCSD -# 2025-01-11 -# -# Acknowledgment -# This file was originally written by Robert Vogel. Claude, -# the AI assistant by Anthropic reviewed this file and provided -# recommendations. - -cmake_minimum_required(VERSION 3.21) -project(hgrm VERSION 0.0.1) - -set(CMAKE_CXX_STANDARD 17) -set(CMAKE_CXX_STANDARD_REQUIRED True) -# set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fexceptions") -set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++17 -O3") - - -find_package(GTest REQUIRED) - -add_library(matrix_lib src/Matrix.cpp) -target_include_directories(matrix_lib PUBLIC include) - -add_library(utils_lib src/utils.cpp) -target_include_directories(utils_lib PUBLIC include) - -add_library(parse_lib src/HaplotypeDataRecord.cpp src/HaplotypeVcfParser.cpp) -target_include_directories(parse_lib PUBLIC include) - - - -# Testing configuration -enable_testing() - -add_executable( - test_matrix - tests/test_matrix.cpp -) - -target_link_libraries( - test_matrix - PRIVATE - matrix_lib - GTest::gtest_main -) - -add_executable( - test_utils - tests/test_utils.cpp -) -target_link_libraries( - test_utils - PRIVATE - utils_lib - GTest::gtest_main -) - - -add_executable( - test_haplotype_data_record - tests/test_haplotype_data_record.cpp -) -target_link_libraries( - test_haplotype_data_record - PRIVATE - parse_lib - matrix_lib - utils_lib - GTest::gtest_main -) - -add_executable( - test_haplotype_vcf_parser - tests/test_haplotype_vcf_parser.cpp -) - - - -target_link_libraries( - test_haplotype_vcf_parser - PRIVATE - parse_lib - matrix_lib - utils_lib - GTest::gtest_main -) - - -add_executable( - hgrm - src/main.cpp -) - -target_link_libraries( - hgrm - PRIVATE - parse_lib - matrix_lib - utils_lib -) - - -include(GoogleTest) -gtest_discover_tests(test_matrix) -gtest_discover_tests(test_utils) -gtest_discover_tests(test_haplotype_data_record) -gtest_discover_tests(test_haplotype_vcf_parser) - diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..e10d716 --- /dev/null +++ b/Makefile @@ -0,0 +1,194 @@ +# +# 2025 Palmer Lab +# +###################################################################### +# machine dependent options +###################################################################### +# CXXFLAGS note: Remember that -g flag is for generating source-level +# debug info. + +# OBJ_OUTPUT_OPTIONS: compiler options. Clang, and I presume +# also gcc, support the creation of dependency files (-MMD) and (-MP) phony +# targets required for constructing an object file. The (-o) and $@ +# are the standard output file designation submitted to the compiler, +# and $@ is an automatic variable storing the rules target. +# library archive program +# + +ifneq ($(shell which clang++),) +CXX = clang++ +CXXFLAGS = -pedantic # -Wextra +else ifneq ($(shell which g++),) +CXX = g++ +CXXFLAGS = -Wpedantic -Wextra +else +$(error "Couldn't establish either clang or gcc compiler availability") +endif + + +CXXFLAGS += -g -std=c++17 -Wall -Werror + +ifndef VIM +CXXFLAGS += -fdiagnostics-color=always +endif + +# Recall that -c flag prevents the compiler linking object files +OBJ_OUTPUT_OPTIONS = -c -MMD -MP -o $@ +AR = ar +AR_FLAGS = crs + +LOCAL_LIB = $(HOME)/.local/lib +LOCAL_LD = $(HOME)/.local/include + +###################################################################### +# define src and obj variables +###################################################################### + +SRC_DIR = src +HEADER_DIR = include +BUILD_DIR = build + +CXXLD := $(PWD)/include $(LOCAL_LD) $(CXXLD) +CXXLDFLAGS = $(addprefix -I, $(CXXLD)) + +CXXLIB += $(LOCAL_LIB) +CXXLIBFLAGS = $(addprefix -L, $(CXXLIB)) + +# APPLICATION FILES +# APP_FILES = matrix.cpp bcfio.cpp +APP_SRC = $(filter-out $(SRC_DIR)/main.cpp, $(wildcard $(SRC_DIR)/*.cpp)) +#APP_SRC = $(addprefix $(SRC_DIR)/, $(APP_FILES)) +APP_OBJS = $(subst $(SRC_DIR), $(BUILD_DIR), $(APP_SRC:.cpp=.o)) +APP_DEPS = $(APP_OBJS:.o=.d) + +# TESTS +TEST_DIR = tests +TEST_SRC = $(wildcard $(TEST_DIR)/test_*.cpp) +TEST_OBJS = $(subst $(TEST_DIR), $(BUILD_DIR), $(TEST_SRC:.cpp=.o)) +TEST_DEPS = $(TEST_OBJS:.o=.d) + + +TEST_DATA_SRC = $(wildcard $(TEST_DIR)/geno_test_data.*) +TEST_DATA_DST = $(subst $(TEST_DIR), $(BUILD_DIR), $(TEST_DATA_SRC)) + +TEST_TARGET_PRG = $(BUILD_DIR)/runtests + + +###################################################################### +# Executable Build Rules +###################################################################### + +TARGET = $(BUILD_DIR)/grm + +.PHONY: all +all: $(TARGET) $(TEST_TARGET_PRG) data + +$(TARGET): $(SRC_DIR)/main.cpp $(APP_OBJS) + $(CXX) $(CXXFLAGS) $(CXXLDFLAGS) $(CXXLIBFLAGS) -o $@ $^ -largparse -lhts + + +$(BUILD_DIR)/%.o: $(SRC_DIR)/%.cpp | $(BUILD_DIR) + $(CXX) $(CXXFLAGS) $(CXXLDFLAGS) $(OBJ_OUTPUT_OPTIONS) $< + +$(BUILD_DIR): + mkdir $@ + +###################################################################### +# Test Build Rules +###################################################################### + + +$(TEST_TARGET_PRG): $(TEST_DIR)/main.cpp $(TEST_OBJS) $(APP_OBJS) | $(TARGET) + $(CXX) $(CXXFLAGS) $(CXXLDFLAGS) $(CXXLIBFLAGS) -o $@ $^ -lgtest -lhts + +$(BUILD_DIR)/test_%.o: $(TEST_DIR)/test_%.cpp + $(CXX) $(CXXFLAGS) $(CXXLDFLAGS) $(OBJ_OUTPUT_OPTIONS) $< + +TEST_GRM_PRG = $(BUILD_DIR)/test_grm +test_grm: $(TEST_GRM_PRG) + ./$(TEST_GRM_PRG) + +$(TEST_GRM_PRG): $(BUILD_DIR)/test_grm.o $(BUILD_DIR)/grm.o | $(BUILD_DIR) + $(CXX) $(CXXFLAGS) $(CXXLDFLAGS) $(CXXLIBFLAGS) -o $@ $^ -lgtest -lgtest_main + +data: | $(TEST_DATA_DST) + +$(BUILD_DIR)/geno_test_data%: $(TEST_DIR)/geno_test_data% + rsync -avz $< $(BUILD_DIR)/ + +# tests: $(BUILD_DIR)/test_log #$(BUILD_DIR)/test_argparse +# +# $(BUILD_DIR)/test_log: $(BUILD_DIR)/test_log.o $(BUILD_DIR)/logger.o ~/.local/lib/libgtest.a +# $(CXX) $(CXXFLAGS) $(CXXLDFLAGS) $(CXXLIBFLAGS) -o $@ $^ +# +# $(BUILD_DIR)/test_argparse: $(BUILD_DIR)/test_argparse.o \ +# $(BUILD_DIR)/argparse.o \ +# ~/.local/lib/libgtest.a +# $(CXX) $(CXXFLAGS) $(CXXLDFLAGS) -I$(LOCAL_INCLUDE) -L$(LOCAL_LIB) -o $@ $^ + +# $(TEST_OBJS): $(TEST_SRC) +# $(CXX) $(CXXFLAGS) $(CXXLDFLAGS) -I$(LOCAL_INCLUDE) $(OBJ_OUTPUT_OPTIONS) $< +# +# $(BUILD_DIR)/test_log.o: $(TEST_DIR)/test_log.cpp +# $(CXX) $(CXXFLAGS) $(CXXLDFLAGS) $(OBJ_OUTPUT_OPTIONS) $^ + + +# $(TEST_OBJS): $(TEST_SRC) +# $(CXX) $(CXXFLAGS) $(CXXLDFLAGS) -I$(LOCAL_INCLUDE) -L$(LOCAL_LIB) -o $@ $^ + +# $(CXX) $(CXXFLAGS) $(CXXLDFLAGS) $(CXXLIBFLAGS) $(OBJ_OUTPUT_OPTIONS) $^ + + +###################################################################### +# +###################################################################### + +check: + ./$(TEST_TARGET_PRG) + +###################################################################### +# +###################################################################### + + +-include $(APP_DEPS) +-include $(TEST_DEPS) + +.PHONY: help +help: + -@echo "build grm" + -@echo "2025 Palmer Lab" + -@echo "" + -@echo "make grm executable" + -@echo "make libargparse" + + +###################################################################### +# install +###################################################################### + +# install: +# dir_header=$${prefix%/}/include/stitchr; \ +# if [ ! -d $${dir_header} ]; then \ +# mkdir -p $${dir_header}; \ +# fi; \ +# for hfile in $$(ls $(HEADER_DIR)); do \ +# cp $$hfile $${dir_header}/$${hfile}; \ +# done; \ +# \ +# dir_lib=$${prefix%/}/lib; \ +# if [ ! -d $${dir_lib} ]; then \ +# mkdir -p $${dir_lib}; \ +# fi; \ +# for libfile in $$(ls $(BUILD_DIR)/*.a); do +# cp $$libfile $${dir_lib}/$${libfile}; \ +# done; \ +# \ +# dir_bin = $${prefix%/}/bin; \ +# if [ ! -d $${dir_bin} ]; then \ +# mkdir -p $${dir_bin}; \ +# fi; \ +# cp $(TARGET) $${dir_bin}/$(notdir $(TARGET)) +# + + diff --git a/README.md b/README.md index c763477..205d532 100644 --- a/README.md +++ b/README.md @@ -1,79 +1,303 @@ -# 🏗️ Being built 🏗️ +# `hwas` a tool for haplotype wide association analyses -# Compute the genetic relationship matrix using expected haplotype counts + 🏗️ **Under construction** 🏗️ +## Table of Contents -The genetic relationship matrix (GRM) is the covariance between samples over -the measured genetic markers. It's utility is in accounting for -relatedness among samples, as random effects in a Linear Mixed Effects -Model, when performing a Genome Wide Association Study (GWAS) [1,2] and computing -the heritability of complex traits [3]. The GRM is traditionally computed using -SNP genotypes, but here we are interested in haplotype based covariance. +1. [About](#about) +2. [Genetic relationship matrices](#grm) +3. [Command line interface](#cli) +4. [Installation and requirements](#install) +5. [`.grm` file format](#grmspec) +4. [Contributing](#contributing) +4. [A.I. Acknowledgement](#ai) +5. [References](#refs) +## About -## Running the software +The genetic relationship matrix (GRM) describes the genetic relationship +between pairs of samples. Its computation is dependent on the random effects +defined by the linear mixed model mapping genetic features to phenotype. For +example, suppose that we are interested in accounting for polygenic SNP effects +using measured genotypes. Let the number of samples be $N$, the set of loci of +contributing to polygenic effects be $\Omega$ with $|\Omega|= M$, and the +quantitative phenotypes of $N$ samples $Y\in\mathbb{R}^{N\times 1}$. Note that +locus $j$ is not a member of the set of loci $\Omega$. Under the +LMM [[1]](#refs) -The program is ran by supplying the path and filename of a VCF. Note, as -of now this software does not support data streams with UNIX pipe or -bgzip, gzip, etc. compression. The output is printed to standard out. +$$ +Y = x_j\beta_j + \mathbf{Z} U + \epsilon.\\ +$$ + +with the fixed effect at locus $j$ being the alternative allele count, denoted +$x_j\in \{0,1,2\}^{N\times 1}$, and fixed effect size $\beta_j$. The +random polygenic and environmental effects have properties + +$$ +\begin{align} +U &\sim \mathcal{N}\left(0, \sigma_g^2 \\: \mathbf{I}_{M\times M}\right)\\ +\epsilon &\sim \mathcal{N}\left(0, \sigma_e^2 \\: \mathbf{I}_{N\times N}\right) +\end{align} +$$ + +with $\mathbf{Z}\in\{0,1,2\}^{N\times M}$ being the matrix of alt allele +counts of the $N$ samples and and $M$ loci in the set $\Omega$, +not a member. + +Under this model the sample phenotype covariance matrix +decomposes into genetic and environmental terms + +$$ +\text{cov}(Y) = \overbrace{ + \sigma^2_g \\; \mathbf{Z}\mathbf{Z}^T +}^{\text{genetic}} ++ +\underbrace{ + \sigma^2_e \\; \mathbf{I}_{N\times N} +}_{\text{environment}} +$$ + +where the genetic component of the phenotype covariance tells us how +to compute the GRM, i.e. $\mathbf{Z}\mathbf{Z}^T$. + +The alt allele count polygenic random effects are one example of genetic +effects. We do not need to limit ourselves to this model, and instead +account for any measurable genetic signals. This program uses several +signals to compute a GRM: + +* genotypes i.e. the alternative allele count, +* the expected alt allele count under a probabilistic model, useful for +imputed genetic signals, +* ancestral haplotypes, i.e. a K dimensional vector of expected haplotype +counts, +* or some combination of the aforementioned signals. + +This program provides an means to compute the GRM of genetic signals in +general. + + +## Genetic relationship matrices + +The genetic relationship matrices that this program computes are as follows: +imputed SNP genotypes considered above, the expected alternative allele +counts under a probabilistic model, the expected ancestral haplotype counts, +and the combined expected alternative allele and expected haplotype counts. +The subsections the follow define the model and the calculation for the +aforementioned GRMs. + +### SNP GRM + +The SNP GRM is presented in the [about](#about) section. Let $A_\text{SNP}$ be the +GRM computed by polygenic SNP effects, then + +$$ +A_\text{SNP} = \mathbf{Z}\mathbf{Z}^T +$$ + +with $\mathbf{Z}$ being the $N\times M$ matrix of alt allele counts. + + +### Expected alternative allele count GRM + +In many cases, as is the case in the Palmer Lab, SNPs are ***imputed***. If +the imputation method estimates the genotype probabilities at each locus of +each sample, as does our method of choice STITCH [[2]](#refs), then it may +be more appropriate consider the expected alternative allele counts (EAC) +under the imputation model rather than the imputed genotype calls. + +Let $Z_{im}\in\{0, 1, 2\}$ be an element of the design matrix of random +polygenic effects for sample $i$ at locus $m\in\Omega$. We consider $Z_{im}$ +to be random as the genotypes are imputed under a probabilistic model. This +makes the expected value $\mathbb{E}[Z_{im} |\mathcal{O}]$ of alternative +allele counts given the experimentally observed reads the an appropriate +genetic signal to consider. Let's call the matrix of expected alternative +allele counts $\mathbf{C} := \mathbb{E}[\mathbf{Z} |\mathcal{O}]$. Given +this, the GRM over expected alternative +allele counts $\mathbf{A}_\text{EAC}$ is + +$$ +\mathbf{A}_\text{EAC} =\mathbf{C}\mathbf{C}^T +$$ + +### Expected haplotype count GRM + +The expected haplotype count GRM $\mathbf{A}_\text{EHC}$ needs +motivation. The reason is that at each locus there is not a single +allele that we are counting, but instead we are counting the +number of copies of each haplotype $k\in\{1, 2, \dots, K\}$ at any +specified locus. Indeed, when $K=2$ we can cast the problem to +SNP case by identifying one of the two haplotypes as an alternative +allele. However when $K>2$ the genetic data is no longer a scalar +but a $K$ dimensional vector $h\in \\{0,1,2\\}^{K\times 1}$ such that +$\sum _{k=1}^K h_k = 2$. Moreover, this implies that at each locus $j$ +there are $K$ effect sizes, that can be expressed as the column +vector $\boldsymbol{\beta}_j\in \mathbb{R}^{K\times 1}$. As we can +see the haplotype model is more complex as the fixed effects are +in a $K$ dimensional space as opposed to a 1 dimensional space. + +An important question when working under the haplotype model is what +type of random genetic effects do we want to account for. If we +only care about the polygenic SNP effects, then we should make use of +GRMs $\mathbf{A}$ or $\mathbf{A}_\text{EHC}$. Another choice would +be to account for poly-haplotype effects, that the LMM mapping genotype +to phenotype at locus $j$ for sample $i$ becomes, + +$$ +\begin{equation} +Y_i = h_{ij}^T\\,\boldsymbol{\beta}_j + \mathbf{W}_1 U_1 + \mathbf{W}_2 U_2 + +\dots + \mathbf{W}_K U_K + \epsilon_i. +\end{equation} +$$ + +Here, the difference between the polygenic SNP and haplotype effects +are explicit. Instead of a single design matrix $\mathbf{Z}$ accounting +for polygenic effects there are now $K$ matrices. Consequently, the +haplotype GRM is + +$$ +\mathbf{A}_\text{EHC} = \sum_{k=1}^K \mathbf{W}_k\mathbf{W}^T_k +$$ + +the sum of the similarity matrices of each haplotype. + +### Leave one chromosome out (loco) GRM + +In a genome wide association study (GWAS) polygenic effects are often +accounted for from all loci except those of the chromosome that fixed +effect sizes are being estimated. This approach is referred to as +"leave-one-chromosome-out", or loco for short, for which we will refer +the resulting GRM as the loco GRM. This matrix can be computed easily +from the set of all chromosome GRMs. + +Let $M_c$ be the number of markers used for computing the GRM $A_c$. +Then it follows that the loco GRM $L_u$ that is used for effect size +estimation of all loci on chromosome $u$ is, + +$$ +\begin{align} +L_u &= \sum_{\forall c \neq u} A_c. +\end{align} +$$ + + +## Command line interface + +The `grm` program consists of two subprograms: + +* `grm contig`: the computation of the GRM of a named contig +* `grm loco`: the aggregation of contig GRMs into a + "leave-one-chromosome-out" matrix (denoted the "loco" matrix). + +to read documentation on the respective subprogram simply ``` -hgrm path/to/my_vcf > grm +grm [contig | loco] --help ``` +The GRM calculation requires tha the genetic information is in the +the bcf family of file formats, i.e. vcf, vcf.gz, or bcf. While +the loco subprogram requires all chromosome matrices to be in the +`.grm` binary file format defined below. + +### Computing a contig GRM + +By default, the GRM is computed +using the expected haplotype counts with FORMAT ID = "HD". + +``` +grm chrm +``` +will produce a binary `.mat` file that stores the GRM and relavent meta data. +Other options include + +### Computing a loco GRM + -## Installation and availability + + +## Installation and requirements The program is only available as source from this repository and requires -* `cmake` (>= 3.31.4) -* `make` +* `GNU make` +* `htslib` https://github.com/samtools/htslib +* `argparse` https://github.com/robert-vogel/argparse * `clang` or `gcc` C++17 compiler -To install navigate to the top level directory of this repository, -and make the `build` directory -```bash -mkdir build -``` -then use `cmake` to generate `make` files, etc., -```bash -cmake -S . -B build/ -``` -followed by navigating to the build directory and running GNU `make` -```bash -make -``` -In the build directory you should now have a binary file called `hgrm`, -this is the executable program. Move it to a directory in your -shell's search path. +## The `.grm` file format + +The `.grm` file format is a binary data format consisting of meta data +and a payload. + + +Assume 64-bit machine, little-endian, e.g. ARM and x86-64. + + +### Meta data + +| offset | type | size | description | +| (bytes) | | (bytes) | | +| --------- | --------- | --------- | ----------------------------------------- | +| 0 | uint32_t | 4 | File signature (0x47524D00 = "GRM\0") | +| 4 | uint32_t | 4 | File version, utils::Version | +| 8 | uint32_t | 4 | Program version, utils::Version | +| 12 | varies | $n_{coords}$ | Genomic coordinates, grm::Coordinates | +| $n_{coords}$ + 12 | varies | $n_{samps}$ | Sample names, grm::Samples | +| $n_{samps} + n_{coords} +12$ | | | Total | -## Contributing + +**Genomic coordinates** + +| offset | type | size | description | +| (bytes) | | (bytes) | | +| --------- | --------- | --------- | ----------------------------------------- | +| 0 | uint32_t | 4 | $l_{contig}=$ strlen(contig name) | +| 4 | char[] | $n_{contig} = l$ | Contig name chars, sizeof(char) = 1 byte | +| $n_{contig} + 4$ | uint32_t | 4 | Number of positions ($l_{pos}$) | +| $n_{contig} + 8$ | uint32_t | $n_{pos} = 4 l_{pos}$ | loci positions on contig | +| $n_{coords} =n_{contig} + n_{pos} + 8$ | | Total | + + +**Sample names** + +| offset | type | size | description | +| (bytes) | | (bytes) | | +| --------- | --------- | --------- | ----------------------------------------- | + +### Payload + +Given $N$ samples, the payload is the upper triangular components of the +$N \times N$ GRM stored as an array of $N(N+1)/2$ 32-bit floating point +numbers in row-major order. While the genotype-based GRM will not produce +fractional values, the expected count GRMs will, making `float32` an +acceptable choice for all GRM types. + + + +## Contributing I am using [GoogleTest](https://google.github.io/googletest/) framework for organizing tests. If you contribute, please make tests for your contributions. To run tests, `build` directory and build the project -``` -cmake -S ../ -B . -make -``` -Then use `cmake`'s utility -``` -ctest -``` +make check -## Acknowledgement -Code design and original version completed by Robert Vogel, -reviewed by Claude Sonnet, the AI assistant from Anthropic -(Jan 2025), with minor recommendations incorporated. -## References +## A.I. Acknowledgement -[1] [Kang et al. Genetics 178: 1709-1723 (2008)](https://academic.oup.com/genetics/article/178/3/1709/6061473) +The problem statement and overall design of the code base was by +Robert Vogel. He has made use of Claude for review and Claude Code, +the AI assistant by Anthropic, to implement sum features. -[2] [Kang et al. Nature Genetics 42 348-354 (2010)](https://www.nature.com/articles/ng.548) -[3] [Yang et al. Nature Genetics 42, 565-569 (2010)](https://www.nature.com/articles/ng.608) +## References +[1] [Yang et al. Nature Genetics 42, 565-569 (2010)](https://www.nature.com/articles/ng.608) +[2] [Davies et al. Nature Genetics 48, 965-969 (2016)](https://www.nature.com/articles/ng.3594) + diff --git a/include/HaplotypeVcfParser.h b/include/HaplotypeVcfParser.h deleted file mode 100644 index 511be81..0000000 --- a/include/HaplotypeVcfParser.h +++ /dev/null @@ -1,137 +0,0 @@ -// Parse STITCH vcf file -// -// -// -// By: Robert Vogel -// Affiliation: Palmer Lab at UCSD -// Date: 2025-01-09 -// -// -// Acknowledgment -// -// Code design and original version completed by Robert Vogel, -// reviewed by Claude Sonnet, the AI assistant from Anthropic -// (Jan 2025), with minor recommendations incorporated. -// -#ifndef HEADER_HAPLOTYPEVCFPARSER_H -#define HEADER_HAPLOTYPEVCFPARSER_H - -#include -#include -#include -#include -#include -#include -#include "Matrix.h" -#include "utils.h" - - -// samples are separated by white space -const char HAP_CODE[] { "HD" }; -const char META_PREFIX { '#' }; -const char MEASUREMENT_DELIM { ':' }; -const char HAP_DELIM { ',' }; -const int NUM_VCF_FIELDS { 9 }; -const char SPACE_DELIM { '\t' }; - - -// NOTE: in the future it may be best to test for set membership -static const char* VCF_FIELD_NAMES[NUM_VCF_FIELDS] { - "#CHROM", - "POS", - "ID", - "REF", - "ALT", - "QUAL", - "FILTER", - "INFO", - "FORMAT" -}; - - -// Move semantics, I don't want to copy data -class HaplotypeDataRecord -{ -public: - - HaplotypeDataRecord()=delete; - HaplotypeDataRecord(size_t, size_t); - HaplotypeDataRecord(const HaplotypeDataRecord&)=delete; - HaplotypeDataRecord(HaplotypeDataRecord&&)=delete; - - - const std::string& chrom() const; - const long pos() const; - const std::string& id() const; - const char ref() const; - const char alt() const; - const std::string& qual() const; - const std::string& filter() const; - const std::string& info() const; - const std::string& format() const; - - void parse_vcf_line(const char*); - const double& operator()(size_t, size_t) const; - - std::array dims() const; - - -private: - size_t n_samples_; - size_t k_founders_; - - std::string chrom_ { "" }; - long pos_ { -1 }; - std::string id_ { "" }; - char ref_ { '\0' }; - char alt_ { '\0' }; - std::string qual_ { "" }; - std::string filter_ { "" }; - std::string info_ { "" }; - std::string format_ { "" }; - - std::unique_ptr samples_ { nullptr }; - - StringRecord line_parse_ { SPACE_DELIM }; - StringRecord field_parse_ { MEASUREMENT_DELIM }; - StringRecord hap_parse_ { HAP_DELIM }; -}; - - -class HaplotypeVcfParser -{ -public: - - HaplotypeVcfParser()=delete; // default constructor - HaplotypeVcfParser(char* filename); // constructor - HaplotypeVcfParser(char* filename, size_t buffer_size); // constructor - //HaplotypeVcfParser(std::string filename); // constructor - HaplotypeVcfParser(const HaplotypeVcfParser&)=delete; // copy constructor - HaplotypeVcfParser(const HaplotypeVcfParser&&)=delete; // move constructor - HaplotypeVcfParser& operator=(const HaplotypeVcfParser&)=delete; // copy assignment - // ~HaplotypeVcfParser(); // descructor - - size_t n_samples() const; - size_t k_founders() const; - - bool load_record(HaplotypeDataRecord&); - -private: - const std::string fname_; - BufferedRead file_io_; - - CharBuffer line_buffer_; - size_t line_buffer_size_ { 0 }; - - size_t n_cols_ { 0 }; - size_t n_samples_ { 0 }; - size_t k_founders_ { 0 }; - size_t fpos_record_one_ { 0 }; - - - void pos_(size_t); - size_t get_line_num_char_(); - void set_params_(); -}; - -#endif diff --git a/include/Matrix.h b/include/Matrix.h deleted file mode 100644 index 821232e..0000000 --- a/include/Matrix.h +++ /dev/null @@ -1,46 +0,0 @@ -// -// By: Robert Vogel -// Affiliation: Palmer Lab at UCSD -// Date: 2025-01-09 -// -// -// Acknowledgment -// -// Code design and original version completed by Robert Vogel, -// reviewed by Claude Sonnet, the AI assistant from Anthropic -// (Jan 2025), with minor recommendations incorporated. -// -// -#ifndef HEADER_MATRIX_H -#define HEADER_MATRIX_H - -#include -#include -#include -#include -#include - -class Matrix -{ -public: - Matrix(size_t, size_t); // constructorconstructor - Matrix(const Matrix&); // copy constructor - Matrix(Matrix&&); // move constructor - Matrix& operator=(const Matrix&)=delete; // copy assignment - Matrix& operator=(Matrix&&)=delete; // move assignment - - - double operator()(const size_t&, const size_t&) const; - double& operator()(const size_t&, const size_t&); - - size_t size() const; - std::array dims() const; - -private: - const size_t nrow_; - const size_t mcol_; - std::unique_ptr data_; - size_t mat_idx_to_array_(const size_t&, const size_t&) const; -}; - -#endif diff --git a/include/bcfio.h b/include/bcfio.h new file mode 100644 index 0000000..9ab670c --- /dev/null +++ b/include/bcfio.h @@ -0,0 +1,250 @@ +// Parse STITCH vcf file +// +// By: Robert Vogel +// Affiliation: Palmer Lab at UCSD +// Date: 2025-01-09 +// +// +// Acknowledgment +// +// +#ifndef HEADER_PARSE_HTS_H +#define HEADER_PARSE_HTS_H + +#include +#include +#include +#include +#include + +#include + +namespace htslib { +extern "C" { +#include +#include +} +} + +// samples are separated by white space +// const char HAP_CODE[] { "HD" }; + +namespace bcfio { + + +// @title The meta data on a BCF attribute +// @description BCF, VCF, and VCF.GZ files hold metadata in the +// header that specify the type and format of data in records. +// I call each unique piece of data in a record a record attribute, +// e.g. an INFO column or FORMAT column of a record are attributes +// of that record. HTSLIB encodes attribute information in an +// unsigned 64 bit integer, and to access any value one needs to +// correctly implement bit shifting and masking. This struct contains +// bit-fields representing each value stored in the uint64_t. +// @bitfield number: the number of distinct values required to +// specify a sample +// record at loci i. For example, a SNP genotype is specified by a +// single +// string, e.g. 0/1, while the posterior genotype (0/0, 0/1, 1/1) +// probabilities requires three numbers. +// @bitfield vl_type: Specifies whether a variable is fixed length +// (BCF_VL_FIXED, in htslib/vcf.h line 68), variable length, etc. +// @bitfield type: the type of variable: binary flag (BCF_HT_FLAG), +// integer, real number, string, and 64 bit integers. Note that HT +// is header type. +// @bitfield coltype: +struct BcfHdrAttr { uint64_t number : 20, vl_type : 4, type : 4, coltype : 4; }; + + +// @title: Interface and manager of htslib bcf_hdr_t +// @description: The bcf header C-struct requires manual allocation +// and release of memory. This class manages applies RAII, reducing +// the chance of a memory leak. +class BcfHeader { +public: + BcfHeader() + : hdr_(nullptr) {}; + + BcfHeader(htslib::htsFile *fid) + : hdr_(fid ? htslib::bcf_hdr_read(fid) : nullptr) {}; + + ~BcfHeader() { if (hdr_) htslib::bcf_hdr_destroy(hdr_); }; + + bool isnull() const { return hdr_ == nullptr; }; + + + // @title: Retreive the set of smaple names + const std::unique_ptr sample_names() const; + + // @title: "get_*" member functions for attribute retrieval + // @description: + // @param id: the id of the formatted data field to retrieve + // @param ptr: the pointer to memory for which the BcfHdrAttr + // data will be copied into memory. + // @return 0 for success < 0 for fail + int get_format_attr(const char *id, BcfHdrAttr *ptr) const; + int get_info_attr(const char *id, BcfHdrAttr *ptr) const; + int get_filter_attr(const char *id, BcfHdrAttr *ptr) const; + + int subset_samples(const char *filename); + // @title: The number of values stored in format id + // @description: Each bcf format field is able to hold unique + // number of values per sample. This function provides a simple + // interface to the bcf file to retrieve that number. + // @param id: the format field id + // @return if an error occured that value returned is < 0, + // otherwise the number of values of fmt field id recorded per + // sample is returned. + int32_t k_fmt(const char *id) const; + + uint32_t n_samples() const { return hdr_->n[BCF_DT_SAMPLE]; }; + + // TODO: what unit test should I do for this? + const htslib::bcf_hdr_t *hts_hdr() const { return hdr_; }; + +private: + htslib::bcf_hdr_t *hdr_; + // BcfHdrAttr attr_ {}; + + // @title: + // @description decoder based upon htslib/vcf.h line 100 in the + // typedef struct bcf_idinfo_t. + // @param name: + // @param bcf_dt_type + // @param ptr + // @return -1 indicates an error has occured and 0 a success + int decode_hts_idinfo_(const char *name, + const int bcf_dt_type, + BcfHdrAttr *ptr) const; +}; + + +// @title: Interface and manage htslib bcf1_t +// @description: The htslib bcf1_t data structure requires manual memory +// management, knowledge of several bit-packed values, knowledge of +// several functions for querying data. This class simplifies +// memory management using C++ RAII idiom and provides a simplified, +// albeit non-comprehensive, interface for loading and querying data +// stored in the bcf1_t struct. +class BcfFloatRecord { +public: + + BcfFloatRecord(): rec_(htslib::bcf_init()) {}; + ~BcfFloatRecord(); + + // provide check-free fast, but unsafe, access to loaded data + float operator[](const size_t idx) const { return *(dst_ + idx); }; + + // provide index checked access to data. + std::optional get(const size_t row_idx, const size_t col_idx) const; + + // @title: Load sample data at the current locus + // @description: Sample data at the current locus, is not made + // available by reading a locus's record and storing in the + // bcf1_t type. Instead, we need to supply a pointer variable + // and format id to make that id's smaple data available. This + // function help simplify this process. + // @param hdr: instance of the bcf header to retreive meta data + // @param tag: the C-string id representing the data we want to + // query. + // @return 0 upon success and != 0 for failure + int load_data(BcfHeader *hdr, const char *tag); + + // the total amount of values per record, n_samples * k_founders + uint64_t size() const { return static_cast(ndst_); }; + uint64_t ncols() const { return col_num_; }; + uint64_t nrows() const { return row_num_; }; + + htslib::bcf1_t *cur_rec() const { return rec_; }; + + bool is_snp() const { return htslib::bcf_is_snp(rec_); } + +private: + htslib::bcf1_t *rec_; + + // These attributes store htslib access points to record data + // + // ndst_: The number of values in memory, the length of dst + // *dst_: an array of length ndst_ with float values of the + // current record + int ndst_ = 0; + float *dst_ = nullptr; + + // data that dst_ point to are stored in row major order, with + // columns being k_fmt and rows being n_samples. + uint64_t col_num_ = 0; + uint64_t row_num_ = 0; +}; + + +// @title Interface with htslib bcf +// @description ReadBCF manages the lifetime of an open htslib file +// and organizes the bcf file header and any one record for easy +// and memory safe parsing. +// @param bcfname: the path and filename to the bcf file to be read. +// @param sample_fname: the path and filename of the text file listing +// samples id's of records to be retreived. If this is not included +// all sample records are retrieved. +class ReadBcf +{ +public: + ReadBcf(); + ReadBcf(const char* filename, htslib::htsFile* fid); + + ReadBcf(const ReadBcf&)=delete; + ReadBcf& operator=(const ReadBcf&)=delete; + + ReadBcf(ReadBcf&&)=delete; + ReadBcf& operator=(ReadBcf&&)=delete; + + ~ReadBcf(); + + bool isopen() const; + + // @title: The number of values stored in format id + // @description: Each bcf format field is able to hold unique + // number of values per sample. This function provides a + // simple interface to the bcf file to retrieve that number. + // @param id: the format field id + // @return if an error occured that value returned is < 0, + // otherwise the number of values of fmt field id recorded per + // sample is returned. + int32_t k_fmt(const char *id) const { return hdr_.k_fmt(id); }; + + // See htslib/vcf.h line 649 + // Remember that n is the number of entries in the triplet of + // dictionaries in the VCF. BCF_DT_SAMPLE, provides the index of n + // that correspondes to the number of samples. + size_t n_samples() const { return hdr_.n_samples(); }; + + int set_samples(const char *filename); + + // TODO: sample_names + const std::unique_ptr sample_names() const { + return hdr_.sample_names(); + } + + int next_record(BcfFloatRecord *rec, const char *id); + +private: + const std::string fname_; + htslib::htsFile *fid_; + BcfHeader hdr_; + +}; + +// mode according to htslib: quoting from htslib/hts.h line 608 +// +// @example +// [rw]b .. compressed BCF, BAM, FAI +// [rw]bu .. uncompressed BCF +// [rw]z .. compressed VCF +// [rw] .. uncompressed VCF +// +// End quote +// +ReadBcf open(const char* filename, const char* mode); + +} + +#endif diff --git a/include/calc.h b/include/calc.h new file mode 100644 index 0000000..284e4b7 --- /dev/null +++ b/include/calc.h @@ -0,0 +1,20 @@ + +#ifndef HEADER_COV_CALC_H +#define HEADER_COV_CALC_H + +#include + +#include +#include +#include + +int compute_genotype_matrix(); + +int compute_eac_matrix(); + +// +int compute_ehc_matrix(Logger *log, bcfio::ReadBcf *bfid, Grm *cov); + +int compute_eac_and_ehc_matrix(); + +#endif diff --git a/include/constants.h b/include/constants.h new file mode 100644 index 0000000..2b8049e --- /dev/null +++ b/include/constants.h @@ -0,0 +1,14 @@ + +#ifndef HEADER_CONSTANTS_H +#define HEADER_CONSTANTS_H + +#include "utils.h" + + +namespace constants { + +constexpr utils::Version PROG_VERSION { 0, 0, 1 }; + +} + +#endif diff --git a/include/grm.h b/include/grm.h new file mode 100644 index 0000000..79d8866 --- /dev/null +++ b/include/grm.h @@ -0,0 +1,283 @@ +// Palmer Lab at UCSD +// +// This library provides the data structure of a genetic relationship matrix +// and functions for file I/O. +// +// GRM BINARY FILE SPECIFICATION +// +// A computed GRM is stored in a custom binary format. The extension ".grm" +// of these files is mandatory. The file is divided into two components, a +// header with meta-data necessary to reproduce the grm calculation and the +// the computed grm values, named the payload. +// +// The .grm file header is defined by the struct Hdr, and contains, at a +// minimum, the following information: +// * program_version: grm program version number +// * data_type: alt_count, expected_alt_count, expected_haplotype_count, +// both expected_alt_count and expected_haplotype_count. +// * coords: Genomic coordinates used in the grm calculation. +// * samples: list of sample id's in order of the grm +// the coords and samples are defined by their own structs with field pointers +// to heap allocated memory addresses. Reading and writing such heap allocated +// structs make use of runtime polymorphism of function "read" and "write" +// +// +// ACKNOWLEDGMENT +// +// Code design and original version completed by Robert Vogel, +// reviewed by Claude Opus 4.6, the AI assistant from Anthropic. +// Some recommendations have been incorporated. +// +#ifndef HEADER_GRM_H +#define HEADER_GRM_H + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "io.h" +#include "constants.h" +#include "utils.h" + + + +// The algorithm for getting the array idx from matrix indexes is simply +// +// idx = i * n_samples - n_skipped_idxs + j +// +// where i is the matrix row index and j is the matrix column index. The +// interesting term is n_skipped_idxs, this is the number of elements +// that referencing (i, j) skip when only storing upper triangle. For +// example, suppose that i = 3 and j = 2. Here three complete rows of +// the matrix has been traversed, therefore the number skipped is +// +// n_skipped_idx >= (i-1) * i / 2 +// +// if j < i, then as we are only storing the upper triangle and +// by symmetry swap the values of i an j. If j > i, then the number +// skipped in row i is i, making the total skipped, +// +// n_skipped_idx = i * (i - 1) / 2 + i +// +// making the equation above read +// +// idx = i * n_samples - i * (i - 1)/2 - i + j +// +// Then let's manually validate +// +// n = 3 +// i j idx_true idx +// 0 0 0 0 +// 0 1 1 1 +// 0 2 2 2 +// 1 0 --- transpose -> 0, 1 --- +// 1 1 3 1*3 - 1*(0)/2 - 1 + 1 = 3 +// 1 2 4 1*3 - 1*(0)/2 - 1 + 2 = 4 +// 2 2 5 2*3 - 2*(1)/2 - 2 + 2 = 6 - 1 = 5 +// +#define MATRIX_IDX_TO_ARRAY(i, j, n) ((i)*(n) - (i)*((i)-1)/2 + (j) - (i)) + + +namespace grm { + + +// Recall that the magic number is simply GRM\0 in hex, +// 0x47 = G, 0x52 = R, and so on. The magic number provides +// a simple means to determine file type when parsing. +constexpr uint32_t FILE_TYPE_SPEC = 0x47524D00; +constexpr utils::Version FILE_VERSION = { 0, 0, 0 }; +constexpr char FILE_SUFFIX[] = ".grm"; + + +enum STATUS { + SUCCESS, + FAILED, + UNKNOWN_FAILURE, + ERROR_IDX_ARR_BOUNDS, + ERROR_FOPEN, + ERROR_EOF_NOT_REACHED, + ERROR_ON_WRITE, + ERROR_ON_READ, + ERROR_FILE_NOT_OPEN, + ERROR_NULLPTR_ARG, + ERROR_INVALID_ARG, + ERROR_NOT_A_GRM_FILE, +}; + + +enum GrmType { + EHC, // Expected Haplotype Count + EAC, // Expected Alternative Allele Count + BOTH, // Both EHC AND EAC + DS, // Dosage, i.e. Called Alternative Allele Count + UNSPECIFIED, +}; + +// @title: Store genomic coordinates used in GRM calculation +struct Coordinates { + Coordinates(): contig(""), len(0), pos(nullptr) {}; + Coordinates(char* contig_in, uint64_t len_in) + : contig(contig_in == nullptr ? "" : contig_in), + len(contig == "" ? 0 : len_in), + pos(len == 0 ? nullptr : std::make_unique(len)) {}; + + Coordinates(const Coordinates&) = delete; + Coordinates& operator=(const Coordinates&) = delete; + + Coordinates(Coordinates&& other); + Coordinates& operator=(Coordinates&& other); + + // Data Fields + std::string contig; + uint64_t len; + std::unique_ptr pos; +}; + + +// Coordinates Storage Layout +// +// type number description +// -------------------------------------------------------------------- +// size_t 1 number of characters (n) in contig name +// char n characters for contig name without null character +// size_t 1 number of genomic positions (npos) +// size_t npos the positions used for computation of the grm +// +STATUS write(io::FileIO* fio, const Coordinates* coords); +STATUS read(io::FileIO* fio, Coordinates* coords); + + +// Samples stores sample id strings and the number of samples +// +struct Samples { + Samples(): len(0), names(nullptr) {}; + Samples(uint64_t n_samples): + len(n_samples), + names(len == 0 ? nullptr : std::make_unique(len)) {}; + + Samples(const Samples&) = delete; + Samples& operator=(const Samples&) = delete; + + Samples(Samples&& other); + Samples& operator=(Samples&& other); + + // Data Fields + uint64_t len; //number of samples + std::unique_ptr names; + +}; + + +// Sample Storage Layout +// +// type number description +// -------------------------------------------------------------------- +// size_t 1 represents number of samples +// size_t 1 the number of characters of longest sample id +// size_t 1 number of characters (n_1) in first sample id +// char n_1 characters of sample id 1 without terminal null '\0' +// size_t 1 number of characters (n_2) in second sample id +// char n_2 characters of sample id 2 without terminal null '\0' +// ... +// size_t 1 number of characters (n_N) in N^{th} sample id +// char n_N characters of sample id N without terminal null '\0' + +STATUS write(io::FileIO* fio, const Samples* samples); +STATUS read(io::FileIO* fio, Samples* samples); + + +// Header +struct Hdr { + + Hdr(); + Hdr(const Hdr&) = delete; + Hdr& operator=(const Hdr&) = delete; + + Hdr(Hdr&&); + Hdr& operator=(Hdr&&); + + // Data Fields + utils::Version prog_version; + utils::Version file_version; + + GrmType grm_type; + std::unique_ptr coords; + std::unique_ptr samples; + +}; + +// Header Storage Layout +// +// type number description +// -------------------------------------------------------------------- +// size_t 1 number of characters (n) in version string +// char n version string without null terminator +// GrmType 1 type of grm +// +// call write for coordinates +// +// call write for samples +// +STATUS write(io::FileIO *fio, const Hdr *hdr); +STATUS read(io::FileIO *fio, Hdr *hdr); + + +// Grm class manages storage and access of GRM matrix +// +// The GRM as an n_sample by n_sample symmetric, positive semi-definite +// matrix. Let Z represent the n_sample by m_marker data genetic data. +// From these data the GRM is computed as GRM = ZZ^T. +// +// @param n_samples of the GRM. +// +struct Grm { + // + Grm(); + Grm(uint64_t n_samps); + + Grm(const Grm&)=delete; + Grm& operator=(const Grm&)=delete; + + Grm(Grm&&); + Grm& operator=(Grm&&); + + // Unchecked indexes when setting and getting of matrix values + float operator()(const uint64_t i, const uint64_t j) const; + float& operator()(const uint64_t i, const uint64_t j); + + // Checked indexes when setting and getting of matrix values + STATUS set(const uint64_t i, const uint64_t j, const float val); + STATUS get(const uint64_t i, const uint64_t j, float *val) const; + + uint64_t size() const; + + STATUS midx_to_arr(const uint64_t i, + const uint64_t j, uint64_t* idx) const; + + uint64_t n_samples; + std::unique_ptr data; +}; + +// @title: Write meta-data and computed grm elements to file +// @description: The binary file written contains a header and payload: +// * Header +// - an instance of grm::Header +// * Payload +// - grm data in row major order +// @param filename: name of file that the data are written +// @param hdr: an instance of grm::Hdr with important meta data +// @return grm::STATUS: +// +STATUS write(io::FileIO* fio, const Hdr *hdr, const Grm *grmatrix); +STATUS read(io::FileIO* fio, Hdr* hdr, Grm* grmatrix); + + +} + +#endif diff --git a/include/io.h b/include/io.h new file mode 100644 index 0000000..0d1b7fd --- /dev/null +++ b/include/io.h @@ -0,0 +1,102 @@ +#ifndef HEADER_IO_H +#define HEADER_IO_H + +#include +#include +#include + +namespace io { + +enum STATUS { + SUCCESS, + FERROR, + FEOF, + INVALID_ARG_ERROR, + FSEEK_ERROR, + FEOF_ERROR, + END_OF_BUF_ERROR +}; + + +struct FileIO { + FileIO(FILE* fid): fid(fid) {}; + ~FileIO() { if (fid) { fclose(fid); fid = nullptr; } }; + + FileIO(const FileIO&) = delete; + FileIO& operator=(const FileIO&) = delete; + + FileIO(FileIO&& other) noexcept : fid(other.fid) { + other.fid = nullptr; + } + + FileIO& operator=(FileIO&& other) { + // protect against self assignment + if (this == &other) return *this; + + if (fid) fclose(fid); + + fid = other.fid; + other.fid = nullptr; + + return *this; + } + + FILE *fid; +}; + + +// @title: open a file +// @description: +// @param filename: name and path of file to open +// @param mode: a mode in the set of those in the C library function fopen +// @return a pointer to opened file +inline FileIO open(const char *filename, const char *mode) { + if (!mode || !filename) + return nullptr; + + FILE *fid = fopen(filename, mode); + if (!fid) + return nullptr; + + // compiler implements elision + return FileIO(fid); +} + + +// @title: file object +// @description: This class manages the lifetime of a C-style file stream +// by RAII. To contruct an instance of the class use the "open" function +// below. +// @param fid: an opened C-style file stream +// int bseek(FileIO *fio); + + + +// @title: File statistics +// @description: This object is returned by any function meant to calculate +// file character statistics. +// struct FileStats { +// size_t nchar = 0; +// size_t nwords = 0; +// size_t nlines = 0; +// size_t nblanklines = 0; +// } + + +// @title: word count +// @description: Similar to the UNIX/Linux wc command line program, wc +// calculates the number of characters, words, lines, etc. that +// the specified file contains. +// @param tio: an instance of TextIO +// @param fs: the structure that the file statistics will be stored +// @return a STATUS code that specifies whether the function was successful +// or failed. +// STATUS wc(TextIO *tio, FileStats *fs); + + + + +// STATUS getline(TextIO *tio, Array linebuf); +} + +#endif diff --git a/include/logger.h b/include/logger.h new file mode 100644 index 0000000..49504a5 --- /dev/null +++ b/include/logger.h @@ -0,0 +1,50 @@ + +#ifndef HEADER_LOGGER_H +#define HEADER_LOGGER_H + +#include +#include +#include +#include +#include + + +class Logger { +public: + Logger(); + + // TODO: right now only accepts a single msg string, I should make + // this arbitrary message elements using va_list, this makes the + // interface match that of sprintf + int info(const char *format, ...); + int warn(const char *format, ...); + int error(const char *format, ...); + + // int info(const char *msg); + // int warn(const char *msg); + // int error(const char *msg); + +private: + time_t t_; + tm *time_point_; + + int status_ { 0 }; + + static constexpr size_t time_buf_len_ { 30 }; + static constexpr size_t str_buf_len_ { 500 }; + static constexpr size_t max_str_ { 450 }; + + char time_buf_[time_buf_len_]; + char str_buf_[str_buf_len_]; + + static constexpr char err_str_[] = { "ERROR" }; + static constexpr char warn_str_[] = { "WARN" }; + static constexpr char info_str_[] = { "INFO" }; + + int load_time_buf_(); + void empty_bufs_(); + void vprintf_(FILE *stream, const char *log_type, + const char *format, va_list arg_ptr); +}; + +#endif diff --git a/include/utils.h b/include/utils.h index f9830b5..adf574d 100644 --- a/include/utils.h +++ b/include/utils.h @@ -1,102 +1,62 @@ #ifndef HEADER_UTILS_H #define HEADER_UTILS_H - -#include -#include -#include -#include -#include -#include -#include -#include - - - -class CharBuffer { -public: - CharBuffer(); - CharBuffer(size_t); - CharBuffer(const CharBuffer&)=delete; - CharBuffer(CharBuffer&&)=delete; - - const char& operator()(size_t) const; - const size_t& size() const; - const size_t& buffer_size() const; - - void append(char s); - void reset(); // set buffer_idx_ to zero - void reset(size_t buffer_size); // set buffer_idx_ to zero - - // remember that it is the caller's responsibility to not - // dereference raw pointer after the CharBuffer instance is - // destructed - const char* data() const; - -private: - size_t buffer_size_; - size_t buffer_idx_; - std::unique_ptr buffer_; -}; - - - -class StringRecord { -public: - // input only delimiters - StringRecord()=delete; - StringRecord(const char); - StringRecord(const char, const size_t); - StringRecord(const char, const char*); - - void update_str(const char* s); - - const char* data() const; - void reset(); - size_t size(); - bool next_field(); - - -private: - size_t size_ { 0 }; - size_t idx_ { 0 }; - const char* str_; - const char delim_; - - CharBuffer buf_; - std::function is_delim_; - bool char_is_delim_(char) const; - +#include + +namespace utils { + +// Claude recommends that I define serialization explicitly so +// that the bit packed field order is not ABI dependent +struct Version { + uint8_t major; + uint16_t minor; + uint16_t micro; + + // pack as 8, 12, 12 for 32 bits in total + uint32_t pack() const { + return (static_cast(major) << 24 + | static_cast(minor & 0x0FFF) << 12 + | static_cast(micro & 0x0FFF)); + } + + static Version unpack(uint32_t vnum) { + return Version { + static_cast(vnum >> 24), + static_cast(vnum >> 12 & 0x0FFF), + static_cast(vnum & 0x0FFF) + }; + } }; +// template +// struct Array { +// Array(size_t size_in): size(size_in), +// data(size > 0 ? new T[size] : nullptr) {}; +// +// ~Array() { if (data) delete[] data; }; +// +// size_t size; +// T *data; +// size_t len = 0; +// +// //unsafe referencing +// T operator[](size_t i) { return data[i]; }; +// T& operator[](size_t i) { return data[i]; }; +// +// STATUS append(T val) { +// if (len >= size-1) +// return END_OF_BUF_ERROR; +// +// data[len++] = val; +// return SUCCESS; +// } +// +// void fill(T val) { +// std::memset(data, val, size); +// len = 0; +// } +// }; + +} -class BufferedRead { -public: - BufferedRead()=delete; - BufferedRead(const BufferedRead&)=delete; - BufferedRead(BufferedRead&&)=delete; - BufferedRead& operator=(const BufferedRead&)=delete; - - BufferedRead(char* filename, size_t buff_size); - ~BufferedRead(); - - size_t get_line(CharBuffer& line_buf); - // size_t get_line(std::unique_ptr line_buf); - char get_char(); - void seek(size_t n); - size_t tell(); - void reset(); - -private: - const char* filename_; - const size_t buff_size_; - - FILE* fid_; - std::unique_ptr buffer_; - - size_t buffer_pos_ { 0 }; - - size_t update_buffer_(); - -}; #endif diff --git a/src/HaplotypeDataRecord.cpp b/src/HaplotypeDataRecord.cpp deleted file mode 100644 index 27ab865..0000000 --- a/src/HaplotypeDataRecord.cpp +++ /dev/null @@ -1,139 +0,0 @@ -// -// -// By: Robert Vogel -// Affiliation: Palmer Lab at UCSD -// Date: 2025-01-09 - -// -// -// Acknowledgment -// -// Code design and original version completed by Robert Vogel, -// reviewed by Claude Sonnet, the AI assistant from Anthropic -// (Jan 2025), with minor recommendations incorporated. -// -#include "HaplotypeVcfParser.h" - - -// Default constructor -HaplotypeDataRecord::HaplotypeDataRecord(size_t n_samples, size_t k_founders) - : n_samples_(n_samples), - k_founders_(k_founders), - samples_(n_samples_ > 0 && k_founders_ > 0 - ? std::make_unique(n_samples_, k_founders_) : nullptr) { - - if (n_samples_ <= 0 || k_founders_ <= 0) - throw std::runtime_error("Data must have more than zero samples and founders"); - - }; - - -// access elements -const std::string& HaplotypeDataRecord::chrom() const { return chrom_; }; -const long HaplotypeDataRecord::pos() const { return pos_; }; -const std::string& HaplotypeDataRecord::id() const { return id_; }; -const char HaplotypeDataRecord::ref() const { return ref_; }; -const char HaplotypeDataRecord::alt() const { return alt_; }; -const std::string& HaplotypeDataRecord::qual() const { return qual_; }; -const std::string& HaplotypeDataRecord::filter() const { return filter_; }; -const std::string& HaplotypeDataRecord::info() const { return info_; }; -const std::string& HaplotypeDataRecord::format() const { return format_; }; - - -void HaplotypeDataRecord::parse_vcf_line(const char* vcf_line) { - - - if (std::isspace(vcf_line[0])) - throw std::runtime_error("No line can begin with spaces"); - - - // A field of a vcf line record is a single string separated from others - // by white space. - // A sample field is a field with the data for a single sample. Sample - // fields have numerous : delimited records - // The counts of the k founders in any one sample field is a comma delimited - // element of a sample field record. - size_t hap_idx { 0 }; // index with hap counts - bool hap_found { false }; // determine whether hap dose is in dataset - size_t sample_idx { 0 }; // sample index - size_t founder_idx { 0 }; - - line_parse_.update_str(vcf_line); - - for (int field_idx = 1; line_parse_.next_field(); field_idx++) { - - if (field_idx == 1) - chrom_ = line_parse_.data(); - else if (field_idx == 2) - pos_ = std::atoi(line_parse_.data()); - else if (field_idx == 3) - id_ = line_parse_.data(); - else if (field_idx == 4) - ref_ = line_parse_.data()[0]; - else if (field_idx == 5) - alt_ = line_parse_.data()[0]; - else if (field_idx == 6) - qual_ = line_parse_.data(); - else if (field_idx == 7) - filter_ = line_parse_.data(); - else if (field_idx == 8) - info_ = line_parse_.data(); - else if (field_idx == 9) { - format_ = line_parse_.data(); - - // verify in the format field that haplotype dose (HD) - // is included in the data. Find the index (hap_idx) - // for which haplotype count data is found in a sample field - // record - field_parse_.update_str(line_parse_.data()); - - hap_found = false; - for (size_t i = 0; field_parse_.next_field(); i++) { - if (std::strcmp(field_parse_.data(), HAP_CODE) == 0 ) { - hap_found = true; - hap_idx = i; - break; - } - hap_idx++; - } - - if (!hap_found) - throw std::runtime_error("Haplotype counts are not specified"); - - } else if (field_idx > 9 && samples_) { - - if (sample_idx < 0 || sample_idx >= n_samples_) - throw std::out_of_range("Index is out of matrix range."); - - field_parse_.update_str(line_parse_.data()); - for (size_t j = 0; field_parse_.next_field(); j++) - if (j == hap_idx) - break; - - hap_parse_.update_str(field_parse_.data()); - - // decompose haplotype counts to respective founders - for (founder_idx = 0; hap_parse_.next_field(); founder_idx++) - (*samples_)(sample_idx, founder_idx) = std::atof(hap_parse_.data()); - - - if (founder_idx != k_founders_) - throw std::runtime_error("Number of founders found for sample is incorrect"); - - sample_idx++; - } - } - - if (sample_idx != n_samples_) - throw std::runtime_error("Number of samples found is not equal to that expected."); - -} - - -const double& HaplotypeDataRecord::operator()(size_t i, size_t j) const { - return (*samples_)(i, j); -} - -std::array HaplotypeDataRecord::dims() const { - return (*samples_).dims(); -} diff --git a/src/HaplotypeVcfParser.cpp b/src/HaplotypeVcfParser.cpp deleted file mode 100644 index affda4a..0000000 --- a/src/HaplotypeVcfParser.cpp +++ /dev/null @@ -1,220 +0,0 @@ -// -// By: Robert Vogel -// Affiliation: Palmer Lab at UCSD -// Date: 2025-01-09 -// -// Input argument -// filename: vcf with haplotpye -// -// -// Acknowledgment -// -// Code design and original version completed by Robert Vogel, -// reviewed by Claude Sonnnet, the AI assistant from Anthropic -// (Jan 2025), with minor recommendations incorporated. - -#include "HaplotypeVcfParser.h" - -const static size_t DEFAULT_BUFFER_SIZE { 100000 }; - - - - -HaplotypeVcfParser::HaplotypeVcfParser(char* filename) - : fname_(filename), - file_io_(BufferedRead(filename, DEFAULT_BUFFER_SIZE)) { - - // get number of characters in data record for line buffer size - size_t nchar { get_line_num_char_() }; - - if (nchar == 0) - throw std::runtime_error("No data to read"); - - // make buffer 10% larger then the number of characters read. - line_buffer_size_ = static_cast(nchar * 1.1); - line_buffer_.reset(line_buffer_size_); - - set_params_(); - pos_(fpos_record_one_); -}; - - -HaplotypeVcfParser::HaplotypeVcfParser(char* filename, size_t buff_size) - : fname_(filename), - file_io_(BufferedRead(filename, buff_size)) { - - // get number of characters in data record for line buffer size - size_t nchar { get_line_num_char_() }; - - if (nchar == 0) - throw std::runtime_error("No data to read"); - - // make buffer 10% larger then the number of characters read. - line_buffer_size_ = static_cast(nchar * 1.1); - line_buffer_.reset(line_buffer_size_); - - set_params_(); - pos_(fpos_record_one_); -}; - - -// HaplotypeVcfParser::HaplotypeVcfParser(std::string filename) -// : fname_(filename), -// fid_(filename) { -// -// if (fid_.bad()) -// throw std::runtime_error("File Access error"); -// else if (fid_.eof()) -// throw std::runtime_error("File is empty"); -// -// // get number of characters in data record for line buffer size -// size_t nchar { get_line_num_char_() }; -// -// -// // make buffer 10% larger then the number of characters read. -// line_buffer_size_ = static_cast(nchar * 1.1); -// line_buffer_ = new char[line_buffer_size_]; -// line_buffer_[0] = '\0'; -// -// pos_(std::ios_base::beg); -// set_params_(); -// }; - - -//HaplotypeVcfParser::~HaplotypeVcfParser() { -// if(fid_) -// fclose(fid_); -// // free(line_buffer_); -//} - - -size_t HaplotypeVcfParser::get_line_num_char_() { - - size_t char_count { 0 }; - size_t max_char_count { 0 }; - char c; - - while ((c = file_io_.get_char()) != '\0') { - char_count++; - - if (c == '\n' && char_count > max_char_count) { - max_char_count = char_count; - char_count = 0; - } - } - - return max_char_count; -} - - -size_t HaplotypeVcfParser::n_samples() const { return n_samples_; } - - -size_t HaplotypeVcfParser::k_founders() const { return k_founders_; } - - -void HaplotypeVcfParser::pos_(size_t n) { - file_io_.reset(); - if (n != 0) - file_io_.seek(n); -} - - -void HaplotypeVcfParser::set_params_() { - pos_(0); - - // skip meta data lines - - size_t n { 0 }; - size_t num_bytes { 0 }; - while ((n = file_io_.get_line(line_buffer_)) > 0) { - - // the +1 is because I don't write newline characters to buffer - num_bytes += line_buffer_.size()+1; - - if (line_buffer_(0) == META_PREFIX && line_buffer_(1) == META_PREFIX) - continue; - - break; - } - - // Store file position of first record - fpos_record_one_ = sizeof(line_buffer_(0)) * num_bytes; - - // if there is no header - if (line_buffer_(0) != META_PREFIX) - return; - - - if (std::isspace(line_buffer_(0))) - throw std::runtime_error("First element of VCF line must not be blank."); - - - // Get column number and sample number - StringRecord line_parser_ { SPACE_DELIM, line_buffer_.data() }; - StringRecord field_parser_ { MEASUREMENT_DELIM }; - StringRecord hap_parser_ { HAP_DELIM }; - - n_cols_ = 0; - n_samples_ = 0; - for (; line_parser_.next_field(); n_cols_++) { - - if (n_cols_ < NUM_VCF_FIELDS - && std::strcmp(line_parser_.data(), VCF_FIELD_NAMES[n_cols_]) != 0) - throw std::runtime_error("File doesn't follow vcf header specification"); - - if (n_cols_ >= NUM_VCF_FIELDS) - n_samples_++; - - } - - - // get k founders from record - if ((n = file_io_.get_line(line_buffer_)) == 0) - throw std::runtime_error("End of file"); - - line_parser_.update_str(line_buffer_.data()); - size_t hap_idx { 0 }; - for (int i = 0; line_parser_.next_field(); i++) { - - if (i == NUM_VCF_FIELDS-1) { - field_parser_.update_str(line_parser_.data()); - - for (;field_parser_.next_field(); hap_idx++) - if (std::strcmp(field_parser_.data(), HAP_CODE) == 0) - break; - - } else if(i == NUM_VCF_FIELDS) { - field_parser_.update_str(line_parser_.data()); - - for(int i = 0; field_parser_.next_field() && i < hap_idx; hap_idx++) - ; - - hap_parser_.update_str(field_parser_.data()); - for (;hap_parser_.next_field(); k_founders_++) - ; - - break; - } - - } - - if (k_founders_ == 0) - throw std::runtime_error("Parse error"); - -} - - -bool HaplotypeVcfParser::load_record(HaplotypeDataRecord& record) { - - size_t n { 0 }; - - if ((n = file_io_.get_line(line_buffer_)) == 0) - return false; - - record.parse_vcf_line(line_buffer_.data()); - - return true; -} - - diff --git a/src/Matrix.cpp b/src/Matrix.cpp deleted file mode 100644 index 32d7df0..0000000 --- a/src/Matrix.cpp +++ /dev/null @@ -1,69 +0,0 @@ -// MAtrix -// -// By: Robert Vogel -// Affiliation: Palmer Lab at UCSD -// Date: 2025-01-10 -// -// -// Acknowledgment -// -// Code design and original version completed by Robert Vogel, -// reviewed by Claude Sonnet, the AI assistant from Anthropic -// (Jan 2025), with minor recommendations incorporated. -// -// - -#include "Matrix.h" - -// default constructor -Matrix::Matrix(size_t nrow, size_t mcol) - : nrow_(nrow), mcol_(mcol), - data_(nrow_ > 0 && mcol_ > 0 ? std::make_unique(size()) : nullptr) { - - if (nrow_ <= 0 || mcol_ <= 0) - throw std::runtime_error("Matrix must have minimum size of 1"); - - // set default values to zero - for (size_t i = 0; i < size(); i++) - data_[i] = 0; - }; - - -// copy constructor -// -Matrix::Matrix(const Matrix& other) - : nrow_(other.nrow_), mcol_(other.mcol_), - data_(std::make_unique(other.size())) { - - // Matrix values have already been validated - for (size_t i = 0; i < size(); i++) - data_[i] = other.data_[i]; -} - - -Matrix::Matrix(Matrix&& other) - : nrow_(other.nrow_), mcol_(other.mcol_), data_(std::move(other.data_)) {}; - - -double Matrix::operator()(const size_t& i, const size_t& j) const { - return data_[mat_idx_to_array_(i, j)]; -} - -double& Matrix::operator()(const size_t& i, const size_t& j) { - return data_[mat_idx_to_array_(i, j)]; -} - -std::array Matrix::dims() const { - return {nrow_, mcol_}; -} - - -size_t Matrix::mat_idx_to_array_(const size_t& i, const size_t& j) const { - if (i >= nrow_ || j >= mcol_) - throw std::runtime_error("Indices must be postive integers or zero."); - - return i*mcol_ + j; -} - - -size_t Matrix::size() const { return nrow_ * mcol_; }; diff --git a/src/bcfio.cpp b/src/bcfio.cpp new file mode 100644 index 0000000..93294ce --- /dev/null +++ b/src/bcfio.cpp @@ -0,0 +1,181 @@ +// +// By: Robert Vogel +// Affiliation: Palmer Lab at UCSD +// Date: 2025-01-09 +// +// Input argument +// filename: vcf with haplotpye +// +// +// + +#include + + +/////////////////////////////////////////////////////////////////// +// BcfHeader +/////////////////////////////////////////////////////////////////// +// +int bcfio::BcfHeader::decode_hts_idinfo_(const char *name, + const int bcf_dt_type, + bcfio::BcfHdrAttr *ptr) const { + + // BCF_DT_ID is the C macro for the ID dictionary index defined + // by htslib see htslib/vcf.h line 86 + int idx = htslib::bcf_hdr_id2int(hdr_, BCF_DT_ID, name); + + if (idx < 0) + return idx; + + uint64_t val = hdr_->id[BCF_DT_ID][idx].val->info[bcf_dt_type]; + + ptr->number = val >> 12 & 0xfffff; + ptr->vl_type = val >> 8 & 0xf; + ptr->type = val >> 4 & 0xf; + ptr->coltype = val & 0xf; + + return 0; +} + +int bcfio::BcfHeader::get_format_attr(const char *id, BcfHdrAttr *ptr) const { + return decode_hts_idinfo_(id, BCF_HL_FMT, ptr); +} + +int bcfio::BcfHeader::get_info_attr(const char *id, BcfHdrAttr *ptr) const { + return decode_hts_idinfo_(id, BCF_HL_INFO, ptr); +} + +int bcfio::BcfHeader::get_filter_attr(const char *id, BcfHdrAttr *ptr) const { + return decode_hts_idinfo_(id, BCF_HL_FLT, ptr); +} + +int32_t bcfio::BcfHeader::k_fmt(const char *id) const { + if (!id) + return -1; + + BcfHdrAttr fmt {}; + + int32_t status { 0 }; + + if ((status = get_format_attr(id, &fmt)) < 0) + return status; + + return static_cast(fmt.number); +} + +int bcfio::BcfHeader::subset_samples(const char *filename) { + return htslib::bcf_hdr_set_samples(hdr_, filename, 1); +} + +const std::unique_ptr bcfio::BcfHeader::sample_names() const { + + std::unique_ptr samp_names = + std::make_unique(n_samples()); + + for (size_t i = 0; i < n_samples(); i++) + samp_names[i] = std::string(*(hdr_->samples + i)); + + return samp_names; +} + +/////////////////////////////////////////////////////////////////// +// BcfFloatRecord +/////////////////////////////////////////////////////////////////// + +bcfio::BcfFloatRecord::~BcfFloatRecord() { + if (rec_) htslib::bcf_destroy(rec_); + if (dst_) free(dst_); + rec_ = nullptr; + dst_ = nullptr; +} + +std::optional bcfio::BcfFloatRecord::get(const size_t row_idx, + const size_t col_idx) const { + size_t idx = row_idx * col_num_ + col_idx; + if (idx >= size()) return std::nullopt; + + return *(dst_ + idx); +} + +int bcfio::BcfFloatRecord::load_data(bcfio::BcfHeader *hdr, const char *id) { + int status { 0 }; + col_num_ = row_num_ = 0; + + status = htslib::bcf_get_format_values(hdr->hts_hdr(), + rec_, + id, + (void**)(&dst_), + &ndst_, + BCF_HT_REAL); + + if (status < 0) + return status; + + int32_t k { 0 }; + if ((k = hdr->k_fmt(id)) < 0) { + col_num_ = row_num_ = 0; + return k; + } + + col_num_ = static_cast(k); + row_num_ = hdr->n_samples(); + + return 0; +} + + +/////////////////////////////////////////////////////////////////// +// BcfRead +/////////////////////////////////////////////////////////////////// +/// + +bcfio::ReadBcf::ReadBcf() + : fname_(""), fid_(nullptr), hdr_() {}; + +bcfio::ReadBcf::ReadBcf(const char *filename, htslib::htsFile* fid) + : fname_(filename), + fid_(fid), + hdr_(fid_) {}; + +bcfio::ReadBcf::~ReadBcf() { + if (fid_) htslib::hts_close(fid_); +} + +bool bcfio::ReadBcf::isopen() const { + return fid_ != nullptr; +} + +// TODO: subset samples by those in sample_fname file +int bcfio::ReadBcf::set_samples(const char *sample_fname) { + + // Subset samples with those found in the file sample_fname + if (!sample_fname || *sample_fname == '\0') { + fprintf(stdout, "No file with sample names detected, retreiving" + " records for all samples.\n"); + return -1; + } + + return hdr_.subset_samples(sample_fname); +}; + +// title: load next record +int bcfio::ReadBcf::next_record(bcfio::BcfFloatRecord *ptr, const char *id) { + int status = htslib::bcf_read(fid_, hdr_.hts_hdr(), ptr->cur_rec()); + if (status != 0) + return status; + + // Unpacking options defined in htslib/vcf.h line 419 + if (htslib::bcf_unpack(ptr->cur_rec(), BCF_UN_ALL) < 0) + return -1; + + return ptr->load_data(&hdr_, id); +} + +bcfio::ReadBcf bcfio::open(const char* filename, const char* mode) { + htslib::htsFile* fid = htslib::hts_open(filename, mode); + if (!fid) + return bcfio::ReadBcf(); + + return bcfio::ReadBcf(filename, fid); +} + diff --git a/src/grm.cpp b/src/grm.cpp new file mode 100644 index 0000000..f0b83b3 --- /dev/null +++ b/src/grm.cpp @@ -0,0 +1,531 @@ +// Palmer Lab at UCSD +// +// +// ACKNOWLEDGMENT +// +// Code design and original version completed by Robert Vogel, +// reviewed by Claude Opus 4.6, the AI assistant from Anthropic +// with minor recommendations incorporated. +// +// + +#include "grm.h" + + +//////////////////////////////////////////////////////////////////// +// COORDINATES CLASS +//////////////////////////////////////////////////////////////////// + +grm::Coordinates::Coordinates(Coordinates&& other) + : contig(std::move(other.contig)), + len(other.len), + pos(std::move(other.pos)) { + + other.len = 0; + other.contig = ""; +} + +grm::Coordinates& grm::Coordinates::operator=(Coordinates&& other) { + if (this == &other) + return *this; + + len = other.len; + contig = std::move(other.contig); + pos = std::move(other.pos); + + other.len = 0; + other.contig = ""; + + return *this; +} + +// remember that Coordinates* should be uninstantiated +grm::STATUS grm::write(io::FileIO* fio, const Coordinates* coords) { + + if (!fio || !fio->fid || !coords) + return grm::ERROR_NULLPTR_ARG; + + size_t nwritten = 0; + + // write contig name to file + uint64_t nchar = coords->contig.size(); + nwritten = fwrite(&nchar, sizeof(nchar), 1, fio->fid); + if (nwritten != 1) + return grm::ERROR_ON_WRITE; + + nwritten = fwrite(coords->contig.c_str(), + sizeof(char), + nchar, + fio->fid); + if (static_cast(nwritten) != nchar) + return grm::ERROR_ON_WRITE; + + // write positions + uint64_t npos = coords->len; + nwritten = fwrite(&npos, sizeof(npos), 1, fio->fid); + if (nwritten != 1) + return grm::ERROR_ON_WRITE; + + nwritten = fwrite(coords->pos.get(), + sizeof(uint64_t), + npos, + fio->fid); + if (static_cast(nwritten) != npos) + return grm::ERROR_ON_WRITE; + + return grm::SUCCESS; +} + + +grm::STATUS grm::read(io::FileIO* fio, Coordinates* coords) { + + if (!fio || !fio->fid || !coords) + return grm::ERROR_NULLPTR_ARG; + + // I create a temporary Coordinates class, because I don't want + // the input coords instance to partial update upon an error + grm::Coordinates tmpc {}; + + // will store the number of bytes read at each step + size_t nread = 0; + + // read contig name + uint64_t size_contig_name = 0; + nread = fread(&size_contig_name, sizeof(size_contig_name), 1, fio->fid); + if (nread != 1) + return grm::ERROR_ON_READ; + + std::unique_ptr buffer = std::make_unique(size_contig_name + 1); + std::memset(buffer.get(), '\0', size_contig_name + 1); + + nread = fread(buffer.get(), sizeof(char), size_contig_name, fio->fid); + if (static_cast(nread) != size_contig_name) + return grm::ERROR_ON_READ; + + tmpc.contig = std::string(buffer.get(), size_contig_name); + + // read in positions + uint64_t npos = 0; + nread = fread(&npos, sizeof(npos), 1, fio->fid); + if (nread != 1) + return grm::ERROR_ON_READ; + + tmpc.len = npos; + + tmpc.pos = std::make_unique(npos); + nread = fread(tmpc.pos.get(), sizeof(uint64_t), npos, fio->fid); + if (static_cast(nread) != npos) + return grm::ERROR_ON_READ; + + *coords = std::move(tmpc); + + return grm::SUCCESS; +} + +//////////////////////////////////////////////////////////////////// +// SAMPLES CLASS +//////////////////////////////////////////////////////////////////// + +grm::Samples::Samples(grm::Samples&& other) + : len(other.len), names(std::move(other.names)) { + other.len = 0; +} + + +grm::Samples& grm::Samples::operator=(grm::Samples&& other) { + if (this == &other) + return *this; + + len = other.len; + names = std::move(other.names); + + other.len = 0; + + return *this; +} + + +grm::STATUS grm::write(io::FileIO* fio, const grm::Samples* samples) { + + if (!fio || !fio->fid || !samples) + return grm::ERROR_NULLPTR_ARG; + + size_t nwritten = 0; + uint64_t nsamps = samples->len; + nwritten = fwrite(&nsamps, sizeof(nsamps), 1, fio->fid); + if (nwritten != 1) + return grm::ERROR_ON_WRITE; + + // When it comes time to read the data, I need to make a character + // buffer to temporarily place the read string. To make this + // buffer, I need to know the length of string with the greatest + // number of characters. Here I find that number and store in + // the binary file. + + uint64_t nchar_max = 0; + uint64_t tmp = 0; + for (uint64_t n = 0; n < nsamps; n++) + if ((tmp = samples->names[n].size()) > nchar_max) nchar_max = tmp; + + if (nchar_max == 0) + return grm::ERROR_INVALID_ARG; + + nwritten = fwrite(&nchar_max, sizeof(nchar_max), 1, fio->fid); + if (nwritten != 1) + return grm::ERROR_ON_WRITE; + + // Write each string to file; + uint64_t nchar = 0; + for (uint64_t n = 0; n < nsamps; n++) { + nchar = samples->names[n].size(); + + nwritten = fwrite(&nchar, sizeof(nchar), 1, fio->fid); + if (nwritten != 1) + return grm::ERROR_ON_WRITE; + + nwritten = fwrite(samples->names[n].c_str(), + sizeof(char), + nchar, + fio->fid); + if (nwritten != nchar) + return grm::ERROR_ON_WRITE; + } + + return grm::SUCCESS; +} + +grm::STATUS grm::read(io::FileIO* fio, grm::Samples* samples) { + + if (!fio || !fio->fid || !samples) + return grm::ERROR_NULLPTR_ARG; + + // I create a temporary Sample class, because I don't want the + // input samples instance to partial update upon an error + grm::Samples tmp_samps {}; + + size_t nread; + uint64_t n_samples = 0; + + nread = fread(&n_samples, sizeof(n_samples), 1, fio->fid); + if (nread != 1) + return grm::ERROR_ON_READ; + + tmp_samps.len = n_samples; + tmp_samps.names = std::make_unique(n_samples); + + // Get the number of characters of the longest string + size_t nchar_max = 0; + nread = fread(&nchar_max, sizeof(size_t), 1, fio->fid); + if (nread != 1) + return grm::ERROR_ON_READ; + + + std::unique_ptr buffer = std::make_unique(nchar_max+1); + std::memset(buffer.get(), '\0', nchar_max + 1); + + // read sample names + uint64_t nchar = 0; + for (uint64_t n = 0; n < n_samples; n++) { + nread = fread(&nchar, sizeof(nchar), 1, fio->fid); + if (nread != 1) + return grm::ERROR_ON_READ; + + nread = fread(buffer.get(), sizeof(char), nchar, fio->fid); + if (static_cast(nread) != nchar) + return grm::ERROR_ON_READ; + + tmp_samps.names[n] = std::string(buffer.get(), nchar); + + std::memset(buffer.get(), '\0', nchar); + nchar = 0; + } + + *samples = std::move(tmp_samps); + + return grm::SUCCESS; +} + +//////////////////////////////////////////////////////////////////// +// HDR CLASS +//////////////////////////////////////////////////////////////////// + +grm::Hdr::Hdr() + : prog_version(constants::PROG_VERSION), + file_version(grm::FILE_VERSION), + grm_type(UNSPECIFIED), + coords(std::make_unique()), + samples(std::make_unique()) {}; + + +grm::Hdr::Hdr(Hdr&& other) + : prog_version(other.prog_version), + file_version(other.file_version), + grm_type(other.grm_type), + coords(nullptr), samples(nullptr) { + + other.prog_version = constants::PROG_VERSION; + other.file_version = grm::FILE_VERSION; + other.grm_type = grm::UNSPECIFIED; + + coords = std::move(other.coords); + samples = std::move(other.samples); +} + +grm::Hdr& grm::Hdr::operator=(Hdr&& other) { + if (this == &other) + return *this; + + prog_version = other.prog_version; + other.prog_version = constants::PROG_VERSION; + + file_version = other.file_version; + other.file_version = FILE_VERSION; + + grm_type = other.grm_type; + other.grm_type = grm::UNSPECIFIED; + + coords = std::move(other.coords); + samples = std::move(other.samples); + return *this; +} + + +grm::STATUS grm::write(io::FileIO* fio, const Hdr* hdr) { + + if (!fio || !fio->fid || !hdr) + return grm::ERROR_NULLPTR_ARG; + + size_t nwritten = 0; + uint32_t tmp_version = hdr->prog_version.pack(); + nwritten = fwrite(&tmp_version, sizeof(tmp_version), 1, fio->fid); + if (nwritten != 1) + return grm::ERROR_ON_WRITE; + + tmp_version = hdr->file_version.pack(); + nwritten = fwrite(&tmp_version, sizeof(tmp_version), 1, fio->fid); + if (nwritten != 1) + return grm::ERROR_ON_WRITE; + + nwritten = fwrite(&hdr->grm_type, sizeof(GrmType), 1, fio->fid); + if (nwritten != 1) + return grm::ERROR_ON_WRITE; + + grm::STATUS status = grm::FAILED; + if ((status = write(fio, hdr->coords.get())) != grm::SUCCESS) + return status; + + if ((status = write(fio, hdr->samples.get())) != grm::SUCCESS) + return status; + + return grm::SUCCESS; +} + +grm::STATUS grm::read(io::FileIO* fio, Hdr* hdr) { + + if (!fio || !fio->fid || !hdr) + return grm::ERROR_NULLPTR_ARG; + + Hdr tmp_hdr {}; + + size_t nread = 0; + uint32_t version = 0; + nread = fread(&version, sizeof(version), 1, fio->fid); + if (nread != 1) + return grm::ERROR_ON_READ; + tmp_hdr.prog_version = utils::Version::unpack(version); + + nread = fread(&version, sizeof(version), 1, fio->fid); + if (nread != 1) + return grm::ERROR_ON_READ; + tmp_hdr.file_version = utils::Version::unpack(version); + + nread = fread(&tmp_hdr.grm_type, sizeof(grm::GrmType), 1, fio->fid); + if (nread != 1) + return grm::ERROR_ON_READ; + + grm::STATUS status = grm::FAILED; + status = read(fio, tmp_hdr.coords.get()); + if (status != grm::SUCCESS) + return status; + + status = read(fio, tmp_hdr.samples.get()); + if (status != grm::SUCCESS) + return status; + + *hdr = std::move(tmp_hdr); + return grm::SUCCESS; +} + + +//////////////////////////////////////////////////////////////////// +// GRM CLASS +//////////////////////////////////////////////////////////////////// +// +// Recall that the GRM is a symmetric matrix, therefore we only need +// to store the upper triagonal and diagonal element values. +// Consequently, the size of the array storing the data is n*(n +1)/2. +// + +grm::Grm::Grm(): n_samples(0), data(nullptr) {}; + +grm::Grm::Grm(uint64_t n_samps) + : n_samples(n_samps), + data(size() != 0 ? std::make_unique(size()) : nullptr) { + + if (data) + std::memset(data.get(), 0, size() * sizeof(float)); +} + +grm::Grm::Grm(grm::Grm&& other) + : n_samples(other.n_samples), data(std::move(other.data)) { + other.n_samples = 0; +}; + +grm::Grm& grm::Grm::operator=(grm::Grm&& other) { + if (this == &other) + return *this; + + n_samples = other.n_samples; + other.n_samples = 0; + + data = std::move(other.data); + + return *this; +} + +// The number of upper diagonal + diagonal elements of the GRM +uint64_t grm::Grm::size() const { + return n_samples * (n_samples + 1) / 2; +}; + + +grm::STATUS grm::Grm::midx_to_arr(const uint64_t i, const uint64_t j, + uint64_t* idx) const { + + if (i >= n_samples || j >= n_samples) + return grm::ERROR_IDX_ARR_BOUNDS; + + // remember that by symmetry, the matrix is equal to its transpose + if (i <= j) + *idx = MATRIX_IDX_TO_ARRAY(i, j, n_samples); + else + *idx = MATRIX_IDX_TO_ARRAY(j, i, n_samples); + + return grm::SUCCESS; +} + + +float grm::Grm::operator()(const uint64_t i, const uint64_t j) const { + if (i > j) + return data[MATRIX_IDX_TO_ARRAY(j, i, n_samples)]; + return data[MATRIX_IDX_TO_ARRAY(i, j, n_samples)]; +} + + +float& grm::Grm::operator()(const uint64_t i, const uint64_t j) { + if (i > j) + return data[MATRIX_IDX_TO_ARRAY(j, i, n_samples)]; + return data[MATRIX_IDX_TO_ARRAY(i, j, n_samples)]; +} + + +grm::STATUS grm::Grm::get(const uint64_t i, const uint64_t j, float *val) const { + uint64_t idx = 0; + grm::STATUS status = grm::FAILED; + if ((status = midx_to_arr(i, j, &idx)) != grm::SUCCESS) + return status; + + *val = data[idx]; + + return status; +} + + +grm::STATUS grm::Grm::set(const uint64_t i, const uint64_t j, const float val) { + uint64_t idx = 0; + grm::STATUS status = grm::FAILED; + if ((status = midx_to_arr(i, j, &idx)) != grm::SUCCESS) + return status; + + data[idx] = val; + + return status; +} + + +grm::STATUS grm::write(io::FileIO *fio, + const grm::Hdr* hdr, const grm::Grm* grmatrix) { + + if (!fio || !fio->fid || !hdr || !grmatrix) + return grm::ERROR_NULLPTR_ARG; + + size_t nwritten = 0; + + // write file type specification + nwritten = fwrite(&FILE_TYPE_SPEC, + sizeof(FILE_TYPE_SPEC), 1, fio->fid); + if (nwritten != 1) + return grm::ERROR_ON_WRITE; + + grm::STATUS status = grm::FAILED; + + // srite meta data stored in header; + if ((status = grm::write(fio, hdr)) != grm::SUCCESS) + return status; + + uint64_t ndata = static_cast(grmatrix->size()); + nwritten = fwrite(&ndata, sizeof(ndata), 1, fio->fid); + if (nwritten != 1) + return grm::ERROR_ON_WRITE; + + nwritten = fwrite(grmatrix->data.get(), + sizeof(grmatrix->data[0]), + ndata, + fio->fid); + if (static_cast(nwritten) != ndata) + return grm::ERROR_ON_WRITE; + + return grm::SUCCESS; +} + + +grm::STATUS grm::read(io::FileIO* fio, + grm::Hdr* hdr, grm::Grm* grmatrix) { + + if (!fio || !fio->fid || !hdr || !grmatrix) + return grm::ERROR_NULLPTR_ARG; + + size_t nread = 0; + uint32_t ftype = 0; + + nread = fread(&ftype, sizeof(ftype), 1, fio->fid); + if (nread != 1) + return grm::ERROR_ON_READ; + + if (ftype != FILE_TYPE_SPEC) + return grm::ERROR_NOT_A_GRM_FILE; + + grm::STATUS status = grm::read(fio, hdr); + if (status != grm::SUCCESS) + return status; + + uint64_t n_samples = hdr->samples->len; + grm::Grm tmp_grm { n_samples }; + uint64_t ndata = 0; + nread = fread(&ndata, sizeof(ndata), 1, fio->fid); + if (nread != 1) + return grm::ERROR_ON_READ; + + if (ndata != n_samples*(n_samples + 1) / 2) + return grm::ERROR_ON_READ; + + nread = fread(tmp_grm.data.get(), + sizeof(tmp_grm.data[0]), ndata, fio->fid); + + if (nread != ndata) + return grm::ERROR_ON_READ; + + *grmatrix = std::move(tmp_grm); + + return grm::SUCCESS; +} diff --git a/src/logger.cpp b/src/logger.cpp new file mode 100644 index 0000000..ed0ffdc --- /dev/null +++ b/src/logger.cpp @@ -0,0 +1,89 @@ + +#include + + +Logger::Logger(): + t_(time(nullptr)), + time_point_(localtime(&t_)) { empty_bufs_(); }; + + +void Logger::empty_bufs_() { + std::memset(time_buf_, '\0', time_buf_len_); + std::memset(str_buf_, '\0', str_buf_len_); +} + + +int Logger::load_time_buf_() { + // get time and format time string + t_ = time(nullptr); + time_point_ = localtime(&t_); + + // strftime returns the number of characters written to buffer, + // a 0 returned indicates an error has occured. + return strftime(time_buf_, time_buf_len_,"%FT%H:%M:%S", time_point_); +} + + +void Logger::vprintf_(FILE *stream, const char *log_type, const char *format, + va_list arg_ptr) { + + if ((status_ = load_time_buf_()) <= 0) + strncpy(str_buf_, + "logger time buffer failure, please notify maintainer.", + max_str_); + else if ((status_ = vsnprintf(str_buf_, max_str_, format, arg_ptr)) <= 0) + strncpy(str_buf_, + "logger msg failure, please notify maintainer.", + max_str_); + else if (status_ >= max_str_) { + status_ = -status_; + strncpy(str_buf_, + "logger msg too long, please shorten msg.", + max_str_); + } + + if (status_ <= 0) + fprintf(stderr, "%s\t%s\t%s\n", time_buf_, err_str_, str_buf_); + else + fprintf(stream, "%s\t%s\t%s\n", time_buf_, log_type, str_buf_); + + empty_bufs_(); +} + + +int Logger::info(const char *format, ...) { + va_list arg_ptr; + va_start(arg_ptr, format); + vprintf_(stdout, info_str_, format, arg_ptr); + va_end(arg_ptr); + return status_; +} + +int Logger::warn(const char *format, ...) { + va_list arg_ptr; + va_start(arg_ptr, format); + vprintf_(stdout, warn_str_, format, arg_ptr); + va_end(arg_ptr); + return status_; +} + +int Logger::error(const char *format, ...) { + va_list arg_ptr; + va_start(arg_ptr, format); + vprintf_(stderr, err_str_, format, arg_ptr); + va_end(arg_ptr); + return status_; +} + +// int Logger::info(const char *msg) { +// return info("%s", msg); +// } +// +// int Logger::warn(const char *msg) { +// return warn("%s", msg); +// } +// +// int Logger::error(const char *msg) { +// return error("%s", msg); +// } + diff --git a/src/main.cpp b/src/main.cpp index 3b81cb4..d440da4 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -1,171 +1,218 @@ // Compute the genomic relationship matrix using haplotypes // -// By: Robert Vogel -// Affiliation: Palmer Lab at UCSD -// Date: 2025-01-09 +// Palmer Lab at UCSD // -// Input argument -// filename: vcf with haplotpye +// This program performs a single-pass computation of the genetic relationship +// matrix (GRM, GR matrix). GR matrices may be constructed using alt allele +// counts, expected alt allele counts, expected haplotype counts, or a +// combination of both expected alt allele and haplotype counts. // -// This program performs a single-pass computation of the -// haplotype based genomic relationship matrix. The approach -// is well defined for the covariance, however under my definition -// of the haplotype based covariance I had to derive the recursion -// relations myself. -// -// -// -// Acknowledgment -// -// Code design and original version completed by Robert Vogel, -// reviewed by Claude Sonnet, the AI assistant from Anthropic -// (Jan 2025), with minor recommendations incorporated. -// -#include -#include -#include "HaplotypeVcfParser.h" +#include +#include +#include +#include +#include +#include -size_t MARKER_PRINT_INTERVAL { 1000 }; -char HELP_LONG_FLAG[] { "--help" }; -char HELP_SHORT_FLAG[] { "-h" }; +#define FAILED_CALC -1 +#define SUCCESS_CALC 0 -int main(int argc, char* argv[]) -{ - if (argc != 2 && argc != 3) - throw std::runtime_error("Must specify vcf"); - - if (argc == 2 - && (strcmp(argv[1], HELP_SHORT_FLAG) == 0 - || strcmp(argv[1], HELP_LONG_FLAG) == 0)) { - printf("hgrm - Compute GRM from expected haplotype counts.\n" - "Usage\n" - "\n" - " hgrm []\n" - "\n" - "Options\n" - " output_matrix_filename Filename to print covariance matrix\n" - "\n" - "Description\n" - " A program to compute a genetic relationship matrix from a vcf\n" - " with expected haplotype counts record per sample per locus.\n"); +const size_t STR_BUF_LEN { 500 }; +char STR_BUF[STR_BUF_LEN]; +int main(int argc, char* argv[]) +{ + // if (argc != 2 && argc != 4) { + // fprintf(stderr, "Incorrect input, see --help for correct usage.\n"); + // exit(EXIT_FAILURE); + // } + + argparse::ArgParser parser { + "grm: Genetic Relationship Matrix", + "This program provides tools for computing the genetic relationship" + " matrix (GRM) and the leave-one-chromosome-out (LOCO) matrices for" + " linear mixed effect based association studies. The GRM may be" + " computed using called genotypes, expected alternative allele counts," + " expected haplotype counts, or both expected alternative allele" + " and haplotype counts. By default the expected alternative allele" + " counts are used." + }; + + argparse::CmdDef *contig_cmd = parser.add_cmd("contig"); + + contig_cmd->add_arg("-o", + argparse::ArgType::STRING, + "the path and filename that the resulting haplotype genetic" + " relationship matrix is printed."); + + contig_cmd->add_arg("--sample_names", + argparse::ArgType::STRING, + "The path and name of the file containing sample names to be" + " included in computing the relationship matrix. The file must" + " include a single sample filename, and if necessary file system" + " path, per line."); + + contig_cmd->add_arg("--gt", + argparse::ArgType::BOOLEAN, + "Use sample genotypes to compute the relationship matrix"); + + contig_cmd->add_arg("--ehc", + argparse::ArgType::BOOLEAN, + "Use sample expected haplotype count to compute the the genetic" + " relationship matrix."); + + contig_cmd->add_arg("-b", + argparse::ArgType::BOOLEAN, + "Use both the expected alternative allele and haplotype counts to" + " compute the genetic relationship matrix"); + + contig_cmd->add_arg("bcf", + argparse::ArgType::STRING, + "The path and filename of the genetic data to compute the GRM. The" + " data may be in any of the htslib supported formats, i.e. vcf," + " vcf.gz, or bcf."); + + + argparse::CmdDef *loco_cmd = parser.add_cmd("loco"); + loco_cmd->add_arg("filename", + argparse::ArgType::STRING, + "Name, and path, of file that stores the name and paths of matrix" + " files used to compute leave-one-chromosome-out (LOCO) relationship" + " matrix."); + + + + Logger log {}; + int status = FAILED_CALC; + argparse::ArgStatus arg_status = parser.parse_args(argc, argv); + + // PARSE ARGUMENTS + if (arg_status == argparse::ArgStatus::HELP) return 0; - } - - char* filename_input { argv[1] }; - char* filename_output { nullptr }; - - if(argc == 3) - filename_output = argv[2]; - - - const std::chrono::time_point timer - { std::chrono::steady_clock::now() }; - - - fprintf(stdout, "Allocating memory\n"); - - // open VCF file and parse meta data and header - HaplotypeVcfParser vcf_data { filename_input, 100000 }; - - - // instantiate matrices to hold calculations - Matrix covariance { vcf_data.n_samples(), vcf_data.n_samples() }; - - // instantiate record object - HaplotypeDataRecord record { vcf_data.n_samples(), vcf_data.k_founders() }; - - // analyze each line, i.e. position, in the VCF - size_t m_markers { 1 }; - double sum { 0 }; - const double* rowi { nullptr }; - const double* rowj { nullptr }; - double* rowi_cov { nullptr }; - const size_t k_founders { vcf_data.k_founders() }; - const size_t n_samples { vcf_data.n_samples() }; - - std::chrono::steady_clock::duration delta_t - { std::chrono::steady_clock::now() - timer }; - - fprintf(stdout, "Computing matrix, elapsed time %lld second(s)\n", - std::chrono::duration_cast(delta_t).count()); - - while(vcf_data.load_record(record)) { - - // for each founder, compute first and second moments - for (size_t i = 0; i < n_samples; i++) { - - rowi = &record(i, 0); - rowi_cov = &covariance(i, 0); - - for (size_t j = i; j < n_samples; j++) { - - rowj = &record(j,0); - sum = 0; + if (arg_status != argparse::ArgStatus::SUCCESS) { + log.error("Error: couldn't parse command line args, exiting\n"); + exit(EXIT_FAILURE); + } - for (int k = 0; k < k_founders; k++) - sum += rowi[k] * rowj[k]; + // PARSE ARGS FOR RESPECTIVE SUBPROGRAMS AND RUN + // + // Compute the GRM for the specified contig + if (parser.is_sub_cmd("contig")) { - rowi_cov[j] += sum; - } + std::optional tmp_str {}; + if((tmp_str = parser.get("bcf")) == std::nullopt) { + log.error("Error retrieving vcf name"); + exit(EXIT_FAILURE); } + std::string bcf_fname { tmp_str.value() }; - if (m_markers % MARKER_PRINT_INTERVAL == 0) { - delta_t = std::chrono::steady_clock::now() - timer; - - fprintf(stdout, "Completed %zu marker loci, elapsed time %lld second(s)\n", - m_markers, - std::chrono::duration_cast(delta_t).count()); + if ((tmp_str = parser.get("o")) == std::nullopt) { + log.error("Error retrieving output name"); + exit(EXIT_FAILURE); } + std::string out_fname { tmp_str.value() }; - m_markers++; - - } + if (out_fname.size() == 0) + out_fname = bcf_fname + ".mat"; + if ((tmp_str = parser.get("sample_names")) == std::nullopt) { + log.error("Error retrieving sample_names file.\n"); + exit(EXIT_FAILURE); + } - FILE* fout = stdout; + std::string samp_fname { tmp_str.value() }; - if (argc == 3 && filename_output != nullptr) { + + std::optional tmp_bool {}; + if ((tmp_bool = parser.get("gt")) == std::nullopt) { + log.error("Error retrieving relationship matrix type.\n"); + exit(EXIT_FAILURE); + } + bool use_gt { tmp_bool.value() }; - if ((fout = fopen(filename_output, "w")) == nullptr) - throw std::runtime_error("Error in opening file for writing."); + if ((tmp_bool = parser.get("ehc")) == std::nullopt) { + log.error("Error retrieving relationship matrix type.\n"); + exit(EXIT_FAILURE); + } + bool use_ehc { tmp_bool.value() }; - delta_t = std::chrono::steady_clock::now() - timer; - fprintf(stdout, "Writing results to file %s, elapsed time %lld second(s)\n", - filename_output, - std::chrono::duration_cast(delta_t).count()); + if ((tmp_bool = parser.get("b")) == std::nullopt) { + log.error("Error retrieving relationship matrix type.\n"); + exit(EXIT_FAILURE); + } + bool use_both { tmp_bool.value() }; - } else if (argc == 3 && filename_output == nullptr) - throw std::runtime_error("Output filename is not specified"); + if ((use_gt && use_both) || (use_gt && use_ehc) || (use_both && use_ehc)) { + log.error("user must specify either use_gt, use_both, use_ehc," + " or omit both options to compute the haplotype based" + " relationship matrix."); + exit(EXIT_FAILURE); + } + + log.info("BCF/VCF file name: %s", bcf_fname.c_str()); + + + bcfio::ReadBcf bfid { bcf_fname.c_str() }; + + int bstatus = 0; + if (samp_fname.size() == 0) + log.info("Sample file: None, use all samples"); + else if ((bstatus = bfid.set_samples(samp_fname.c_str())) == 0) + log.info("Sample file: %s", samp_fname.c_str()); + else if (bstatus < 0) { + log.error("Subsetting by sample file, %s, resulted in error", + samp_fname.c_str()); + exit(EXIT_FAILURE); + } else if (bstatus > 0) { + log.error("One or more samples specified in sample file, %s," + " do not %s", + samp_fname.c_str(), + bcf_fname.c_str()); + exit(EXIT_FAILURE); + } - size_t i { 0 }; - size_t j { 0 }; - for (i = 0; i < n_samples; i++) { + log.info("Output matrix file: %s", out_fname.c_str()); + + grm::Grm grmatrix { bfid.n_samples() }; + + if (use_gt) { + log.info("Relationship matrix: genotype"); + status = compute_genotype_matrix(); + } else if (use_ehc) { + log.info("Relationship matrix: expected haplotype count"); + status = compute_ehc_matrix(&log, &bfid, &grmatrix); + } else if (use_both) { + log.info("Relationship matrix: expected alt allele and haplotype" + " counts"); + status = compute_eac_and_ehc_matrix(); + } else { + log.info("Relationship matrix: expected alternative allele counts"); + status = compute_eac_matrix(); + } - for (j = 0; j < n_samples-1; j++) { - if (j < i) - fprintf(fout, "%0.5f,", covariance(j,i)); - else - fprintf(fout, "%0.5f,", covariance(i,j)); + if (status == FAILED_CALC) + log.error("Computation failed"); - } + log.info("Writing to file"); - fprintf(fout,"%0.5f\n", covariance(i, j)); + grmatrix.write(out_fname); + } else if (parser.is_sub_cmd("loco")) { + // Compute the leave-one-chromosome-out matrix given a set of + // matricies. + // + printf("loco selected\n"); + } else if (parser.is_sub_cmd("assoc")) { + printf("association statistics selected\n"); } - fclose(fout); - - delta_t = std::chrono::steady_clock::now() - timer; - fprintf(stdout, "Done, elapsed time %lld second(s)\n", - std::chrono::duration_cast(delta_t).count()); - return 0; + return status; } diff --git a/src/use_both.cpp b/src/use_both.cpp new file mode 100644 index 0000000..807097e --- /dev/null +++ b/src/use_both.cpp @@ -0,0 +1,7 @@ +// Compute GRM with expected alt allele and haplotype counts +// +#include + +int compute_eac_and_ehc_matrix() { + return -1; +} diff --git a/src/use_eac.cpp b/src/use_eac.cpp new file mode 100644 index 0000000..984ebdf --- /dev/null +++ b/src/use_eac.cpp @@ -0,0 +1,6 @@ + +#include + +int compute_eac_matrix() { + return -1; +} diff --git a/src/use_ehc.cpp b/src/use_ehc.cpp new file mode 100644 index 0000000..79d0081 --- /dev/null +++ b/src/use_ehc.cpp @@ -0,0 +1,141 @@ + + +#include + +int compute_ehc_matrix(Logger *log, bcfio::ReadBcf *bfid, Grm *cov) { + + int output_status = 0; + + // instantiate matrices to hold calculations + const size_t n_samples { bfid->n_samples() }; + int32_t k { 0 }; + if ((k = bfid->k_fmt("HD")) < 0) { + log->error("%s\n", "Wrong format id tag"); + exit(EXIT_FAILURE); + } + const size_t k_haps { static_cast(k) }; + + size_t idx_row { 0 }; + // size_t idx_col { 0 }; + size_t idx_hap { 0 }; + size_t idx_rec { 0 }; + + bcfio::BcfFloatRecord rec {}; + + std::optional val { 0 }; + + while (bfid->next_record(&rec, "HD") == 0) { + + for (idx_row = 0; idx_row < n_samples; idx_row++) { + + for (idx_hap = 0; idx_hap < k_haps; idx_hap++) { + + if ((val = rec.get(idx_row, idx_hap)) == std::nullopt) { + printf("IDX: (%zu, %zu) = null\n", idx_row, idx_hap); + return -1; + } + + (*cov)(idx_row, idx_row) += val.value(); + + printf("%f\t ", val.value()); + } + + printf("\n"); + } + + log->info("Processed %s records", + std::to_string(++idx_rec).c_str()); + + + } + + return output_status; +} + +// size_t m_markers { 1 }; +// +// double sum { 0 }; +// const double* rowi { nullptr }; +// const double* rowj { nullptr }; +// double* rowi_cov { nullptr }; +// const size_t k_founders { vcf_data.k_founders() }; +// const size_t n_samples { vcf_data.n_samples() }; +// +// std::chrono::steady_clock::duration delta_t +// { std::chrono::steady_clock::now() - timer }; +// +// fprintf(stdout, "Computing matrix, elapsed time %lld second(s)\n", +// std::chrono::duration_cast(delta_t).count()); +// +// while(vcf_data.load_record(record)) { +// +// // for each founder, compute first and second moments +// for (size_t i = 0; i < n_samples; i++) { +// +// rowi = &record(i, 0); +// rowi_cov = &covariance(i, 0); +// +// for (size_t j = i; j < n_samples; j++) { +// +// rowj = &record(j,0); +// sum = 0; +// +// for (int k = 0; k < k_founders; k++) +// sum += rowi[k] * rowj[k]; +// +// rowi_cov[j] += sum; +// } +// } +// +// if (m_markers % MARKER_PRINT_INTERVAL == 0) { +// delta_t = std::chrono::steady_clock::now() - timer; +// +// fprintf(stdout, "Completed %zu marker loci, elapsed time %lld second(s)\n", +// m_markers, +// std::chrono::duration_cast(delta_t).count()); +// } +// +// m_markers++; +// +// } +// +// +// FILE* fout = stdout; +// +// if (argc == 3 && filename_output != nullptr) { +// +// if ((fout = fopen(filename_output, "w")) == nullptr) +// throw std::runtime_error("Error in opening file for writing."); +// +// delta_t = std::chrono::steady_clock::now() - timer; +// fprintf(stdout, "Writing results to file %s, elapsed time %lld second(s)\n", +// filename_output, +// std::chrono::duration_cast(delta_t).count()); +// +// } else if (argc == 3 && filename_output == nullptr) +// throw std::runtime_error("Output filename is not specified"); +// +// +// size_t i { 0 }; +// size_t j { 0 }; +// for (i = 0; i < n_samples; i++) { +// +// for (j = 0; j < n_samples-1; j++) { +// if (j < i) +// fprintf(fout, "%0.5f,", covariance(j,i)); +// else +// fprintf(fout, "%0.5f,", covariance(i,j)); +// +// } +// +// fprintf(fout,"%0.5f\n", covariance(i, j)); +// } +// +// fclose(fout); +// +// +// delta_t = std::chrono::steady_clock::now() - timer; +// +// fprintf(stdout, "Done, elapsed time %lld second(s)\n", +// std::chrono::duration_cast(delta_t).count()); +// diff --git a/src/use_gt.cpp b/src/use_gt.cpp new file mode 100644 index 0000000..5abc58f --- /dev/null +++ b/src/use_gt.cpp @@ -0,0 +1,7 @@ +// Compute GRM with called genotypes + +#include + +int compute_genotype_matrix() { + return -1; +} diff --git a/src/utils.cpp b/src/utils.cpp deleted file mode 100644 index 46aebc2..0000000 --- a/src/utils.cpp +++ /dev/null @@ -1,301 +0,0 @@ -#include "../include/utils.h" - -const static size_t DEFAULT_BUFFER_SIZE { 1000 }; - - - - -CharBuffer::CharBuffer() - : buffer_size_(0), - buffer_idx_(0), - buffer_(nullptr) {} - - -CharBuffer::CharBuffer(size_t buffer_size) - : buffer_size_(buffer_size), - buffer_idx_(0), - buffer_(buffer_size > 0 ? std::make_unique(buffer_size+1) : nullptr) { - - if (buffer_ == nullptr) - throw std::runtime_error("CharBuffer needs to have length > 0"); - - buffer_[buffer_idx_] = '\0'; - buffer_[buffer_size_] = '\0'; -} - - -const char& CharBuffer::operator()(size_t idx) const { - if (idx >= size() || idx >= buffer_idx_) - throw std::out_of_range("index too large for buffer."); - - if (!buffer_) - throw std::runtime_error("error"); - return buffer_[idx]; -} - - -const size_t& CharBuffer::buffer_size() const { return buffer_size_; } -const size_t& CharBuffer::size() const { return buffer_idx_; } - - -void CharBuffer::append(char s) { - if (buffer_ == nullptr) - throw std::runtime_error("No data stored in buffer."); - - if (buffer_idx_ >= buffer_size_) - throw std::out_of_range("CharBuffer full"); - - buffer_[buffer_idx_++] = s; - buffer_[buffer_idx_] = '\0'; -} - - -void CharBuffer::reset(size_t buffer_size) { - buffer_size_ = buffer_size; - buffer_ = std::make_unique(buffer_size_); - buffer_[0] = '\0'; - buffer_[buffer_size_] = '\0'; - buffer_idx_ = 0; -} - - -void CharBuffer::reset() { - if (buffer_ == nullptr) - throw std::runtime_error("No data stored in buffer."); - - buffer_idx_ = 0; - buffer_[buffer_idx_] = '\0'; -} - - -const char* CharBuffer::data() const { - if (buffer_ == nullptr) - throw std::runtime_error("No data stored in buffer."); - - return buffer_.get(); -} - - - -StringRecord::StringRecord(const char delim) - : str_(nullptr), - delim_(delim), - buf_(DEFAULT_BUFFER_SIZE) { - - if (std::isspace(delim)) - is_delim_ = [](char c){ return std::isspace(c); }; - else - is_delim_ = [this](char c){ return char_is_delim_(c); }; -}; - -StringRecord::StringRecord(const char delim, const size_t buf_size) - : str_(nullptr), - delim_(delim), - buf_(buf_size) { - - if (std::isspace(delim)) - is_delim_ = [](char c){ return std::isspace(c); }; - else - is_delim_ = [this](char c){ return char_is_delim_(c); }; -}; - -StringRecord::StringRecord(const char delim, const char* str) - : str_(str), - delim_(delim), - buf_(DEFAULT_BUFFER_SIZE) { - - if (std::isspace(delim)) - is_delim_ = [](char c){ return std::isspace(c); }; - else - is_delim_ = [this](char c){ return char_is_delim_(c); }; - - for (;str_[size_] != '\0'; size_++) - ; -}; - - -bool StringRecord::char_is_delim_(char c) const { - if (c == delim_) - return true; - return false; -} - -// size including null character -void StringRecord::update_str(const char* str) { - reset(); - str_ = str; - - for (;str_[size_] != '\0'; size_++) - ; - - return ; -} - - -void StringRecord::reset() { - idx_ = 0; - size_ = 0; - str_ = nullptr; - buf_.reset(); -}; - - -size_t StringRecord::size() { - if (!str_) - return 0; - - return size_; -} - - -const char* StringRecord::data() const { - return buf_.data(); -} - - -bool StringRecord::next_field() { - if (!str_ || idx_ == size()) - return false; - - // eliminate preceeding white space - for(; idx_ < size(); idx_++) - if (!is_delim_(str_[idx_])) - break; - - buf_.reset(); - for (; idx_ < size(); idx_++) { - if (!is_delim_(str_[idx_])) - buf_.append(str_[idx_]); - - if (is_delim_(str_[idx_]) && !is_delim_(str_[idx_-1])) { - idx_++; - break; - } - } - - return true; -} - - -// bool StringRecord::is_delim_(char s) { -// return s == delim_; -// } - - -// buff_size is the count of characters to load -BufferedRead::BufferedRead(char* filename, size_t buff_size) - : filename_(filename), - buff_size_(buff_size), - fid_(std::fopen(filename, "r")), - buffer_(buff_size > 0 ? std::make_unique(buff_size+1) : nullptr) { - - if (!fid_) - throw std::runtime_error("File Access error"); - - if (std::feof(fid_)) - throw std::runtime_error("File is empty"); - - if (buffer_ == nullptr) - throw std::runtime_error("Buffer wasn't properly set"); - - buffer_pos_ = 0; - buffer_[0] = '\0'; -} - - -BufferedRead::~BufferedRead() { - if (fid_) - fclose(fid_); -} - -// Note: -// Return count == 0 indicates that we have reached the end of file. -// -size_t BufferedRead::get_line(CharBuffer& line_buffer) { - - line_buffer.reset(); - - size_t count { 0 }; - size_t buff_length { 1 }; - - if (buffer_[buffer_pos_] == '\0' - || buffer_pos_ == buff_size_) - buff_length = update_buffer_(); - - - if (buff_length == 0) - return count; - - while (buffer_[buffer_pos_] != '\n') { - - line_buffer.append(buffer_[buffer_pos_++]); - count++; - - if (buffer_[buffer_pos_] == '\0' || buffer_pos_ == buff_size_) - buff_length = update_buffer_(); - - if (buff_length == 0) - break; - } - - if (buffer_[buffer_pos_] == '\n') - buffer_pos_++; - - return count; -} - - -char BufferedRead::get_char() { - size_t buff_length { 1 }; - - - if (buffer_[buffer_pos_] == '\0') - buff_length = update_buffer_(); - - // when all entries of a file are read, detected by - // buff_length == 0, return the null character - if (buff_length == 0) - return '\0'; - - return buffer_[buffer_pos_++]; -} - - -size_t BufferedRead::update_buffer_() { - - if (std::ferror(fid_)) - throw std::runtime_error("File read error"); - - // note, if end of file, then we set n = 0; - size_t n { 0 }; - if (!std::feof(fid_)) - n = fread(buffer_.get(), sizeof(buffer_[0]), buff_size_, fid_); - - buffer_pos_ = 0; - buffer_[n] = '\0'; - - return n; -} - - -void BufferedRead::seek(size_t n) { - - int c = 0; - if ((c = std::fseek(fid_, n, SEEK_SET)) != 0) - throw std::runtime_error("Failed to relocate file stream to position."); - -} - - -size_t BufferedRead::tell() { - return std::ftell(fid_); -} - - -void BufferedRead::reset() { - buffer_pos_ = 0; - buffer_[buffer_pos_] = '\0'; - seek(0); -} - - diff --git a/tests/geno_test_data.bcf b/tests/geno_test_data.bcf new file mode 100644 index 0000000..6c77e17 Binary files /dev/null and b/tests/geno_test_data.bcf differ diff --git a/tests/test.vcf b/tests/geno_test_data.vcf similarity index 100% rename from tests/test.vcf rename to tests/geno_test_data.vcf diff --git a/tests/geno_test_data.vcf.gz b/tests/geno_test_data.vcf.gz new file mode 100644 index 0000000..43958ea Binary files /dev/null and b/tests/geno_test_data.vcf.gz differ diff --git a/tests/hd_01.csv b/tests/hd_01.csv new file mode 100644 index 0000000..8fd1fc2 --- /dev/null +++ b/tests/hd_01.csv @@ -0,0 +1,11 @@ +1.004,0,0.002,0.001,0,0.991,0.001,0.002 +0.001,0,0,0.989,0.005,0.005,1,0 +0.998,0,0.001,0,0,0,1,0 +0.84,0,0,0,0,0,1.159,0 +0.005,0.535,0.193,0.069,0.597,0.004,0.001,0.597 +0.001,0,0.001,1.002,0.02,0.013,0,0.962 +0.001,0.001,0,1.004,0.024,0.97,0,0.001 +0.001,0.955,0.907,0.125,0.003,0.004,0,0.005 +0.002,0.001,0.001,0.991,0.988,0.007,0.002,0.008 +0.249,0.001,0,1,0.748,0,0,0.001 +0.592,0,1,0,0,0,0,0.407 diff --git a/tests/hd_02.csv b/tests/hd_02.csv new file mode 100644 index 0000000..9e05481 --- /dev/null +++ b/tests/hd_02.csv @@ -0,0 +1,11 @@ +1.004,0,0.002,0.001,0,0.991,0.001,0.002 +0.001,0,0,0.989,0.005,0.005,1,0 +0.998,0,0.001,0,0,0,1.001,0 +0.84,0,0,0,0,0,1.159,0 +0.005,0.535,0.193,0.069,0.597,0.004,0.001,0.597 +0.001,0,0.001,1.002,0.02,0.013,0,0.962 +0.001,0.001,0,1.004,0.024,0.97,0,0.001 +0,0.955,0.907,0.125,0.003,0.004,0,0.005 +0.002,0.001,0.001,0.991,0.988,0.007,0.002,0.008 +0.249,0.001,0,1,0.748,0,0,0.001 +0.592,0,1,0,0,0,0,0.407 diff --git a/tests/hd_03.csv b/tests/hd_03.csv new file mode 100644 index 0000000..9e05481 --- /dev/null +++ b/tests/hd_03.csv @@ -0,0 +1,11 @@ +1.004,0,0.002,0.001,0,0.991,0.001,0.002 +0.001,0,0,0.989,0.005,0.005,1,0 +0.998,0,0.001,0,0,0,1.001,0 +0.84,0,0,0,0,0,1.159,0 +0.005,0.535,0.193,0.069,0.597,0.004,0.001,0.597 +0.001,0,0.001,1.002,0.02,0.013,0,0.962 +0.001,0.001,0,1.004,0.024,0.97,0,0.001 +0,0.955,0.907,0.125,0.003,0.004,0,0.005 +0.002,0.001,0.001,0.991,0.988,0.007,0.002,0.008 +0.249,0.001,0,1,0.748,0,0,0.001 +0.592,0,1,0,0,0,0,0.407 diff --git a/tests/hd_04.csv b/tests/hd_04.csv new file mode 100644 index 0000000..1305763 --- /dev/null +++ b/tests/hd_04.csv @@ -0,0 +1,11 @@ +1.004,0,0.002,0.001,0,0.991,0.001,0.002 +0.001,0,0,0.989,0.005,0.005,1,0 +0.998,0,0.001,0,0,0,1.001,0 +0.84,0,0,0,0,0,1.159,0 +.005,0.535,0.193,0.069,0.597,0.004,0.001,0.597 +0.001,0,0.001,1.002,0.02,0.013,0,0.962 +0.001,0.001,0,1.004,0.024,0.97,0,0.001 +0,0.955,0.907,0.125,0.003,0.004,0,0.005 +0.002,0.001,0.001,0.991,0.988,0.007,0.002,0.008 +0.249,0.001,0,1,0.748,0,0,0.001 +0.592,0,1,0,0,0,0,0.407 diff --git a/tests/hd_05.csv b/tests/hd_05.csv new file mode 100644 index 0000000..9e05481 --- /dev/null +++ b/tests/hd_05.csv @@ -0,0 +1,11 @@ +1.004,0,0.002,0.001,0,0.991,0.001,0.002 +0.001,0,0,0.989,0.005,0.005,1,0 +0.998,0,0.001,0,0,0,1.001,0 +0.84,0,0,0,0,0,1.159,0 +0.005,0.535,0.193,0.069,0.597,0.004,0.001,0.597 +0.001,0,0.001,1.002,0.02,0.013,0,0.962 +0.001,0.001,0,1.004,0.024,0.97,0,0.001 +0,0.955,0.907,0.125,0.003,0.004,0,0.005 +0.002,0.001,0.001,0.991,0.988,0.007,0.002,0.008 +0.249,0.001,0,1,0.748,0,0,0.001 +0.592,0,1,0,0,0,0,0.407 diff --git a/tests/hd_06.csv b/tests/hd_06.csv new file mode 100644 index 0000000..9e05481 --- /dev/null +++ b/tests/hd_06.csv @@ -0,0 +1,11 @@ +1.004,0,0.002,0.001,0,0.991,0.001,0.002 +0.001,0,0,0.989,0.005,0.005,1,0 +0.998,0,0.001,0,0,0,1.001,0 +0.84,0,0,0,0,0,1.159,0 +0.005,0.535,0.193,0.069,0.597,0.004,0.001,0.597 +0.001,0,0.001,1.002,0.02,0.013,0,0.962 +0.001,0.001,0,1.004,0.024,0.97,0,0.001 +0,0.955,0.907,0.125,0.003,0.004,0,0.005 +0.002,0.001,0.001,0.991,0.988,0.007,0.002,0.008 +0.249,0.001,0,1,0.748,0,0,0.001 +0.592,0,1,0,0,0,0,0.407 diff --git a/tests/hd_07.csv b/tests/hd_07.csv new file mode 100644 index 0000000..9e05481 --- /dev/null +++ b/tests/hd_07.csv @@ -0,0 +1,11 @@ +1.004,0,0.002,0.001,0,0.991,0.001,0.002 +0.001,0,0,0.989,0.005,0.005,1,0 +0.998,0,0.001,0,0,0,1.001,0 +0.84,0,0,0,0,0,1.159,0 +0.005,0.535,0.193,0.069,0.597,0.004,0.001,0.597 +0.001,0,0.001,1.002,0.02,0.013,0,0.962 +0.001,0.001,0,1.004,0.024,0.97,0,0.001 +0,0.955,0.907,0.125,0.003,0.004,0,0.005 +0.002,0.001,0.001,0.991,0.988,0.007,0.002,0.008 +0.249,0.001,0,1,0.748,0,0,0.001 +0.592,0,1,0,0,0,0,0.407 diff --git a/tests/hd_08.csv b/tests/hd_08.csv new file mode 100644 index 0000000..4bf0c43 --- /dev/null +++ b/tests/hd_08.csv @@ -0,0 +1,11 @@ +1.004,0,0.002,0.001,0,0.991,0.001,0.002 +0,0,0,0.989,0.005,0.005,1,0 +0.998,0,0.001,0,0,0,1.001,0 +0.84,0,0,0,0,0,1.159,0 +0.005,0.535,0.193,0.068,0.597,0.004,0.001,0.597 +0.001,0,0.001,1.002,0.02,0.013,0,0.962 +0.001,0.001,0,1.004,0.024,0.97,0,0.001 +0,0.955,0.908,0.125,0.003,0.004,0,0.005 +0.002,0.001,0.001,0.991,0.988,0.007,0.002,0.008 +0.249,0.001,0,1,0.748,0,0,0.001 +0.592,0,1,0,0,0,0,0.407 diff --git a/tests/main.cpp b/tests/main.cpp new file mode 100644 index 0000000..9ece626 --- /dev/null +++ b/tests/main.cpp @@ -0,0 +1,6 @@ +#include + +int main(int argc, char *argv[]) { + testing::InitGoogleTest(&argc, argv) ; + return RUN_ALL_TESTS(); +} diff --git a/tests/samples_test b/tests/samples_test new file mode 100644 index 0000000..fc07ba5 --- /dev/null +++ b/tests/samples_test @@ -0,0 +1,5 @@ +S03 +S05 +S06 +S08 +S11 diff --git a/tests/test_bcfio.cpp b/tests/test_bcfio.cpp new file mode 100644 index 0000000..1818ac8 --- /dev/null +++ b/tests/test_bcfio.cpp @@ -0,0 +1,499 @@ + +#include +#include +#include +#include +#include + +namespace htslib { +extern "C" { +#include +#include +} +} + +#include + + +char VCF_NAME[] { "build/geno_test_data.vcf" }; +char VCFGZ_NAME[] { "build/geno_test_data.vcf.gz" }; +char BCF_NAME[] { "build/geno_test_data.bcf" }; +uint8_t K_FOUNDERS = 8; +uint8_t N_SAMPS = 11; + + + +/////////////////////////////////////////////////////////////////////////// +// Test bcfio::BcfHeader +/////////////////////////////////////////////////////////////////////////// + +TEST(TestBcfHeader, ConstructorVcfHdr) { + htslib::htsFile* fid = htslib::hts_open(VCF_NAME, "r"); + bcfio::BcfHeader hdr { fid }; + + EXPECT_FALSE(hdr.isnull()); + + bcfio::BcfHdrAttr attr {}; + + int status = hdr.get_format_attr("HD", &attr); + EXPECT_EQ(status, 0); + EXPECT_EQ(attr.number, static_cast(K_FOUNDERS)); + + // This is a bit tricky. htslib/vcf.h encodes BCF_VL_FIXED + // and BCF_HT_REAL as integer macros. The values are bit packed + // in a uint64_t type. Consequently I need to cast the htslib + // macros to unsigned int + EXPECT_EQ(attr.vl_type, static_cast(BCF_VL_FIXED)); + EXPECT_EQ(attr.type, static_cast(BCF_HT_REAL)); + + if (fid) htslib::hts_close(fid); +} + +TEST(TestBcfHeader, ConstructorVcfGzHdr) { + htslib::htsFile* fid = htslib::hts_open(VCFGZ_NAME, "r"); + bcfio::BcfHeader hdr { fid }; + + EXPECT_FALSE(hdr.isnull()); + + bcfio::BcfHdrAttr attr {}; + + int status = hdr.get_format_attr("HD", &attr); + EXPECT_EQ(status, 0); + EXPECT_EQ(attr.number, static_cast(K_FOUNDERS)); + EXPECT_EQ(attr.vl_type, static_cast(BCF_VL_FIXED)); + EXPECT_EQ(attr.type, static_cast(BCF_HT_REAL)); + + if (fid) htslib::hts_close(fid); +} + +TEST(TestBcfHeader, ConstructorBcfHdr) { + htslib::htsFile *fid = htslib::hts_open(BCF_NAME, "r"); + bcfio::BcfHeader hdr { fid }; + + EXPECT_FALSE(hdr.isnull()); + + bcfio::BcfHdrAttr attr {}; + + int status = hdr.get_format_attr("HD", &attr); + EXPECT_EQ(status, 0); + EXPECT_EQ(attr.number, static_cast(K_FOUNDERS)); + EXPECT_EQ(attr.vl_type, static_cast(BCF_VL_FIXED)); + EXPECT_EQ(attr.type, static_cast(BCF_HT_REAL)); + + if (fid) htslib::hts_close(fid); +} + + +TEST(TestBcfHeader, BcfHdrFmtGt) { + htslib::htsFile *fid = htslib::hts_open(BCF_NAME, "r"); + bcfio::BcfHeader hdr { fid }; + + EXPECT_FALSE(hdr.isnull()); + + bcfio::BcfHdrAttr attr {}; + + int status = hdr.get_format_attr("GT", &attr); + EXPECT_EQ(status, 0); + EXPECT_EQ(attr.number, static_cast(1)); + EXPECT_EQ(attr.vl_type, static_cast(BCF_VL_FIXED)); + EXPECT_EQ(attr.type, static_cast(BCF_HT_STR)); + + if (fid) htslib::hts_close(fid); +} + + +TEST(TestBcfHeader, BcfHdrFmtGp) { + htslib::htsFile *fid = htslib::hts_open(BCF_NAME, "r"); + bcfio::BcfHeader hdr { fid }; + + EXPECT_FALSE(hdr.isnull()); + + bcfio::BcfHdrAttr attr {}; + + int status = hdr.get_format_attr("GP", &attr); + EXPECT_EQ(status, 0); + EXPECT_EQ(attr.number, static_cast(3)); + EXPECT_EQ(attr.vl_type, static_cast(BCF_VL_FIXED)); + EXPECT_EQ(attr.type, static_cast(BCF_HT_REAL)); + + if (fid) htslib::hts_close(fid); +} + +TEST(TestBcfHeader, BcfHdrFmtDs) { + htslib::htsFile *fid = htslib::hts_open(BCF_NAME, "r"); + bcfio::BcfHeader hdr { fid }; + + EXPECT_FALSE(hdr.isnull()); + + bcfio::BcfHdrAttr attr {}; + + int status = hdr.get_format_attr("DS", &attr); + EXPECT_EQ(status, 0); + EXPECT_EQ(attr.number, static_cast(1)); + EXPECT_EQ(attr.vl_type, static_cast(BCF_VL_FIXED)); + EXPECT_EQ(attr.type, static_cast(BCF_HT_REAL)); + + if (fid) htslib::hts_close(fid); +} + + +TEST(TestBcfHeader, BcfHdrFmtErr) { + htslib::htsFile *fid = htslib::hts_open(BCF_NAME, "r"); + bcfio::BcfHeader hdr { fid }; + + EXPECT_FALSE(hdr.isnull()); + + bcfio::BcfHdrAttr attr {}; + + int status = hdr.get_format_attr("DOESNOTEXIST", &attr); + EXPECT_NE(status, 0); + + if (fid) htslib::hts_close(fid); +} + + +TEST(TestBcfHeader, BcfHdrFilter) { + htslib::htsFile *fid = htslib::hts_open(BCF_NAME, "r"); + bcfio::BcfHeader hdr { fid }; + + EXPECT_FALSE(hdr.isnull()); + + bcfio::BcfHdrAttr attr {}; + + int status = hdr.get_filter_attr("PASS", &attr); + EXPECT_EQ(status, 0); + + status = hdr.get_filter_attr("PASSING", &attr); + EXPECT_NE(status, 0); + + if (fid) htslib::hts_close(fid); +} + + +TEST(TestBcfHeader, BcfHdrInfoEaf) { + htslib::htsFile *fid = htslib::hts_open(BCF_NAME, "r"); + bcfio::BcfHeader hdr { fid }; + + EXPECT_FALSE(hdr.isnull()); + + bcfio::BcfHdrAttr attr {}; + + int status = hdr.get_info_attr("EAF", &attr); + EXPECT_EQ(status, 0); + EXPECT_EQ(attr.type, static_cast(BCF_HT_REAL)); + EXPECT_EQ(attr.vl_type, static_cast(BCF_VL_VAR)); + + if (fid) htslib::hts_close(fid); +} + + +TEST(TestBcfHeader, BcfHdrInfoErc) { + htslib::htsFile *fid = htslib::hts_open(BCF_NAME, "r"); + bcfio::BcfHeader hdr { fid }; + + EXPECT_FALSE(hdr.isnull()); + + bcfio::BcfHdrAttr attr {}; + + int status = hdr.get_info_attr("ERC", &attr); + EXPECT_EQ(status, 0); + EXPECT_EQ(attr.type, static_cast(BCF_HT_REAL)); + EXPECT_EQ(attr.vl_type, static_cast(BCF_VL_VAR)); + + if (fid) htslib::hts_close(fid); +} + + +TEST(TestBcfHeader, BcfHdrInfoErr) { + htslib::htsFile *fid = htslib::hts_open(BCF_NAME, "r"); + bcfio::BcfHeader hdr { fid }; + + EXPECT_FALSE(hdr.isnull()); + + bcfio::BcfHdrAttr attr {}; + + int status = hdr.get_info_attr("NOTAINFOMEMBER", &attr); + EXPECT_NE(status, 0); + + if (fid) htslib::hts_close(fid); +} + +TEST(TestBcfHeader, BcfHdrNull) { + // htslib::htsFile *fid = htslib::hts_open("doesnotexist", "r"); + htslib::htsFile *fid = nullptr; + bcfio::BcfHeader hdr { fid }; + + EXPECT_TRUE(hdr.isnull()); + if (fid) htslib::hts_close(fid); +} + + +TEST(TestBcfHeader, Kfmt) { + htslib::htsFile *fid = htslib::hts_open(BCF_NAME, "r"); + bcfio::BcfHeader hdr { fid }; + + // DS is alt allele dosage, which is more clearly defined as the expected + // count of alt alleles under the trained HMM + EXPECT_EQ(hdr.k_fmt("DS"), 1); + EXPECT_EQ(hdr.k_fmt("HD"), K_FOUNDERS); + + // error detection + EXPECT_TRUE(hdr.k_fmt("WRONG_ID") < 0); + EXPECT_TRUE(hdr.k_fmt("") < 0); + EXPECT_TRUE(hdr.k_fmt(nullptr) < 0); +} + +TEST(TestBcfHeader, Nsamples) { + htslib::htsFile *fid = htslib::hts_open(BCF_NAME, "r"); + bcfio::BcfHeader hdr { fid }; + + // DS is alt allele dosage, which is more clearly defined as the + // expected count of alt alleles under the trained HMM + EXPECT_EQ(hdr.n_samples(), N_SAMPS); +} + +TEST(TestBcfHeader, VcfSampNames) { + htslib::htsFile *fid = htslib::hts_open(VCF_NAME, "r"); + bcfio::BcfHeader hdr { fid }; + + const std::unique_ptr s = hdr.sample_names(); + + char samp_name[] = "S011"; + for (uint8_t i = 0; i < hdr.n_samples(); i++) { + snprintf(samp_name, 5, "S%02u", i+1); + EXPECT_STREQ(s[i].c_str(), samp_name); + } +} + + +/////////////////////////////////////////////////////////////////////////// +// Test bcfio::BcfFloatRecord +/////////////////////////////////////////////////////////////////////////// + +TEST(TestBcfFloatRecord, Constructor) { + bcfio::BcfFloatRecord brec {}; + + EXPECT_EQ(brec.size(), static_cast(0)); + EXPECT_EQ(brec.get(1, 3), std::nullopt); +} + +TEST(TestBcfFloatRecord, Load) { + htslib::htsFile *fid = htslib::hts_open(BCF_NAME, "r"); + bcfio::BcfHeader hdr { fid }; + +} + + +/////////////////////////////////////////////////////////////////////////// +// Test bcfio::ReadBcf +/////////////////////////////////////////////////////////////////////////// + + +TEST(TestReadBcf, DefaultConstructor) { + bcfio::ReadBcf bcf {}; + + EXPECT_FALSE(bcf.isopen()); +} + +TEST(TestReadBcf, Constructor) { + htslib::htsFile* fid = htslib::hts_open(VCF_NAME, "r"); + bcfio::ReadBcf bcf { VCF_NAME, fid }; + + EXPECT_TRUE(bcf.isopen()); + + EXPECT_EQ(bcf.n_samples(), N_SAMPS); + EXPECT_EQ(bcf.k_fmt("HD"), K_FOUNDERS); +} + +TEST(TestReadBcf, OpenFailure) { + bcfio::ReadBcf bcf = bcfio::open("", "r"); + + EXPECT_FALSE(bcf.isopen()); +} + +TEST(TestReadBcf, OpenVCFSuccess) { + bcfio::ReadBcf bcf = bcfio::open(VCF_NAME, "r"); + EXPECT_TRUE(bcf.isopen()); +} + +TEST(TestReadBcf, OpenVCFGZSuccess) { + bcfio::ReadBcf bcf = bcfio::open(VCFGZ_NAME, "r"); + EXPECT_TRUE(bcf.isopen()); +} + +TEST(TestReadBcf, OpenBCFSuccess) { + bcfio::ReadBcf bcf = bcfio::open(BCF_NAME, "r"); + EXPECT_TRUE(bcf.isopen()); +} + +TEST(TestReadBcf, Kfmt) { + bcfio::ReadBcf bcf = bcfio::open( VCF_NAME, "r"); + + // DS is alt allele dosage, which is more clearly defined as the expected + // count of alt alleles under the trained HMM + EXPECT_EQ(bcf.k_fmt("DS"), 1); + EXPECT_EQ(bcf.k_fmt("HD"), K_FOUNDERS); + + // TODO: what happens if I submit "GT", it exists but is a string + // not float + // error detection + EXPECT_TRUE(bcf.k_fmt("WRONG_ID") < 0); + EXPECT_TRUE(bcf.k_fmt("") < 0); + EXPECT_TRUE(bcf.k_fmt(nullptr) < 0); +} + +TEST(TestReadBcf, VcfGzSampNames) { + bcfio::ReadBcf bcf = bcfio::open( VCFGZ_NAME, "r"); + + const std::unique_ptr s = bcf.sample_names(); + + char samp_name[] = "S011"; + + for (uint8_t i = 0; i < bcf.n_samples(); i++) { + snprintf(samp_name, 5, "S%02u", i+1); + EXPECT_STREQ(s[i].c_str(), samp_name); + } +} + + +TEST(TestReadBcf, BcfSampNames) { + bcfio::ReadBcf bcf = bcfio::open( BCF_NAME, "r"); + + const std::unique_ptr s = bcf.sample_names(); + + char samp_name[] = "S011"; + + for (uint8_t i = 0; i < bcf.n_samples(); i++) { + snprintf(samp_name, 5, "S%02u", i+1); + EXPECT_STREQ(s[i].c_str(), samp_name); + } +} + + +TEST(TestReadBcf, LoadRecord) { + + bcfio::ReadBcf bcf = bcfio::open(VCF_NAME, "r"); + bcfio::BcfFloatRecord rec {}; + + bcf.next_record(&rec, "HD"); + + int32_t k_founders = bcf.k_fmt("HD"); + EXPECT_FALSE(k_founders <= 0); + EXPECT_EQ(rec.size(), bcf.n_samples() * static_cast(k_founders)); + EXPECT_EQ(bcf.n_samples(), rec.nrows()); + EXPECT_EQ(static_cast(k_founders), rec.ncols()); + EXPECT_TRUE(rec.is_snp()); + + // EXPECT_EQ(record.chrom(), "chr12"); + // EXPECT_EQ(record.pos(), 788); + // EXPECT_EQ(record.id(), "."); + // EXPECT_EQ(record.ref(), 'A'); + // EXPECT_EQ(record.alt(), 'G'); + // EXkECT_EQ(record.qual(), "."); + // EXPECT_EQ(record.filter(), "PASS"); + // EXPECT_EQ(record.info(), "EAF=0.00228;INFO_SCORE=1;HWE=1;ERC=0.01949;EAC=7.94153;PAF=0.00245;REF_PANEL=0"); + // EXPECT_EQ(record.format(), "GT:GP:DS:HD"); +} + + + +struct Buff { + Buff(const uint32_t size_in) + : size(size_in), array(new char[size]) { reset(); }; + ~Buff() { if (array) delete[] array; }; + + void reset() { + std::memset(array, '\0', size); + } + + uint32_t size; + char* array; +}; + +struct FloatArray { + FloatArray(const uint32_t size_in) + : size(size_in), array(new float[size]) { reset(); }; + ~FloatArray() { if (array) delete[] array; }; + + void reset() { + for (uint32_t i = 0; i < size; i++) + array[i] = static_cast(0); + } + + uint32_t size; + float *array; +}; + + +int get_sample_truth_vals(FILE* fid, + FloatArray* data, + Buff* buff) { + + size_t buff_idx = 0; + size_t data_idx = 0; + int c; + while ((c = fgetc(fid)) != EOF) { + + if (c == ',' || c == '\n') { + if (buff_idx >= buff->size-1) + return -1; + buff->array[buff_idx] = '\0'; + + if (data_idx >= data->size) + return -1; + + data->array[data_idx++] = atof(buff->array); + buff_idx = 0; + buff->reset(); + + if (c == '\n') break; + + continue; + } + + buff->array[buff_idx++] = c; + } + + return 0; +} + + + +TEST(TestReadBcf, HDRecordValue) { + + bcfio::ReadBcf bcf = bcfio::open(VCF_NAME, "r"); + bcfio::BcfFloatRecord rec {}; + + Buff buff_fname { 100 }; + Buff buff_data { 1000 }; + FloatArray data { static_cast(bcf.k_fmt("HD")) }; + + FILE* fid; + // iterate positions + size_t pos = 1; + int status; + while (bcf.next_record(&rec, "HD") == 0) { + + snprintf(buff_fname.array, + buff_fname.size, + "tests/hd_%02zu.csv", pos++); + + fid = fopen(buff_fname.array, "r"); + + EXPECT_EQ(rec.nrows(), bcf.n_samples()); + EXPECT_EQ(rec.ncols(), static_cast(bcf.k_fmt("HD"))); + + // loop over samples + for (uint64_t i = 0; i < rec.nrows(); i++) { + + status = get_sample_truth_vals(fid, &data, &buff_data); + if (status != 0) + printf("\n\nERROR\n\n"); + + // loop over haplotypes + for (uint64_t j = 0; j < rec.ncols(); j++) + EXPECT_EQ(data.array[j], rec.get(i, j).value()); + } + fclose(fid); + } +} diff --git a/tests/test_grm.cpp b/tests/test_grm.cpp new file mode 100644 index 0000000..7917529 --- /dev/null +++ b/tests/test_grm.cpp @@ -0,0 +1,1208 @@ + +#include +#include +#include +#include +#include + + +//////////////////////////////////////////////////////////////////// +// COORDINATES TESTS +//////////////////////////////////////////////////////////////////// + +TEST(TestCoords, DefaultConstructor) { + // verify default values + grm::Coordinates coords {}; + + std::string contig = std::string(""); + + EXPECT_EQ(coords.contig.size(), static_cast(0)); + EXPECT_EQ(coords.contig, contig); + + EXPECT_EQ(coords.len, static_cast(0)); + EXPECT_EQ(coords.pos, nullptr); +} + + +TEST(TestCoords, ConstructorValidInput) { + char contig_in[] = "chr12"; + std::string contig { contig_in }; + uint64_t len_in = 1000; + + grm::Coordinates coords { contig_in, len_in }; + + EXPECT_EQ(coords.contig, contig); + EXPECT_EQ(coords.len, len_in); + EXPECT_NE(coords.pos, nullptr); +} + +// If input contig name is nullptr, resort to the default values +// of the default constructor. +TEST(TestCoords, ConstructorInvalidInput) { + uint64_t len_in = 10; + grm::Coordinates coords { nullptr, len_in }; + + EXPECT_EQ(coords.contig, std::string("")); + EXPECT_EQ(coords.len, static_cast(0)); + EXPECT_EQ(coords.pos, nullptr); +} + + +TEST(TestCoords, ConstructorZeroLength) { + char contig_in[] = "chr1"; + grm::Coordinates coords { contig_in, 0 }; + + EXPECT_EQ(coords.contig, std::string("chr1")); + EXPECT_EQ(coords.len, static_cast(0)); + EXPECT_EQ(coords.pos, nullptr); +} + + +TEST(TestCoords, MoveConstructor) { + char contig_in[] = "chr7"; + uint64_t len_in = 5; + grm::Coordinates src { contig_in, len_in }; + + // set positions to known values + for (uint64_t i = 0; i < len_in; i++) + src.pos[i] = (i + 1) * 100; + + grm::Coordinates dst { std::move(src) }; + + // destination should have source's data + EXPECT_EQ(dst.contig, std::string("chr7")); + EXPECT_EQ(dst.len, len_in); + EXPECT_NE(dst.pos, nullptr); + for (uint64_t i = 0; i < len_in; i++) + EXPECT_EQ(dst.pos[i], (i + 1) * 100); + + // source should be in moved-from state + EXPECT_EQ(src.contig, std::string("")); + EXPECT_EQ(src.len, static_cast(0)); + EXPECT_EQ(src.pos, nullptr); +} + + +TEST(TestCoords, MoveAssignment) { + char contig_in[] = "chr3"; + uint64_t len_in = 3; + grm::Coordinates src { contig_in, len_in }; + src.pos[0] = 10; + src.pos[1] = 20; + src.pos[2] = 30; + + grm::Coordinates dst {}; + dst = std::move(src); + + EXPECT_EQ(dst.contig, std::string("chr3")); + EXPECT_EQ(dst.len, static_cast(3)); + EXPECT_NE(dst.pos, nullptr); + EXPECT_EQ(dst.pos[0], static_cast(10)); + EXPECT_EQ(dst.pos[1], static_cast(20)); + EXPECT_EQ(dst.pos[2], static_cast(30)); + + EXPECT_EQ(src.contig, std::string("")); + EXPECT_EQ(src.len, static_cast(0)); + EXPECT_EQ(src.pos, nullptr); +} + + +TEST(TestCoords, WriteReadRoundTrip) { + char contig_in[] = "chr22"; + uint64_t len_in = 4; + grm::Coordinates src { contig_in, len_in }; + src.pos[0] = 100; + src.pos[1] = 200; + src.pos[2] = 500; + src.pos[3] = 1000; + + // write to a temporary file + io::FileIO fio_w { tmpfile() }; + ASSERT_NE(fio_w.fid, nullptr); + + grm::STATUS status = grm::write(&fio_w, &src); + ASSERT_EQ(status, grm::SUCCESS); + + // rewind and read back + rewind(fio_w.fid); + grm::Coordinates dst {}; + status = grm::read(&fio_w, &dst); + ASSERT_EQ(status, grm::SUCCESS); + + EXPECT_EQ(dst.contig, std::string("chr22")); + EXPECT_EQ(dst.len, len_in); + ASSERT_NE(dst.pos, nullptr); + for (uint64_t i = 0; i < len_in; i++) + EXPECT_EQ(dst.pos[i], src.pos[i]); +} + + +TEST(TestCoords, WriteNullArgs) { + grm::Coordinates coords {}; + io::FileIO fio { tmpfile() }; + ASSERT_NE(fio.fid, nullptr); + + EXPECT_EQ(grm::write(nullptr, &coords), grm::ERROR_NULLPTR_ARG); + EXPECT_EQ(grm::write(&fio, static_cast(nullptr)), + grm::ERROR_NULLPTR_ARG); +} + + +TEST(TestCoords, ReadNullArgs) { + grm::Coordinates coords {}; + io::FileIO fio { tmpfile() }; + ASSERT_NE(fio.fid, nullptr); + + EXPECT_EQ(grm::read(nullptr, &coords), grm::ERROR_NULLPTR_ARG); + EXPECT_EQ(grm::read(&fio, static_cast(nullptr)), + grm::ERROR_NULLPTR_ARG); +} + + +TEST(TestCoords, WriteReadSinglePosition) { + char contig_in[] = "chrX"; + grm::Coordinates src { contig_in, 1 }; + src.pos[0] = 42; + + io::FileIO fio { tmpfile() }; + ASSERT_NE(fio.fid, nullptr); + + ASSERT_EQ(grm::write(&fio, &src), grm::SUCCESS); + rewind(fio.fid); + + grm::Coordinates dst {}; + ASSERT_EQ(grm::read(&fio, &dst), grm::SUCCESS); + + EXPECT_EQ(dst.contig, "chrX"); + EXPECT_EQ(dst.len, static_cast(1)); + ASSERT_NE(dst.pos, nullptr); + EXPECT_EQ(dst.pos[0], static_cast(42)); +} + + +// A FileIO with a null fid should be treated like a null arg +TEST(TestCoords, WriteNullFid) { + grm::Coordinates coords {}; + io::FileIO fio { nullptr }; + + EXPECT_EQ(grm::write(&fio, &coords), grm::ERROR_NULLPTR_ARG); +} + + +TEST(TestCoords, ReadNullFid) { + grm::Coordinates coords {}; + io::FileIO fio { nullptr }; + + EXPECT_EQ(grm::read(&fio, &coords), grm::ERROR_NULLPTR_ARG); +} + + +// Overwriting a previously populated Coordinates with read +TEST(TestCoords, ReadOverwritesExisting) { + char contig_in[] = "chr9"; + grm::Coordinates src { contig_in, 2 }; + src.pos[0] = 10; + src.pos[1] = 20; + + io::FileIO fio { tmpfile() }; + ASSERT_NE(fio.fid, nullptr); + ASSERT_EQ(grm::write(&fio, &src), grm::SUCCESS); + rewind(fio.fid); + + // dst starts with different data + char other_contig[] = "chr1"; + grm::Coordinates dst { other_contig, 5 }; + ASSERT_EQ(grm::read(&fio, &dst), grm::SUCCESS); + + EXPECT_EQ(dst.contig, "chr9"); + EXPECT_EQ(dst.len, static_cast(2)); + EXPECT_EQ(dst.pos[0], static_cast(10)); + EXPECT_EQ(dst.pos[1], static_cast(20)); +} + + +//////////////////////////////////////////////////////////////////// +// SAMPLES TESTS +//////////////////////////////////////////////////////////////////// + +TEST(TestSamples, DefaultConstructor) { + grm::Samples samps {}; + + EXPECT_EQ(samps.len, static_cast(0)); + EXPECT_EQ(samps.names, nullptr); +} + + +TEST(TestSamples, ConstructorValidInput) { + uint64_t n = 5; + grm::Samples samps { n }; + + EXPECT_EQ(samps.len, n); + EXPECT_NE(samps.names, nullptr); +} + + +TEST(TestSamples, ConstructorZero) { + grm::Samples samps { 0 }; + + EXPECT_EQ(samps.len, static_cast(0)); + EXPECT_EQ(samps.names, nullptr); +} + + +TEST(TestSamples, MoveConstructor) { + grm::Samples src { 3 }; + src.names[0] = "sample_A"; + src.names[1] = "sample_B"; + src.names[2] = "sample_C"; + + grm::Samples dst { std::move(src) }; + + EXPECT_EQ(dst.len, static_cast(3)); + ASSERT_NE(dst.names, nullptr); + EXPECT_EQ(dst.names[0], "sample_A"); + EXPECT_EQ(dst.names[1], "sample_B"); + EXPECT_EQ(dst.names[2], "sample_C"); + + EXPECT_EQ(src.len, static_cast(0)); + EXPECT_EQ(src.names, nullptr); +} + + +TEST(TestSamples, MoveAssignment) { + grm::Samples src { 2 }; + src.names[0] = "id_1"; + src.names[1] = "id_2"; + + grm::Samples dst {}; + dst = std::move(src); + + EXPECT_EQ(dst.len, static_cast(2)); + ASSERT_NE(dst.names, nullptr); + EXPECT_EQ(dst.names[0], "id_1"); + EXPECT_EQ(dst.names[1], "id_2"); + + EXPECT_EQ(src.len, static_cast(0)); + EXPECT_EQ(src.names, nullptr); +} + + +TEST(TestSamples, WriteReadRoundTrip) { + grm::Samples src { 3 }; + src.names[0] = "alpha"; + src.names[1] = "beta"; + src.names[2] = "gamma"; + + io::FileIO fio { tmpfile() }; + ASSERT_NE(fio.fid, nullptr); + + grm::STATUS status = grm::write(&fio, &src); + ASSERT_EQ(status, grm::SUCCESS); + + rewind(fio.fid); + + grm::Samples dst {}; + status = grm::read(&fio, &dst); + ASSERT_EQ(status, grm::SUCCESS); + + EXPECT_EQ(dst.len, static_cast(3)); + ASSERT_NE(dst.names, nullptr); + EXPECT_EQ(dst.names[0], "alpha"); + EXPECT_EQ(dst.names[1], "beta"); + EXPECT_EQ(dst.names[2], "gamma"); +} + + +TEST(TestSamples, WriteReadVaryingLengthNames) { + grm::Samples src { 3 }; + src.names[0] = "a"; + src.names[1] = "longer_sample_name"; + src.names[2] = "xy"; + + io::FileIO fio { tmpfile() }; + ASSERT_NE(fio.fid, nullptr); + + ASSERT_EQ(grm::write(&fio, &src), grm::SUCCESS); + rewind(fio.fid); + + grm::Samples dst {}; + ASSERT_EQ(grm::read(&fio, &dst), grm::SUCCESS); + + EXPECT_EQ(dst.len, static_cast(3)); + ASSERT_NE(dst.names, nullptr); + EXPECT_EQ(dst.names[0], "a"); + EXPECT_EQ(dst.names[1], "longer_sample_name"); + EXPECT_EQ(dst.names[2], "xy"); +} + + +TEST(TestSamples, WriteNullArgs) { + grm::Samples samps {}; + io::FileIO fio { tmpfile() }; + ASSERT_NE(fio.fid, nullptr); + + EXPECT_EQ(grm::write(nullptr, &samps), grm::ERROR_NULLPTR_ARG); + EXPECT_EQ(grm::write(&fio, static_cast(nullptr)), + grm::ERROR_NULLPTR_ARG); +} + + +TEST(TestSamples, ReadNullArgs) { + grm::Samples samps {}; + io::FileIO fio { tmpfile() }; + ASSERT_NE(fio.fid, nullptr); + + EXPECT_EQ(grm::read(nullptr, &samps), grm::ERROR_NULLPTR_ARG); + EXPECT_EQ(grm::read(&fio, static_cast(nullptr)), + grm::ERROR_NULLPTR_ARG); +} + + +TEST(TestSamples, WriteReadSingleSample) { + grm::Samples src { 1 }; + src.names[0] = "only_sample"; + + io::FileIO fio { tmpfile() }; + ASSERT_NE(fio.fid, nullptr); + + ASSERT_EQ(grm::write(&fio, &src), grm::SUCCESS); + rewind(fio.fid); + + grm::Samples dst {}; + ASSERT_EQ(grm::read(&fio, &dst), grm::SUCCESS); + + EXPECT_EQ(dst.len, static_cast(1)); + ASSERT_NE(dst.names, nullptr); + EXPECT_EQ(dst.names[0], "only_sample"); +} + + +// Writing samples where all names are empty strings should fail +// because nchar_max would be 0, triggering ERROR_INVALID_ARG +TEST(TestSamples, WriteEmptyNamesReturnsError) { + grm::Samples samps { 2 }; + // names are default-constructed empty strings + + io::FileIO fio { tmpfile() }; + ASSERT_NE(fio.fid, nullptr); + + EXPECT_EQ(grm::write(&fio, &samps), grm::ERROR_INVALID_ARG); +} + + +TEST(TestSamples, WriteNullFid) { + grm::Samples samps { 1 }; + samps.names[0] = "s1"; + io::FileIO fio { nullptr }; + + EXPECT_EQ(grm::write(&fio, &samps), grm::ERROR_NULLPTR_ARG); +} + + +TEST(TestSamples, ReadNullFid) { + grm::Samples samps {}; + io::FileIO fio { nullptr }; + + EXPECT_EQ(grm::read(&fio, &samps), grm::ERROR_NULLPTR_ARG); +} + + +// Overwriting a previously populated Samples with read +TEST(TestSamples, ReadOverwritesExisting) { + grm::Samples src { 2 }; + src.names[0] = "new_a"; + src.names[1] = "new_b"; + + io::FileIO fio { tmpfile() }; + ASSERT_NE(fio.fid, nullptr); + ASSERT_EQ(grm::write(&fio, &src), grm::SUCCESS); + rewind(fio.fid); + + // dst starts with different data + grm::Samples dst { 3 }; + dst.names[0] = "old_x"; + dst.names[1] = "old_y"; + dst.names[2] = "old_z"; + + ASSERT_EQ(grm::read(&fio, &dst), grm::SUCCESS); + + EXPECT_EQ(dst.len, static_cast(2)); + ASSERT_NE(dst.names, nullptr); + EXPECT_EQ(dst.names[0], "new_a"); + EXPECT_EQ(dst.names[1], "new_b"); +} + + +//////////////////////////////////////////////////////////////////// +// HDR TESTS +//////////////////////////////////////////////////////////////////// + +TEST(TestHdr, DefaultConstructor) { + grm::Hdr hdr {}; + + EXPECT_EQ(hdr.prog_version.major, constants::PROG_VERSION.major); + EXPECT_EQ(hdr.prog_version.minor, constants::PROG_VERSION.minor); + EXPECT_EQ(hdr.prog_version.micro, constants::PROG_VERSION.micro); + + EXPECT_EQ(hdr.file_version.major, grm::FILE_VERSION.major); + EXPECT_EQ(hdr.file_version.minor, grm::FILE_VERSION.minor); + EXPECT_EQ(hdr.file_version.micro, grm::FILE_VERSION.micro); + + EXPECT_EQ(hdr.grm_type, grm::UNSPECIFIED); + EXPECT_NE(hdr.coords, nullptr); + EXPECT_NE(hdr.samples, nullptr); +} + + +TEST(TestHdr, MoveConstructor) { + grm::Hdr src {}; + src.grm_type = grm::EHC; + src.coords->contig = "chr1"; + + utils::Version saved_prog = src.prog_version; + utils::Version saved_file = src.file_version; + + grm::Hdr dst { std::move(src) }; + + EXPECT_EQ(dst.grm_type, grm::EHC); + EXPECT_EQ(dst.prog_version.major, saved_prog.major); + EXPECT_EQ(dst.file_version.major, saved_file.major); + EXPECT_NE(dst.coords, nullptr); + EXPECT_EQ(dst.coords->contig, "chr1"); + + // source should be reset + EXPECT_EQ(src.grm_type, grm::UNSPECIFIED); + EXPECT_EQ(src.coords, nullptr); + EXPECT_EQ(src.samples, nullptr); +} + + +TEST(TestHdr, MoveAssignment) { + grm::Hdr src {}; + src.grm_type = grm::DS; + + grm::Hdr dst {}; + dst = std::move(src); + + EXPECT_EQ(dst.grm_type, grm::DS); + EXPECT_NE(dst.coords, nullptr); + EXPECT_NE(dst.samples, nullptr); + + EXPECT_EQ(src.grm_type, grm::UNSPECIFIED); + EXPECT_EQ(src.coords, nullptr); + EXPECT_EQ(src.samples, nullptr); +} + + +TEST(TestHdr, WriteReadRoundTrip) { + grm::Hdr src {}; + src.grm_type = grm::EAC; + + // set up coordinates + char contig[] = "chr5"; + *src.coords = grm::Coordinates { contig, 3 }; + src.coords->pos[0] = 100; + src.coords->pos[1] = 200; + src.coords->pos[2] = 300; + + // set up samples + *src.samples = grm::Samples { 2 }; + src.samples->names[0] = "samp1"; + src.samples->names[1] = "samp2"; + + io::FileIO fio { tmpfile() }; + ASSERT_NE(fio.fid, nullptr); + + ASSERT_EQ(grm::write(&fio, &src), grm::SUCCESS); + + rewind(fio.fid); + + grm::Hdr dst {}; + ASSERT_EQ(grm::read(&fio, &dst), grm::SUCCESS); + + EXPECT_EQ(dst.prog_version.major, src.prog_version.major); + EXPECT_EQ(dst.prog_version.minor, src.prog_version.minor); + EXPECT_EQ(dst.prog_version.micro, src.prog_version.micro); + EXPECT_EQ(dst.file_version.major, src.file_version.major); + EXPECT_EQ(dst.grm_type, grm::EAC); + + ASSERT_NE(dst.coords, nullptr); + EXPECT_EQ(dst.coords->contig, "chr5"); + EXPECT_EQ(dst.coords->len, static_cast(3)); + EXPECT_EQ(dst.coords->pos[0], static_cast(100)); + EXPECT_EQ(dst.coords->pos[1], static_cast(200)); + EXPECT_EQ(dst.coords->pos[2], static_cast(300)); + + ASSERT_NE(dst.samples, nullptr); + EXPECT_EQ(dst.samples->len, static_cast(2)); + EXPECT_EQ(dst.samples->names[0], "samp1"); + EXPECT_EQ(dst.samples->names[1], "samp2"); +} + + +TEST(TestHdr, WriteNullArgs) { + grm::Hdr hdr {}; + io::FileIO fio { tmpfile() }; + ASSERT_NE(fio.fid, nullptr); + + EXPECT_EQ(grm::write(nullptr, &hdr), grm::ERROR_NULLPTR_ARG); + EXPECT_EQ(grm::write(&fio, static_cast(nullptr)), + grm::ERROR_NULLPTR_ARG); +} + + +TEST(TestHdr, ReadNullArgs) { + grm::Hdr hdr {}; + io::FileIO fio { tmpfile() }; + ASSERT_NE(fio.fid, nullptr); + + EXPECT_EQ(grm::read(nullptr, &hdr), grm::ERROR_NULLPTR_ARG); + EXPECT_EQ(grm::read(&fio, static_cast(nullptr)), + grm::ERROR_NULLPTR_ARG); +} + + +// Verify that every GrmType enum value survives a write/read round-trip +TEST(TestHdr, WriteReadAllGrmTypes) { + grm::GrmType types[] = { + grm::EHC, grm::EAC, grm::BOTH, grm::DS, grm::UNSPECIFIED + }; + + for (grm::GrmType t : types) { + grm::Hdr src {}; + src.grm_type = t; + + // need valid samples for write to succeed + *src.samples = grm::Samples { 1 }; + src.samples->names[0] = "s"; + + io::FileIO fio { tmpfile() }; + ASSERT_NE(fio.fid, nullptr); + ASSERT_EQ(grm::write(&fio, &src), grm::SUCCESS); + + rewind(fio.fid); + + grm::Hdr dst {}; + ASSERT_EQ(grm::read(&fio, &dst), grm::SUCCESS); + EXPECT_EQ(dst.grm_type, t); + } +} + + +// Verify version fields survive pack/unpack through file I/O +TEST(TestHdr, VersionRoundTrip) { + grm::Hdr src {}; + // prog_version and file_version are set by defaults; + // verify they round-trip exactly + + *src.samples = grm::Samples { 1 }; + src.samples->names[0] = "s"; + + io::FileIO fio { tmpfile() }; + ASSERT_NE(fio.fid, nullptr); + ASSERT_EQ(grm::write(&fio, &src), grm::SUCCESS); + rewind(fio.fid); + + grm::Hdr dst {}; + ASSERT_EQ(grm::read(&fio, &dst), grm::SUCCESS); + + EXPECT_EQ(dst.prog_version.major, constants::PROG_VERSION.major); + EXPECT_EQ(dst.prog_version.minor, constants::PROG_VERSION.minor); + EXPECT_EQ(dst.prog_version.micro, constants::PROG_VERSION.micro); + + EXPECT_EQ(dst.file_version.major, grm::FILE_VERSION.major); + EXPECT_EQ(dst.file_version.minor, grm::FILE_VERSION.minor); + EXPECT_EQ(dst.file_version.micro, grm::FILE_VERSION.micro); +} + + +TEST(TestHdr, WriteNullFid) { + grm::Hdr hdr {}; + io::FileIO fio { nullptr }; + + EXPECT_EQ(grm::write(&fio, &hdr), grm::ERROR_NULLPTR_ARG); +} + + +TEST(TestHdr, ReadNullFid) { + grm::Hdr hdr {}; + io::FileIO fio { nullptr }; + + EXPECT_EQ(grm::read(&fio, &hdr), grm::ERROR_NULLPTR_ARG); +} + + +//////////////////////////////////////////////////////////////////// +// GRM STRUCT TESTS +//////////////////////////////////////////////////////////////////// + +TEST(TestGrm, DefaultConstructor) { + grm::Grm g {}; + + EXPECT_EQ(g.n_samples, static_cast(0)); + EXPECT_EQ(g.data, nullptr); + EXPECT_EQ(g.size(), static_cast(0)); +} + + +TEST(TestGrm, ConstructorValidInput) { + grm::Grm g { 4 }; + + EXPECT_EQ(g.n_samples, static_cast(4)); + EXPECT_NE(g.data, nullptr); + EXPECT_EQ(g.size(), static_cast(4 * 5 / 2)); // n*(n+1)/2 = 10 + + // data should be zero-initialized + for (uint64_t i = 0; i < g.size(); i++) + EXPECT_FLOAT_EQ(g.data[i], 0.0f); +} + + +TEST(TestGrm, ConstructorZero) { + grm::Grm g { 0 }; + + EXPECT_EQ(g.n_samples, static_cast(0)); + EXPECT_EQ(g.data, nullptr); + EXPECT_EQ(g.size(), static_cast(0)); +} + + +TEST(TestGrm, Size) { + EXPECT_EQ(grm::Grm(0).size(), static_cast(0)); + EXPECT_EQ(grm::Grm(1).size(), static_cast(1)); + EXPECT_EQ(grm::Grm(2).size(), static_cast(3)); + EXPECT_EQ(grm::Grm(3).size(), static_cast(6)); + EXPECT_EQ(grm::Grm(4).size(), static_cast(10)); + EXPECT_EQ(grm::Grm(5).size(), static_cast(15)); +} + + +TEST(TestGrm, MidxToArr) { + // Verify the manual example from grm.h comments for n=3 + grm::Grm g { 3 }; + uint64_t idx = 0; + + // Upper triangle and diagonal + EXPECT_EQ(g.midx_to_arr(0, 0, &idx), grm::SUCCESS); + EXPECT_EQ(idx, static_cast(0)); + + EXPECT_EQ(g.midx_to_arr(0, 1, &idx), grm::SUCCESS); + EXPECT_EQ(idx, static_cast(1)); + + EXPECT_EQ(g.midx_to_arr(0, 2, &idx), grm::SUCCESS); + EXPECT_EQ(idx, static_cast(2)); + + EXPECT_EQ(g.midx_to_arr(1, 1, &idx), grm::SUCCESS); + EXPECT_EQ(idx, static_cast(3)); + + EXPECT_EQ(g.midx_to_arr(1, 2, &idx), grm::SUCCESS); + EXPECT_EQ(idx, static_cast(4)); + + EXPECT_EQ(g.midx_to_arr(2, 2, &idx), grm::SUCCESS); + EXPECT_EQ(idx, static_cast(5)); + + // Lower triangle should map to same idx by symmetry + EXPECT_EQ(g.midx_to_arr(1, 0, &idx), grm::SUCCESS); + EXPECT_EQ(idx, static_cast(1)); // same as (0,1) + + EXPECT_EQ(g.midx_to_arr(2, 0, &idx), grm::SUCCESS); + EXPECT_EQ(idx, static_cast(2)); // same as (0,2) + + EXPECT_EQ(g.midx_to_arr(2, 1, &idx), grm::SUCCESS); + EXPECT_EQ(idx, static_cast(4)); // same as (1,2) +} + + +TEST(TestGrm, MidxToArrBoundsCheck) { + grm::Grm g { 3 }; + uint64_t idx = 0; + + EXPECT_EQ(g.midx_to_arr(3, 0, &idx), grm::ERROR_IDX_ARR_BOUNDS); + EXPECT_EQ(g.midx_to_arr(0, 3, &idx), grm::ERROR_IDX_ARR_BOUNDS); + EXPECT_EQ(g.midx_to_arr(3, 3, &idx), grm::ERROR_IDX_ARR_BOUNDS); +} + + +TEST(TestGrm, OperatorParensSetAndGet) { + grm::Grm g { 3 }; + + // Set via operator() + g(0, 0) = 1.0f; + g(0, 1) = 2.0f; + g(0, 2) = 3.0f; + g(1, 1) = 4.0f; + g(1, 2) = 5.0f; + g(2, 2) = 6.0f; + + // Read back via operator() const + const grm::Grm& cg = g; + EXPECT_FLOAT_EQ(cg(0, 0), 1.0f); + EXPECT_FLOAT_EQ(cg(0, 1), 2.0f); + EXPECT_FLOAT_EQ(cg(0, 2), 3.0f); + EXPECT_FLOAT_EQ(cg(1, 1), 4.0f); + EXPECT_FLOAT_EQ(cg(1, 2), 5.0f); + EXPECT_FLOAT_EQ(cg(2, 2), 6.0f); +} + + +TEST(TestGrm, OperatorParensSymmetry) { + grm::Grm g { 3 }; + + g(0, 1) = 7.5f; + g(2, 0) = 3.3f; + + const grm::Grm& cg = g; + + // (i,j) should equal (j,i) due to symmetry + EXPECT_FLOAT_EQ(cg(0, 1), cg(1, 0)); + EXPECT_FLOAT_EQ(cg(0, 2), cg(2, 0)); + EXPECT_FLOAT_EQ(cg(0, 1), 7.5f); + EXPECT_FLOAT_EQ(cg(0, 2), 3.3f); +} + + +TEST(TestGrm, SetAndGet) { + grm::Grm g { 3 }; + + EXPECT_EQ(g.set(0, 0, 1.1f), grm::SUCCESS); + EXPECT_EQ(g.set(1, 2, 2.2f), grm::SUCCESS); + + float val = 0.0f; + EXPECT_EQ(g.get(0, 0, &val), grm::SUCCESS); + EXPECT_FLOAT_EQ(val, 1.1f); + + EXPECT_EQ(g.get(1, 2, &val), grm::SUCCESS); + EXPECT_FLOAT_EQ(val, 2.2f); + + // symmetry + EXPECT_EQ(g.get(2, 1, &val), grm::SUCCESS); + EXPECT_FLOAT_EQ(val, 2.2f); +} + + +TEST(TestGrm, SetGetBoundsCheck) { + grm::Grm g { 3 }; + + EXPECT_EQ(g.set(3, 0, 1.0f), grm::ERROR_IDX_ARR_BOUNDS); + EXPECT_EQ(g.set(0, 3, 1.0f), grm::ERROR_IDX_ARR_BOUNDS); + + float val = 0.0f; + EXPECT_EQ(g.get(3, 0, &val), grm::ERROR_IDX_ARR_BOUNDS); + EXPECT_EQ(g.get(0, 3, &val), grm::ERROR_IDX_ARR_BOUNDS); +} + + +TEST(TestGrm, MoveConstructor) { + grm::Grm src { 3 }; + src(0, 0) = 1.0f; + src(1, 2) = 5.0f; + + grm::Grm dst { std::move(src) }; + + EXPECT_EQ(dst.n_samples, static_cast(3)); + EXPECT_NE(dst.data, nullptr); + EXPECT_FLOAT_EQ(dst(0, 0), 1.0f); + EXPECT_FLOAT_EQ(dst(1, 2), 5.0f); + + EXPECT_EQ(src.n_samples, static_cast(0)); + EXPECT_EQ(src.data, nullptr); +} + + +TEST(TestGrm, MoveAssignment) { + grm::Grm src { 2 }; + src(0, 0) = 1.0f; + src(0, 1) = 2.0f; + src(1, 1) = 3.0f; + + grm::Grm dst {}; + dst = std::move(src); + + EXPECT_EQ(dst.n_samples, static_cast(2)); + EXPECT_FLOAT_EQ(dst(0, 0), 1.0f); + EXPECT_FLOAT_EQ(dst(0, 1), 2.0f); + EXPECT_FLOAT_EQ(dst(1, 1), 3.0f); + + EXPECT_EQ(src.n_samples, static_cast(0)); + EXPECT_EQ(src.data, nullptr); +} + + +TEST(TestGrm, SingleSampleMatrix) { + grm::Grm g { 1 }; + + EXPECT_EQ(g.n_samples, static_cast(1)); + EXPECT_EQ(g.size(), static_cast(1)); + EXPECT_NE(g.data, nullptr); + + g(0, 0) = 2.5f; + EXPECT_FLOAT_EQ(g(0, 0), 2.5f); + + float val = 0.0f; + EXPECT_EQ(g.get(0, 0, &val), grm::SUCCESS); + EXPECT_FLOAT_EQ(val, 2.5f); + + EXPECT_EQ(g.set(0, 0, 3.0f), grm::SUCCESS); + EXPECT_EQ(g.get(0, 0, &val), grm::SUCCESS); + EXPECT_FLOAT_EQ(val, 3.0f); +} + + +// Set values via the lower triangle (i > j), confirm they are +// accessible via the upper triangle (i < j) +TEST(TestGrm, SetViaLowerTriangleGetViaUpper) { + grm::Grm g { 4 }; + + // set off-diagonal values using lower triangle indices + g.set(1, 0, 1.1f); + g.set(2, 0, 2.2f); + g.set(2, 1, 3.3f); + g.set(3, 0, 4.4f); + g.set(3, 1, 5.5f); + g.set(3, 2, 6.6f); + + // get via upper triangle + float val = 0.0f; + g.get(0, 1, &val); EXPECT_FLOAT_EQ(val, 1.1f); + g.get(0, 2, &val); EXPECT_FLOAT_EQ(val, 2.2f); + g.get(1, 2, &val); EXPECT_FLOAT_EQ(val, 3.3f); + g.get(0, 3, &val); EXPECT_FLOAT_EQ(val, 4.4f); + g.get(1, 3, &val); EXPECT_FLOAT_EQ(val, 5.5f); + g.get(2, 3, &val); EXPECT_FLOAT_EQ(val, 6.6f); +} + + +// Validate indexing for a 5x5 matrix +TEST(TestGrm, LargerMatrixIndexing) { + uint64_t n = 5; + grm::Grm g { n }; + + // fill every element with a unique value + float counter = 1.0f; + for (uint64_t i = 0; i < n; i++) + for (uint64_t j = i; j < n; j++) + g.set(i, j, counter++); + + // verify all values + counter = 1.0f; + float val = 0.0f; + for (uint64_t i = 0; i < n; i++) { + for (uint64_t j = i; j < n; j++) { + g.get(i, j, &val); + EXPECT_FLOAT_EQ(val, counter++); + } + } +} + + +// Verify that move assignment replaces a non-empty Grm +TEST(TestGrm, MoveAssignmentReplacesExisting) { + grm::Grm dst { 2 }; + dst(0, 0) = 99.0f; + + grm::Grm src { 4 }; + src(0, 0) = 1.0f; + src(3, 3) = 42.0f; + + dst = std::move(src); + + EXPECT_EQ(dst.n_samples, static_cast(4)); + EXPECT_FLOAT_EQ(dst(0, 0), 1.0f); + EXPECT_FLOAT_EQ(dst(3, 3), 42.0f); +} + + +//////////////////////////////////////////////////////////////////// +// MACRO TEST +//////////////////////////////////////////////////////////////////// + +TEST(TestMatrixMacro, ManualValidation) { + // Validate MATRIX_IDX_TO_ARRAY against the worked example + // in grm.h for a 3x3 matrix + uint64_t n = 3; + EXPECT_EQ(MATRIX_IDX_TO_ARRAY(0, 0, n), static_cast(0)); + EXPECT_EQ(MATRIX_IDX_TO_ARRAY(0, 1, n), static_cast(1)); + EXPECT_EQ(MATRIX_IDX_TO_ARRAY(0, 2, n), static_cast(2)); + EXPECT_EQ(MATRIX_IDX_TO_ARRAY(1, 1, n), static_cast(3)); + EXPECT_EQ(MATRIX_IDX_TO_ARRAY(1, 2, n), static_cast(4)); + EXPECT_EQ(MATRIX_IDX_TO_ARRAY(2, 2, n), static_cast(5)); +} + + +TEST(TestMatrixMacro, LargerMatrix) { + // 4x4 matrix: upper triangle has 10 elements (0..9) + uint64_t n = 4; + uint64_t expected = 0; + for (uint64_t i = 0; i < n; i++) + for (uint64_t j = i; j < n; j++) + EXPECT_EQ(MATRIX_IDX_TO_ARRAY(i, j, n), expected++); +} + + +//////////////////////////////////////////////////////////////////// +// FULL GRM FILE WRITE/READ TESTS +//////////////////////////////////////////////////////////////////// + +// Helper to build a complete Hdr + Grm for file I/O tests +static void build_test_data(grm::Hdr* hdr, grm::Grm* g) { + hdr->grm_type = grm::EHC; + + char contig[] = "chr1"; + *hdr->coords = grm::Coordinates { contig, 2 }; + hdr->coords->pos[0] = static_cast(50); + hdr->coords->pos[1] = static_cast(150); + + *hdr->samples = grm::Samples { 3 }; + hdr->samples->names[0] = "s1"; + hdr->samples->names[1] = "s2"; + hdr->samples->names[2] = "s3"; + + *g = grm::Grm { 3 }; + g->set(0, 0, 1.0f); + g->set(0, 1, 0.5f); + g->set(0, 2, 0.2f); + g->set(1, 1, 1.0f); + g->set(1, 2, 0.3f); + g->set(2, 2, 1.0f); +} + + +TEST(TestGrmFile, WriteReadRoundTrip) { + grm::Hdr hdr_w {}; + grm::Grm grm_w {}; + build_test_data(&hdr_w, &grm_w); + + io::FileIO fio { tmpfile() }; + ASSERT_NE(fio.fid, nullptr); + + ASSERT_EQ(grm::write(&fio, &hdr_w, &grm_w), grm::SUCCESS); + + rewind(fio.fid); + + grm::Hdr hdr_r {}; + grm::Grm grm_r {}; + ASSERT_EQ(grm::read(&fio, &hdr_r, &grm_r), grm::SUCCESS); + + // verify header + EXPECT_EQ(hdr_r.grm_type, grm::EHC); + EXPECT_EQ(hdr_r.coords->contig, "chr1"); + EXPECT_EQ(hdr_r.coords->len, static_cast(2)); + EXPECT_EQ(hdr_r.samples->len, static_cast(3)); + EXPECT_EQ(hdr_r.samples->names[0], "s1"); + EXPECT_EQ(hdr_r.samples->names[1], "s2"); + EXPECT_EQ(hdr_r.samples->names[2], "s3"); + + // verify grm data + EXPECT_EQ(grm_r.n_samples, static_cast(3)); + float val = 0.0f; + grm_r.get(0, 0, &val); EXPECT_FLOAT_EQ(val, 1.0f); + grm_r.get(0, 1, &val); EXPECT_FLOAT_EQ(val, 0.5f); + grm_r.get(0, 2, &val); EXPECT_FLOAT_EQ(val, 0.2f); + grm_r.get(1, 1, &val); EXPECT_FLOAT_EQ(val, 1.0f); + grm_r.get(1, 2, &val); EXPECT_FLOAT_EQ(val, 0.3f); + grm_r.get(2, 2, &val); EXPECT_FLOAT_EQ(val, 1.0f); +} + + +TEST(TestGrmFile, ReadBadMagicNumber) { + io::FileIO fio { tmpfile() }; + ASSERT_NE(fio.fid, nullptr); + + // write a wrong magic number + uint32_t bad_magic = 0x12345678; + fwrite(&bad_magic, sizeof(bad_magic), 1, fio.fid); + rewind(fio.fid); + + grm::Hdr hdr {}; + grm::Grm g {}; + EXPECT_EQ(grm::read(&fio, &hdr, &g), grm::ERROR_NOT_A_GRM_FILE); +} + + +TEST(TestGrmFile, WriteNullArgs) { + grm::Hdr hdr {}; + grm::Grm g { 2 }; + io::FileIO fio { tmpfile() }; + ASSERT_NE(fio.fid, nullptr); + + EXPECT_EQ(grm::write(nullptr, &hdr, &g), grm::ERROR_NULLPTR_ARG); + EXPECT_EQ(grm::write(&fio, static_cast(nullptr), &g), + grm::ERROR_NULLPTR_ARG); + EXPECT_EQ(grm::write(&fio, &hdr, static_cast(nullptr)), + grm::ERROR_NULLPTR_ARG); +} + + +// The read function for grm::read(fio, hdr, grm) should return +// ERROR_NULLPTR_ARG for null pointer arguments (consistent with write). +TEST(TestGrmFile, ReadNullArgs) { + grm::Hdr hdr {}; + grm::Grm g {}; + io::FileIO fio { tmpfile() }; + ASSERT_NE(fio.fid, nullptr); + + EXPECT_EQ(grm::read(nullptr, &hdr, &g), grm::ERROR_NULLPTR_ARG); + EXPECT_EQ(grm::read(&fio, static_cast(nullptr), &g), + grm::ERROR_NULLPTR_ARG); + EXPECT_EQ(grm::read(&fio, &hdr, static_cast(nullptr)), + grm::ERROR_NULLPTR_ARG); +} + + +// Reading from an empty file should fail +TEST(TestGrmFile, ReadEmptyFile) { + io::FileIO fio { tmpfile() }; + ASSERT_NE(fio.fid, nullptr); + + grm::Hdr hdr {}; + grm::Grm g {}; + EXPECT_EQ(grm::read(&fio, &hdr, &g), grm::ERROR_ON_READ); +} + + +TEST(TestGrmFile, WriteNullFid) { + grm::Hdr hdr {}; + grm::Grm g { 1 }; + io::FileIO fio { nullptr }; + + EXPECT_EQ(grm::write(&fio, &hdr, &g), grm::ERROR_NULLPTR_ARG); +} + + +TEST(TestGrmFile, ReadNullFid) { + grm::Hdr hdr {}; + grm::Grm g {}; + io::FileIO fio { nullptr }; + + EXPECT_EQ(grm::read(&fio, &hdr, &g), grm::ERROR_NULLPTR_ARG); +} + + +// Verify the magic number (FILE_TYPE_SPEC) is the first 4 bytes written +TEST(TestGrmFile, MagicNumberWrittenFirst) { + grm::Hdr hdr {}; + grm::Grm g {}; + build_test_data(&hdr, &g); + + io::FileIO fio { tmpfile() }; + ASSERT_NE(fio.fid, nullptr); + ASSERT_EQ(grm::write(&fio, &hdr, &g), grm::SUCCESS); + + rewind(fio.fid); + + uint32_t magic = 0; + size_t nread = fread(&magic, sizeof(magic), 1, fio.fid); + ASSERT_EQ(nread, static_cast(1)); + EXPECT_EQ(magic, grm::FILE_TYPE_SPEC); +} + + +// Single sample GRM file round-trip +TEST(TestGrmFile, WriteReadSingleSample) { + grm::Hdr hdr {}; + hdr.grm_type = grm::DS; + + char contig[] = "chr1"; + *hdr.coords = grm::Coordinates { contig, 1 }; + hdr.coords->pos[0] = 500; + + *hdr.samples = grm::Samples { 1 }; + hdr.samples->names[0] = "lone_sample"; + + grm::Grm g { 1 }; + g.set(0, 0, 0.75f); + + io::FileIO fio { tmpfile() }; + ASSERT_NE(fio.fid, nullptr); + ASSERT_EQ(grm::write(&fio, &hdr, &g), grm::SUCCESS); + + rewind(fio.fid); + + grm::Hdr hdr_r {}; + grm::Grm g_r {}; + ASSERT_EQ(grm::read(&fio, &hdr_r, &g_r), grm::SUCCESS); + + EXPECT_EQ(hdr_r.grm_type, grm::DS); + EXPECT_EQ(hdr_r.samples->len, static_cast(1)); + EXPECT_EQ(hdr_r.samples->names[0], "lone_sample"); + + EXPECT_EQ(g_r.n_samples, static_cast(1)); + float val = 0.0f; + g_r.get(0, 0, &val); + EXPECT_FLOAT_EQ(val, 0.75f); +} + + +// Read from a file that has the correct magic number but is +// truncated after it (no header data follows) +TEST(TestGrmFile, ReadTruncatedAfterMagic) { + io::FileIO fio { tmpfile() }; + ASSERT_NE(fio.fid, nullptr); + + fwrite(&grm::FILE_TYPE_SPEC, sizeof(grm::FILE_TYPE_SPEC), 1, fio.fid); + rewind(fio.fid); + + grm::Hdr hdr {}; + grm::Grm g {}; + EXPECT_EQ(grm::read(&fio, &hdr, &g), grm::ERROR_ON_READ); +} + + +// Write then read two independent GRM files sequentially to the same +// file, verifying both are recovered correctly +TEST(TestGrmFile, WriteReadTwoSequential) { + grm::Hdr hdr1 {}; + hdr1.grm_type = grm::EHC; + char c1[] = "chr1"; + *hdr1.coords = grm::Coordinates { c1, 1 }; + hdr1.coords->pos[0] = 10; + *hdr1.samples = grm::Samples { 2 }; + hdr1.samples->names[0] = "a"; + hdr1.samples->names[1] = "b"; + + grm::Grm g1 { 2 }; + g1.set(0, 0, 1.0f); + g1.set(0, 1, 0.5f); + g1.set(1, 1, 1.0f); + + grm::Hdr hdr2 {}; + hdr2.grm_type = grm::EAC; + char c2[] = "chr2"; + *hdr2.coords = grm::Coordinates { c2, 1 }; + hdr2.coords->pos[0] = 20; + *hdr2.samples = grm::Samples { 2 }; + hdr2.samples->names[0] = "x"; + hdr2.samples->names[1] = "y"; + + grm::Grm g2 { 2 }; + g2.set(0, 0, 2.0f); + g2.set(0, 1, 0.8f); + g2.set(1, 1, 2.0f); + + io::FileIO fio { tmpfile() }; + ASSERT_NE(fio.fid, nullptr); + + ASSERT_EQ(grm::write(&fio, &hdr1, &g1), grm::SUCCESS); + ASSERT_EQ(grm::write(&fio, &hdr2, &g2), grm::SUCCESS); + + rewind(fio.fid); + + // read first + grm::Hdr r1 {}; + grm::Grm rg1 {}; + ASSERT_EQ(grm::read(&fio, &r1, &rg1), grm::SUCCESS); + EXPECT_EQ(r1.grm_type, grm::EHC); + EXPECT_EQ(r1.samples->names[0], "a"); + float val = 0.0f; + rg1.get(0, 1, &val); + EXPECT_FLOAT_EQ(val, 0.5f); + + // read second + grm::Hdr r2 {}; + grm::Grm rg2 {}; + ASSERT_EQ(grm::read(&fio, &r2, &rg2), grm::SUCCESS); + EXPECT_EQ(r2.grm_type, grm::EAC); + EXPECT_EQ(r2.samples->names[0], "x"); + rg2.get(0, 1, &val); + EXPECT_FLOAT_EQ(val, 0.8f); +} diff --git a/tests/test_haplotype_data_record.cpp b/tests/test_haplotype_data_record.cpp deleted file mode 100644 index 4f0ce29..0000000 --- a/tests/test_haplotype_data_record.cpp +++ /dev/null @@ -1,117 +0,0 @@ - - -#include "../include/HaplotypeVcfParser.h" -#include - - - - -TEST(TestConstructorAssignment, Constructor) { - - size_t num_vcf_columns { 13 }; - size_t num_samps { 4 }; - size_t num_founders { 3 }; - static constexpr int num_records { 3 }; - - - char vcf_record[] { "chr12 1 . A T Q1 F1 INFO1 GT:AB:HD 0/0:2:1,0,1 0/0:1:0,2,0 1/0:0:0,0,2 1/1:1:2,0,0\n" }; - - - // Test first record - HaplotypeDataRecord hap_record { num_samps, num_founders }; - - std::array dims { hap_record.dims() }; - - EXPECT_EQ(dims[0], num_samps); - EXPECT_EQ(dims[1], num_founders); - - hap_record.parse_vcf_line(vcf_record); - - EXPECT_EQ(hap_record(0,0), 1); - EXPECT_EQ(hap_record(1,0), 0); - EXPECT_EQ(hap_record(2,0), 0); - EXPECT_EQ(hap_record(3,0), 2); - - EXPECT_EQ(hap_record(0,1), 0); - EXPECT_EQ(hap_record(1,1), 2); - EXPECT_EQ(hap_record(2,1), 0); - EXPECT_EQ(hap_record(3,1), 0); - - EXPECT_EQ(hap_record(0,2), 1); - EXPECT_EQ(hap_record(1,2), 0); - EXPECT_EQ(hap_record(2,2), 2); - EXPECT_EQ(hap_record(3,2), 0); - - // Test second record - - char vcf_record2[] { "chr12 2 . G T Q2 F2 INFO2 GT:AB:HD 0/0:2:1,0,1 0/1:1:0,2,0 1/0:0:0,0,2 0/1:1:0,0,2\n" }; - HaplotypeDataRecord hap_record2 { num_samps, num_founders }; - - hap_record2.parse_vcf_line(vcf_record2); - - EXPECT_EQ(hap_record2(0,0), 1); - EXPECT_EQ(hap_record2(1,0), 0); - EXPECT_EQ(hap_record2(2,0), 0); - EXPECT_EQ(hap_record2(3,0), 0); - - EXPECT_EQ(hap_record2(0,1), 0); - EXPECT_EQ(hap_record2(1,1), 2); - EXPECT_EQ(hap_record2(2,1), 0); - EXPECT_EQ(hap_record2(3,1), 0); - - EXPECT_EQ(hap_record2(0,2), 1); - EXPECT_EQ(hap_record2(1,2), 0); - EXPECT_EQ(hap_record2(2,2), 2); - EXPECT_EQ(hap_record2(3,2), 2); - -} - - - -// TEST(TestConstructorAssignment, CopyConstructor) { -// -// size_t num_vcf_columns { 13 }; -// size_t num_founders { 3 }; -// size_t num_samps { 4 }; -// -// -// char vcf_record[] { "chr12 1 . A T Q1 F1 INFO1 GT:AB:HD 0/0:2:1,0,1 0/0:1:0,2,0 1/0:0:0,0,2 1/1:1:2,0,0\n" }; -// -// HaplotypeDataRecord origin { num_founders, num_samps }; -// HaplotypeDataRecord replicate { origin }; -// -// // test that the objects are indeed different -// EXPECT_NE(&origin, &replicate); -// -// std::array o_dims { origin.dims() }; -// std::array rep_dims { replicate.dims() }; -// -// EXPECT_EQ(o_dims[0], rep_dims[0]); -// EXPECT_EQ(o_dims[1], rep_dims[1]); -// -// -// // make sure that copy constructor correct -// EXPECT_EQ(origin(0,0), replicate(0,0)); -// EXPECT_EQ(origin(0,1), replicate(0,1)); -// EXPECT_EQ(origin(0,2), replicate(0,2)); -// EXPECT_EQ(origin(0,3), replicate(0,3)); -// -// -// EXPECT_EQ(origin(1,0), replicate(1,0)); -// EXPECT_EQ(origin(1,1), replicate(1,1)); -// EXPECT_EQ(origin(1,2), replicate(1,2)); -// EXPECT_EQ(origin(1,3), replicate(1,3)); -// -// -// EXPECT_EQ(origin(2,0), replicate(2,0)); -// EXPECT_EQ(origin(2,1), replicate(2,1)); -// EXPECT_EQ(origin(2,2), replicate(2,2)); -// EXPECT_EQ(origin(2,3), replicate(2,3)); -// -// } - -// std::array meta { "##version=4.2\n", -// "##ID=\n", -// "##FORMAT=\n" }; -// -// std::string h { "#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT S1 S2 S3 S4\n" }; diff --git a/tests/test_haplotype_vcf_parser.cpp b/tests/test_haplotype_vcf_parser.cpp deleted file mode 100644 index 0f8e5ce..0000000 --- a/tests/test_haplotype_vcf_parser.cpp +++ /dev/null @@ -1,116 +0,0 @@ - -#include "../include/HaplotypeVcfParser.h" -#include "../include/utils.h" -#include - - - -char VCF_NAME[] { "../tests/test.vcf" }; - - -TEST(TestHaplotypeVCFParser, Constructor) { - - HaplotypeVcfParser vcf { VCF_NAME }; - - EXPECT_EQ(vcf.n_samples(), 11); - EXPECT_EQ(vcf.k_founders(), 8); -} - - -TEST(TestHaplotypeVCFParser, LoadRecord) { - - HaplotypeVcfParser vcf { VCF_NAME }; - - HaplotypeDataRecord record { vcf.n_samples(), vcf.k_founders() }; - - bool record_loaded { false }; - record_loaded = vcf.load_record(record); - - EXPECT_TRUE(record_loaded); - - EXPECT_EQ(record.chrom(), "chr12"); - EXPECT_EQ(record.pos(), 788); - EXPECT_EQ(record.id(), "."); - EXPECT_EQ(record.ref(), 'A'); - EXPECT_EQ(record.alt(), 'G'); - EXPECT_EQ(record.qual(), "."); - EXPECT_EQ(record.filter(), "PASS"); - EXPECT_EQ(record.info(), "EAF=0.00228;INFO_SCORE=1;HWE=1;ERC=0.01949;EAC=7.94153;PAF=0.00245;REF_PANEL=0"); - EXPECT_EQ(record.format(), "GT:GP:DS:HD"); - - EXPECT_EQ(record(0,0), 1.004); - EXPECT_EQ(record(0,1), 0); - EXPECT_EQ(record(0,2),0.002); - EXPECT_EQ(record(0,3),0.001); - EXPECT_EQ(record(0,4),0); - EXPECT_EQ(record(0,5),0.991); - EXPECT_EQ(record(0,6), 0.001); - EXPECT_EQ(record(0,7), 0.002); - - EXPECT_EQ(record(1,0), 0.001); - EXPECT_EQ(record(1,1), 0); - EXPECT_EQ(record(1,2),0); - EXPECT_EQ(record(1,3),0.989); - EXPECT_EQ(record(1,4),0.005); - EXPECT_EQ(record(1,5),0.005); - EXPECT_EQ(record(1,6), 1); - EXPECT_EQ(record(1,7), 0); - - EXPECT_EQ(record(2,0), 0.998); - EXPECT_EQ(record(2,1), 0); - EXPECT_EQ(record(2,2),0.001); - EXPECT_EQ(record(2,3),0); - EXPECT_EQ(record(2,4),0); - EXPECT_EQ(record(2,5),0); - EXPECT_EQ(record(2,6), 1); - EXPECT_EQ(record(2,7), 0); - - EXPECT_EQ(record(3,0), 0.84); - EXPECT_EQ(record(3,1), 0); - EXPECT_EQ(record(3,2),0); - EXPECT_EQ(record(3,3),0); - EXPECT_EQ(record(3,4),0); - EXPECT_EQ(record(3,5),0); - EXPECT_EQ(record(3,6), 1.159); - EXPECT_EQ(record(3,7), 0); - - EXPECT_EQ(record(10,0), 0.592); - EXPECT_EQ(record(10,1), 0); - EXPECT_EQ(record(10,2),1); - EXPECT_EQ(record(10,3),0); - EXPECT_EQ(record(10,4),0); - EXPECT_EQ(record(10,5),0); - EXPECT_EQ(record(10,6),0); - EXPECT_EQ(record(10,7),0.407); - - // load second record - record_loaded = vcf.load_record(record); - EXPECT_EQ(record.chrom(), "chr12"); - EXPECT_EQ(record.pos(), 1321); - EXPECT_EQ(record.id(), "."); - EXPECT_EQ(record.ref(), 'A'); - EXPECT_EQ(record.alt(), 'C'); - EXPECT_EQ(record.qual(), "."); - EXPECT_EQ(record.filter(), "PASS"); - EXPECT_EQ(record.info(), "EAF=0.01487;INFO_SCORE=0.17212;HWE=1;ERC=1.33325;EAC=116.998;PAF=0.01127;REF_PANEL=0"); - EXPECT_EQ(record.format(), "GT:GP:DS:HD"); - - EXPECT_EQ(record(0,0),1.004); - EXPECT_EQ(record(0,1), 0); - EXPECT_EQ(record(0,2),0.002); - EXPECT_EQ(record(0,3),0.001); - EXPECT_EQ(record(0,4),0); - EXPECT_EQ(record(0,5),0.991); - EXPECT_EQ(record(0,6), 0.001); - EXPECT_EQ(record(0,7), 0.002); - - EXPECT_EQ(record(10,0), 0.592); - EXPECT_EQ(record(10,1), 0); - EXPECT_EQ(record(10,2),1); - EXPECT_EQ(record(10,3),0); - EXPECT_EQ(record(10,4),0); - EXPECT_EQ(record(10,5),0); - EXPECT_EQ(record(10,6),0); - EXPECT_EQ(record(10,7),0.407); - -} diff --git a/tests/test_matrix.cpp b/tests/test_matrix.cpp deleted file mode 100644 index bbc5e9f..0000000 --- a/tests/test_matrix.cpp +++ /dev/null @@ -1,87 +0,0 @@ - -#include "../include/Matrix.h" -#include - - -TEST(TestMatrix, initialize) { - size_t n_row { 3 }; - size_t m_col { 2 }; - - Matrix a { n_row, m_col }; - - std::array dims { a.dims() }; - EXPECT_EQ(dims[0], n_row); - EXPECT_EQ(dims[1], m_col); - - for (size_t i = 0; i < n_row; i++) - for (size_t j = 0; j < m_col; j++) - EXPECT_EQ(a(i, j), 0); - - n_row = 0; - - EXPECT_THROW({ - size_t n_row = 0; - size_t m_col = 2; - Matrix b(n_row, m_col); - }, - std::runtime_error); - - EXPECT_ANY_THROW({ - size_t n_row = -1; - size_t m_col = 2; - Matrix b(n_row, m_col); - }); - - EXPECT_THROW({Matrix b(1, 0);}, std::runtime_error); - EXPECT_ANY_THROW({Matrix b(1, -1);}); -} - - - -TEST(TestMatrix, val) { - size_t n_row { 3 }; - size_t m_col { 5 }; - - Matrix a { n_row, m_col }; - - std::array dims { a.dims() }; - - double x = { 1 }; - for (size_t i = 0; i < dims[0]; i++) - for (size_t j = 0; j < dims[1]; j++) - a(i, j) = x++; - - x = 1; - for (size_t i = 0; i < dims[0]; i++) - for (size_t j = 0; j < dims[1]; j++) - EXPECT_FLOAT_EQ(a(i, j), x++); -} - - -TEST(TestMatrix, out_of_bounds) { - size_t n_row { 3 }; - size_t m_col { 5 }; - - Matrix a { n_row, m_col }; - - EXPECT_THROW({ double tmp = a(4, 3); }, std::runtime_error); - EXPECT_THROW({ double tmp = a(3, 5); }, std::runtime_error); - EXPECT_THROW({ double tmp = a(3, 2); }, std::runtime_error); - EXPECT_THROW({ double tmp = a(2, 5); }, std::runtime_error); - EXPECT_THROW({ double tmp = a(-2, 4); }, std::runtime_error); - -} - - -TEST(TestMatrix, dim_and_size) { - size_t n_row { 3 }; - size_t m_col { 5 }; - - Matrix a { n_row, m_col }; - - std::array dims { a.dims() }; - EXPECT_EQ(dims[0], n_row); - EXPECT_EQ(dims[1], m_col); - - EXPECT_EQ(a.size(), n_row * m_col); -} diff --git a/tests/test_utils.cpp b/tests/test_utils.cpp deleted file mode 100644 index 50043fe..0000000 --- a/tests/test_utils.cpp +++ /dev/null @@ -1,97 +0,0 @@ - -#include "../include/utils.h" -#include - - - -TEST(TestCharBuffer, Constructor) { - size_t buff_size { 100 }; - CharBuffer buff { buff_size }; - - EXPECT_EQ(buff.buffer_size(), buff_size); -} - - -TEST(TestCharBuffer, Append) { - size_t buff_size { 5 }; - CharBuffer buff { buff_size }; - - buff.append('t'); - EXPECT_EQ(buff(0), 't'); - - buff.append('h'); - buff.append('e'); - EXPECT_EQ(buff(0), 't'); - EXPECT_EQ(buff(1), 'h'); - EXPECT_EQ(buff(2), 'e'); - - EXPECT_THROW(buff(6), std::out_of_range); - EXPECT_THROW(buff(5), std::out_of_range); - EXPECT_THROW(buff(-1), std::out_of_range); - - buff.append('r'); - buff.append('e'); - EXPECT_THROW(buff.append('e'), std::out_of_range); -} - - -TEST(TestCharBuffer, Reset) { - size_t buff_size { 5 }; - CharBuffer buff { buff_size }; - - buff.append('t'); - EXPECT_EQ(buff(0), 't'); - - buff.reset(); - EXPECT_THROW(buff(0), std::out_of_range); -} - - -TEST(TestCharBuffer, Data) { - size_t buff_size { 5 }; - CharBuffer buff { buff_size }; - buff.append('t'); - buff.append('h'); - buff.append('e'); - - EXPECT_EQ(buff.data()[0], 't'); - - std::string s { buff.data() }; - EXPECT_EQ(s, "the"); - - buff.append('e'); - s = buff.data(); - EXPECT_EQ(s, "thee"); - -} - - -TEST(TestStringRecord, ConstructorDelim) { - CharBuffer buf { 10 }; - char s[] { "the\tcat\n" }; - StringRecord record { '\t' }; - - record.update_str(s); - - EXPECT_TRUE(record.next_field()); - EXPECT_EQ(static_cast(record.data()), "the"); - - EXPECT_TRUE(record.next_field()); - EXPECT_EQ(static_cast(record.data()), "cat"); - - EXPECT_FALSE(record.next_field()); - - EXPECT_EQ(record.size(), 8); - - record.reset(); - - EXPECT_EQ(record.size(), 0); -} - -TEST(TestStringRecord, SingleArgConstructor) { - StringRecord record { '\t' }; - - EXPECT_EQ(record.size(), 0); - EXPECT_FALSE(record.next_field()); - -}