[med-svn] [Git][med-team/gemma][master] 5 commits: routine-update: New upstream version
Nilesh Patra
gitlab at salsa.debian.org
Tue Feb 2 20:11:32 GMT 2021
Nilesh Patra pushed to branch master at Debian Med / gemma
Commits:
75edcfd0 by Nilesh Patra at 2021-02-03T01:24:05+05:30
routine-update: New upstream version
- - - - -
21f321c9 by Nilesh Patra at 2021-02-03T01:24:12+05:30
New upstream version 0.98.4+dfsg
- - - - -
6785eb98 by Nilesh Patra at 2021-02-03T01:26:04+05:30
Update upstream source from tag 'upstream/0.98.4+dfsg'
Update to upstream version '0.98.4+dfsg'
with Debian dir f49962a5842ce7593c7eff0d268331d0f4929eb0
- - - - -
04648604 by Nilesh Patra at 2021-02-03T01:28:59+05:30
Refresh patches
- - - - -
8cf6fd2b by Nilesh Patra at 2021-02-03T01:29:26+05:30
routine-update: Ready to upload to unstable
- - - - -
22 changed files:
- + .guix-dev-static
- .travis.yml
- INSTALL.md
- Makefile
- README.md
- RELEASE-NOTES.md
- VERSION
- debian/changelog
- debian/patches/Makefile.patch
- debian/patches/Reproducible_builds.patch
- debian/patches/no-lquadmath.patch
- doc/manual.pdf
- doc/manual.tex
- + guix.scm
- src/gemma.cpp
- src/gemma.h
- src/lmm.cpp
- src/lmm.h
- src/param.cpp
- src/param.h
- src/version.h
- test/dev_test_suite.sh
Changes:
=====================================
.guix-dev-static
=====================================
@@ -0,0 +1,2 @@
+# Typical Guix container invocation
+~/opt/guix/bin/guix environment -C guix --ad-hoc gcc-toolchain:static gdb gsl openblas zlib:static bash ld-wrapper perl vim which gfortran-toolchain
=====================================
.travis.yml
=====================================
@@ -2,6 +2,8 @@ language: C++
matrix:
include:
- os: linux
+ arch:
+ - amd64
compiler: gcc
addons:
apt:
@@ -21,9 +23,21 @@ matrix:
- libopenblas-dev
- os: osx
compiler: clang
+ - os: linux
+ arch:
+ - arm64
+ compiler: gcc
+ addons:
+ apt:
+ sources:
+ - ubuntu-toolchain-r-test
+ packages:
+ - libgsl-dev
+ - libopenblas-dev
script:
- echo $MATRIX_EVAL
- eval "${MATRIX_EVAL}"
+ - if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then brew install gsl ; fi
- $CXX --version
- if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then brew install gsl ; fi
- if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then OPENBLAS="$(brew --prefix openblas)" ; fi
=====================================
INSTALL.md
=====================================
@@ -10,8 +10,8 @@ and it should give you the version.
## GEMMA dependencies
-GEMMA runs on Linux and MAC OSX and the runtime has the following
-dependencies:
+GEMMA runs on Linux, MAC OSX and Windows (with Docker). The runtime
+has the following dependencies:
* C++ tool chain >= 5.5.0 (see Travis CI and we test with file .guix-dev-gcc-older)
* GNU Science library (GSL) 2.x (GEMMA dropped support for GSL 1.x)
@@ -52,41 +52,37 @@ using the following command
guix package -i gemma
-To build GEMMA from source you can opt to install the build tools with GNU Guix
+A more recent version may be found in the guix-bioinformatics channel
+which is maintained by the authors. See the
+[README](http://git.genenetwork.org/guix-bioinformatics/guix-bioinformatics), e.g.
- guix package -i make gcc linux-libre-headers gsl openblas lapack glibc ld-wrapper
+ env GUIX_PACKAGE_PATH=./guix-bioinformatics guix package -A gemma
-The current build container is in [guix-dev](../guix-dev)
+To build GEMMA from source you can opt to install the build tools with
+GNU Guix, the current build container is in [guix-dev](./.guix-dev)
- guix environment -C guix --ad-hoc gcc-toolchain gdb gsl openblas zlib bash ld-wrapper perl vim which
+ source .guix-dev
+ make
-To build with an older gcc, for example:
+Guix allows for easy versioning. To build with an older gcc, for
+example:
guix environment -C guix --ad-hoc gcc-toolchain at 9.3.0 gdb gsl openblas zlib bash ld-wrapper perl vim which
-### Install from source
+### Install with Docker
+
+Recent version of GEMMA come with a 64-bit Docker image that should run
+on Linux, Windows and MacOS.
-Note: Eigen is no longer required!
+### Install from source
-Install listed dependencies and run
+Install listed dependencies (you may want to take hints from
+the Travis-CI [tests](./.travis.yml)) and run
make -j 4
(the -j switch builds on 4 cores).
-if you get an Eigen error you may need to override the include
-path. E.g. to build GEMMA on GNU Guix with shared libs the following
-may work
-
- make EIGEN_INCLUDE_PATH=~/.guix-profile/include/eigen3 WITH_OPENBLAS=1
-
-another example overriding optimization and LIB flags (so as to link
-against gslv1) would be
-
- make EIGEN_INCLUDE_PATH=~/.guix-profile/include/eigen3 WITH_OPENBLAS=1 GCC_FLAGS="-Wall" LIBS="$HOME/opt/gsl1/lib/libgsl.a $HOME/opt/gsl1/lib/libgslcblas.a -L$HOME/.guix-profile/lib -pthread -llapack -lblas -lz"
-
-to run GEMMA tests
-
time make check
You can run gemma in the debugger with, for example
@@ -99,27 +95,23 @@ You can run gemma in the debugger with, for example
Note that if you get <optimized out> warnings on inspecting variables you
should compile with GCC_FLAGS="" to disable optimizations (-O3). E.g.
- make EIGEN_INCLUDE_PATH=~/.guix-profile/include/eigen3 WITH_OPENBLAS=1 GCC_FLAGS=
+ make WITH_OPENBLAS=1 GCC_FLAGS=
Other options, such as compiling with warnings, are listed in the
Makefile.
### GNU Guix commands used
-With git bisect build the older versions of gemma with openblas
+Some development examples. With git bisect build the older versions
+of gemma with openblas
- ~/.config/guix/current/bin/guix environment -C guix --ad-hoc gcc gdb gfortran:lib gsl eigen lapack openblas zlib bash ld-wrapper perl ldc
- make clean ; make EIGEN_INCLUDE_PATH=$GUIX_ENVIRONMENT/include/eigen3 WITH_OPENBLAS=1 FORCE_DYNAMIC=1 -j 8
+ ~/.config/guix/current/bin/guix environment -C guix --ad-hoc gcc gdb gfortran:lib gsl lapack openblas zlib bash ld-wrapper perl ldc
+ make clean ; make WITH_OPENBLAS=1 FORCE_DYNAMIC=1 -j 8
or with atlas
- ~/.config/guix/current/bin/guix environment -C guix --ad-hoc gcc gdb gfortran:lib gsl eigen lapack atlas zlib bash ld-wrapper perl ldc
- make clean ; make EIGEN_INCLUDE_PATH=$GUIX_ENVIRONMENT/include/eigen3/Eigen/ WITH_OPENBLAS= FORCE_DYNAMIC=1 -j 25
-
-You may need to symlink Eigen in some older versions
-
- ln -s $GUIX_ENVIRONMENT/include/eigen3/Eigen src/Eigen
-
+ ~/.config/guix/current/bin/guix environment -C guix --ad-hoc gcc gdb gfortran:lib gsl lapack atlas zlib bash ld-wrapper perl ldc
+ make clean ; make WITH_OPENBLAS= FORCE_DYNAMIC=1 -j 25
## Run tests
@@ -130,3 +122,31 @@ GEMMA includes the shunit2 test framework (version 2.0).
or
./run_tests.sh
+
+## Releases
+
+### Docker release
+
+To distribute GEMMA I made static versions of the binary. A container
+can be made instead with, for example
+
+```sh
+env GUIX_PACKAGE_PATH=~/guix-bioinformatics ~/.config/guix/current/bin/guix \
+ pack -f docker gemma-gn2 -S /bin=bin
+```
+
+which created a container in of size 51MB. Tiny! For more information
+see
+[GUIX-NOTES](http://git.genenetwork.org/guix-bioinformatics/guix-notes/CONTAINERS.org).
+
+
+### Static release
+
+To create a static release, locate the gfortran lib and use
+
+ source .guix-dev-static
+ make WITH_GFORTRAN=1 EXTRA_FLAGS=-L/gnu/store/741057r2x06zwg6zcmqmdyv51spm6n9i-gfortran-7.5.0-lib/lib static
+
+otherwise OpenBlas will complain with
+
+ undefined reference to `_gfortran_concat_string'
=====================================
Makefile
=====================================
@@ -23,6 +23,10 @@
#
# make static
#
+# Static build may require some extra flags
+#
+# make WITH_GFORTRAN=1 EXTRA_FLAGS=-L/gnu/store/741057r2x06zwg6zcmqmdyv51spm6n9i-gfortran-7.5.0-lib/lib static
+#
# Run tests with
#
# make check
@@ -72,6 +76,7 @@ WITH_GFORTRAN = # Add -lgfortran (if OpenBlas does not
OPENBLAS_LEGACY = # Using older OpenBlas
FORCE_STATIC = # Static linking of libraries
GCC_FLAGS = -DHAVE_INLINE -pthread -Wall -std=gnu++11 # extra flags -Wl,--allow-multiple-definition
+EXTRA_FLAGS =
GSL_INCLUDE_PATH =
ifeq ($(SYS), WIN)
@@ -172,6 +177,9 @@ endif
ifdef WITH_GFORTRAN
LIBS += -lgfortran -lquadmath
endif
+ifdef EXTRA_FLAGS
+ LIBS += $(EXTRA_FLAGS)
+endif
.PHONY: all test
=====================================
README.md
=====================================
@@ -3,28 +3,32 @@ Nat. Genet., 2016)](cfw.gif)
# GEMMA: Genome-wide Efficient Mixed Model Association
-[](https://travis-ci.org/genetics-statistics/GEMMA) [](https://anaconda.org/bioconda/gemma) [](https://anaconda.org/bioconda/gemma) [](https://github.com/brewsci/homebrew-bio) [](https://www.gnu.org/software/guix/packages/G/) [](https://packages.debian.org/testing/gemma)
+[](https://travis-ci.com/genetics-statistics/GEMMA) [](https://anaconda.org/bioconda/gemma) [](https://anaconda.org/bioconda/gemma) [](https://github.com/brewsci/homebrew-bio) [](https://www.gnu.org/software/guix/packages/G/) [](https://packages.debian.org/search?keywords=gemma&searchon=names&suite=all§ion=all)
GEMMA is a software toolkit for fast application of linear mixed
models (LMMs) and related models to genome-wide association studies
(GWAS) and other large-scale data sets.
-Check out [NEWS.md](NEWS.md) to see what's new in each GEMMA release.
+Check out [RELEASE-NOTES.md](RELEASE-NOTES.md) to see what's new in
+each GEMMA release.
-Please post feature requests or suspected bugs to
+Please post suspected bugs to
[Github issues](https://github.com/genetics-statistics/GEMMA/issues). For
questions or other discussion, please post to the
[GEMMA Google Group](https://groups.google.com/group/gemma-discussion). We
also encourage contributions, for example, by forking the repository,
making your changes to the code, and issuing a pull request.
-Currently, GEMMA is supported for 64-bit Mac OS X and Linux
-platforms. *Windows is not currently supported.* though you can
-run GEMMA in a Linux VM or [container](https://docs.docker.com/docker-for-windows/). If you are interested
-in helping to make GEMMA available on Windows platforms (e.g., by
-providing installation instructions for Windows, or by contributing
-Windows binaries) please post a note in the
-[Github issues](https://github.com/genetics-statistics/GEMMA/issues).
+Currently, GEMMA provides a runnable Docker container for 64-bit
+MacOS, Windows and Linux platforms. GEMMA can be installed with
+Debian, Conda, Homebrew and GNU Guix. With Guix you find the latest
+version
+[here](http://git.genenetwork.org/guix-bioinformatics/guix-bioinformatics)
+as it is the version we use every day on http://genenetwork.org. For
+installation instructions see also [INSTALL.md](INSTALL.md). We use
+continous integration builds on Travis-CI for Linux (amd64 & arm64)
+and MacOS (amd64). GEMMA builds on multiple architectures, see the
+[Debian build farm](https://buildd.debian.org/status/package.php?p=gemma).
*(The above image depicts physiological and behavioral trait
loci identified in CFW mice using GEMMA, from [Parker et al, Nature
@@ -50,19 +54,19 @@ Genetics, 2016](https://doi.org/10.1038/ng.3609).)
1. Fast assocation tests implemented using the univariate linear mixed
model (LMM). In GWAS, this can correct for population structure and
-sample nonexchangeability. It also provides estimates of the
+sample non-exchangeability. It also provides estimates of the
proportion of variance in phenotypes explained by available genotypes
(PVE), often called "chip heritability" or "SNP heritability".
2. Fast association tests for multiple phenotypes implemented using a
multivariate linear mixed model (mvLMM). In GWAS, this can correct for
-populations tructure and sample nonexchangeability jointly in multiple
-complex phenotypes.
+population structure and sample (non)exchangeability - jointly in
+multiple complex phenotypes.
3. Bayesian sparse linear mixed model (BSLMM) for estimating PVE,
phenotype prediction, and multi-marker modeling in GWAS.
-4. Estimation of variance components ("chip heritability") partitioned
+4. Estimation of variance components ("chip/SNP heritability") partitioned
by different SNP functional categories from raw (individual-level)
data or summary data. For raw data, HE regression or the REML AI
algorithm can be used to estimate variance components when
@@ -73,7 +77,8 @@ MQS algorithm to estimate variance components.
To install GEMMA you can
-1. Download the precompiled binaries (64-bit Linux and Mac only)
+1. Download the precompiled or Docker binaries
+ from [releases](https://github.com/genetics-statistics/GEMMA/releases).
2. Use existing package managers, see [INSTALL.md](INSTALL.md).
@@ -88,20 +93,16 @@ numerical libraries.
1. Fetch the [latest stable release][latest_release] and download the
file appropriate for your platform.
-2. For .tar.bz2 files unpack the tar ball
+2. For Docker images, install Docker, load the image into Docker and
+ run with something like
- tar xvjf gemma-$version-installer.tar.bz2
-
- run the installer
-
- ./install.sh ~/gemma
-
- and run gemma
-
- ~/gemma/bin/gemma
+ docker run -w /run -v ${PWD}:/run ed5bf7499691 gemma -gk -bfile example/mouse_hs1940
3. For .gz files run `gunzip gemma.linux.gz` or `gunzip
-gemma.linux.gz` to unpack the file.
+gemma.linux.gz` to unpack the file. And make sure it is executable with
+
+ chmod u+x gemma-linux
+ ./gemma-linux
## Run GEMMA
@@ -131,24 +132,25 @@ Above example files can be downloaded from
GEMMA has a wide range of debugging options which can be viewed with
```
-gemma -h 14
-
DEBUG OPTIONS
+
-check enable checks (slower)
-no-fpe-check disable hardware floating point checking
-strict strict mode will stop when there is a problem
-silence silent terminal display
-debug debug output
-debug-data debug data output
+ -nind [num] read up to num individuals
+ -issue [num] enable tests relevant to issue tracker
-legacy run gemma in legacy mode
```
-typically when running gemma you should use -debug which includes relevant
-checks.
+typically when running gemma you should use -debug which includes
+relevant checks. When compiled for debugging the debug version of
+GEMMA gives more information.
-For performances you may want to use the -no-check option
-instead. Also check the build optimization notes in
-[INSTALL.md](INSTALL.md).
+For performance you may want to use the -no-check option. Also check
+the build optimization notes in [INSTALL.md](INSTALL.md).
## Help
@@ -191,7 +193,7 @@ studies.](https://doi.org/10.1101/042846) *Annals of Applied Statistics*, in pre
## License
-Copyright (C) 2012–2018, Xiang Zhou and team.
+Copyright (C) 2012–2020, Xiang Zhou and team.
The *GEMMA* source code repository is free software: you can
redistribute it under the terms of the
=====================================
RELEASE-NOTES.md
=====================================
@@ -4,6 +4,15 @@ see
and
[commits](https://github.com/genetics-statistics/GEMMA/commits/master).
+## ChangeLog v0.98.4 (2021/01/29)
+
+* Fix error on free with randomizer, see #239
+* Moved travis-ci.org to travis-ci.com
+* GEMMA builds on ARM and other architectures, see #189 and https://buildd.debian.org/status/package.php?p=gemma (thanks @tillea)
+* Fixed static build with 00480e8549987b6cae7100b28bcead2a2d501177 - requires gfortran path for OpenBLAS
+* Updated README's and Manual
+* Added `-lmm 9` switch which shows ~beta/se~ with ~lmle~ and ~plrt~, see #237
+
## ChangeLog v0.98.3 (2020/11/28)
Maintenance release
@@ -21,6 +30,15 @@ Maintenance release
[docs](https://www.gnu.org/software/gsl/doc/html/rng.html).
* The tests now use a fixed seed for the randomizer
+A docker binary that runs on Linux, MaxOS and Windows can be downloaded from
+
+http://ipfs.genenetwork.org/ipfs/Qmaq1q73ox53ykKdRF6tYDXL9bEKJQfnGCqBxFdo1fcYPb/gemma-0.98.3-AMD64-Guix-docker-release.tgz
+
+After loading the image into Docker, run with something like
+
+ docker run -w /run -v ${PWD}:/run ed5bf7499691 gemma -gk -bfile example/mouse_hs1940
+
+
## ChangeLog v0.98.2 (2019/05/28)
GCC 10.1 fix release
=====================================
VERSION
=====================================
@@ -1 +1 @@
-0.98.3
+0.98.4
=====================================
debian/changelog
=====================================
@@ -1,3 +1,10 @@
+gemma (0.98.4+dfsg-1) unstable; urgency=medium
+
+ * Team upload.
+ * New upstream version
+
+ -- Nilesh Patra <npatra974 at gmail.com> Wed, 03 Feb 2021 01:29:26 +0530
+
gemma (0.98.3+dfsg-1) unstable; urgency=medium
* Team upload.
=====================================
debian/patches/Makefile.patch
=====================================
@@ -4,7 +4,7 @@ Last-Update: 2018-04-26
--- a/Makefile
+++ b/Makefile
-@@ -219,7 +219,7 @@ $(OBJS): $(HDR)
+@@ -225,7 +225,7 @@
.SUFFIXES : .cpp .c .o $(SUFFIXES)
@@ -13,7 +13,7 @@ Last-Update: 2018-04-26
$(CPP) $(CPPFLAGS) $(TEST_SRC_DIR)/unittests-main.o $(TEST_SRC_DIR)/unittests-math.o $(filter-out src/main.o, $(OBJS)) $(LIBS) -o ./bin/unittests-gemma
unittests: all ./bin/unittests-gemma
-@@ -228,12 +228,12 @@ unittests: all ./bin/unittests-gemma
+@@ -234,12 +234,12 @@
fast-check: all unittests
rm -vf test/output/*
cd test && ./dev_test_suite.sh | tee ../dev_test.log
=====================================
debian/patches/Reproducible_builds.patch
=====================================
@@ -3,11 +3,9 @@ Description: Patch the makefile in order to enable reproducible builds (stable o
Last-Update: 2018-06-14
Forwarded: TODO
-Index: gemma-0.98+dfsg/Makefile
-===================================================================
---- gemma-0.98+dfsg.orig/Makefile
-+++ gemma-0.98+dfsg/Makefile
-@@ -215,8 +215,8 @@ ifdef WITH_LAPACK
+--- a/Makefile
++++ b/Makefile
+@@ -201,8 +201,8 @@
endif
endif
=====================================
debian/patches/no-lquadmath.patch
=====================================
@@ -6,12 +6,12 @@ Bug-Debian: https://bugs.debian.org/902153
--- a/Makefile
+++ b/Makefile
-@@ -172,7 +172,7 @@ ifdef WITH_GSLCBLAS
+@@ -175,7 +175,7 @@
LIBS += -lgslcblas
endif
ifdef WITH_GFORTRAN
- LIBS += -lgfortran -lquadmath
+ LIBS += -lgfortran
endif
-
- .PHONY: all test
+ ifdef EXTRA_FLAGS
+ LIBS += $(EXTRA_FLAGS)
=====================================
doc/manual.pdf
=====================================
Binary files a/doc/manual.pdf and b/doc/manual.pdf differ
=====================================
doc/manual.tex
=====================================
@@ -75,7 +75,7 @@ association studies (GWAS). It fits a univariate linear mixed model
(LMM) for marker association tests with a single phenotype to account
for population stratification and sample structure, and for estimating
the proportion of variance in phenotypes explained (PVE) by typed
-genotypes (i.e. "chip heritability") \cite{Zhou:2012}. It fits a
+genotypes (i.e. ``chip heritability'' or ``SNP heritability'') \cite{Zhou:2012}. It fits a
multivariate linear mixed model (mvLMM) for testing marker
associations with multiple phenotypes simultaneously while controlling
for population stratification, and for estimating genetic correlations
@@ -124,7 +124,8 @@ labels) for $n$ individuals; $\bW=(\bw_1, \cdots, \bw_c)$ is an
$n\times c$ matrix of covariates (fixed effects) including a column of
1s; $\boldsymbol \alpha$ is a $c$-vector of the corresponding
coefficients including the intercept; $\bx$ is an $n$-vector of marker
-genotypes; $\beta$ is the effect size of the marker; $\bu$ is an
+genotypes; $\beta$ is the effect size of the marker and is an estimate
+ of the marker/SNP additive effect; $\bu$ is an
$n$-vector of random effects; $\bepsilon$ is an $n$-vector of errors;
$\tau^{-1}$ is the variance of the residual errors; $\lambda$ is the
ratio between the two variance components; $\bK$ is a known $n\times
@@ -139,8 +140,8 @@ score). GEMMA obtains either the maximum likelihood estimate (MLE) or
the restricted maximum likelihood estimate (REML) of $\lambda$ and
$\beta$, and outputs the corresponding $p$ value.
-In addition, GEMMA estimates the PVE by typed genotypes or ``chip
-heritability".
+In addition, GEMMA estimates the PVE by typed genotypes or ``chip or
+SNP heritability''.
\subsubsection{Multivariate Linear Mixed Model}
GEMMA can fit a multivariate linear mixed model in the following form:
@@ -307,19 +308,12 @@ platform.
The binary executable of GEMMA works well for a reasonably large
number of individuals (say, for example, the ``-eigen " option works
-for at least 45,000 individuals). Due to the outdated computation
-environment the software was compiled on, however, for larger sample
-size and for improved computation efficiency, it is recommended to
-compile GEMMA on user's own modern computer system.
+for at least 45,000 individuals).
If you want to compile GEMMA by yourself, you will need to download
the source code, and you will need a standard C/C++ compiler such as
-GNU gcc, as well as the GSL and LAPACK libraries. You will need to
-change the library paths in the Makefile accordingly. A sample
-Makefile is provided along with the source code. For details on
-installing GSL library, please refer to
-\url{http://www.gnu.org/s/gsl/}. For details on installing LAPACK
-library, please refer to \url{http://www.netlib.org/lapack/}.
+GNU gcc, as well as GSL and OpenBLAS libraries. A sample
+Makefile is provided along with the source code.
\newpage
@@ -334,7 +328,7 @@ genotypes and using BIMBAM files for phenotypes) will result in
unwanted errors. BIMBAM format is particularly useful for imputed
genotypes, as PLINK codes genotypes using 0/1/2, while BIMBAM can
accommodate any real values between 0 and 2 (and any real values if
-paired with ``-notsnp" option). In addition, to estimate variance
+paired with ``-notsnp'' option). In addition, to estimate variance
components using summary statistics, GEMMA requires two other input
files: one contains marginal z-scores and the other contains SNP
category.
@@ -809,7 +803,11 @@ chr rs ps n_mis n_obs allele1 allele0 af beta se
\end{verbatim}
%
-The 11 columns are: chromosome numbers, snp ids, base pair positions on the chromosome, number of missing individuals for a given snp, number of non-missing individuals for a given snp, minor allele, major allele, allele frequency, beta estimates, standard errors for beta, and $p$ values from the Wald test.
+The 11 columns are: chromosome numbers, snp ids, base pair positions
+on the chromosome, number of missing individuals for a given snp,
+number of non-missing individuals for a given snp, minor allele, major
+allele, allele frequency, beta estimates (additive effect), standard
+errors for beta, and $p$ values from the Wald test.
\subsection{Estimate Relatedness Matrix from Genotypes}
=====================================
guix.scm
=====================================
@@ -0,0 +1,61 @@
+;; To use this file to build HEAD of gemma:
+;;
+;; guix build -f guix.scm
+;;
+;; To get a development container (e.g., run in emacs shell).
+;;
+;; guix environment -C -l guix.scm
+
+(use-modules
+ ((guix licenses) #:prefix license:)
+ (guix gexp)
+ (guix packages)
+ (guix git-download)
+ (guix build-system meson)
+ (gnu packages algebra)
+ (gnu packages base)
+ (gnu packages compression)
+ (gnu packages bioinformatics)
+ (gnu packages build-tools)
+ (gnu packages curl)
+ (gnu packages llvm)
+ (gnu packages maths)
+ (gnu packages ninja)
+ (gnu packages parallel)
+ (gnu packages perl)
+ (gnu packages perl6)
+ (gnu packages pkg-config)
+ ;; (gnu packages shell) ;; for shunit2
+ (srfi srfi-1)
+ (ice-9 popen)
+ (ice-9 rdelim))
+
+(define %source-dir (dirname (current-filename)))
+
+(define %git-commit
+ (read-string (open-pipe "git show HEAD | head -1 | cut -d ' ' -f 2" OPEN_READ)))
+
+(define-public gemma-git
+ (package
+ (name "gemma-git")
+ (version (git-version "0.98.4" "HEAD" %git-commit))
+ (source (local-file %source-dir #:recursive? #t))
+ (build-system meson-build-system)
+ (inputs `(
+ ("gsl" ,gsl)
+ ;; ("shunit2" ,shunit2) ;; comes with gemma
+ ("openblas" ,openblas)
+ ("zlib" ,zlib)
+ ))
+ (native-inputs ; for running tests
+ `(("perl" ,perl)
+ ("which" ,which)
+ ))
+ (home-page "https://github.com/genetics-statistics")
+ (synopsis "Tool for genome-wide efficient mixed model association")
+ (description "Genome-wide Efficient Mixed Model Association (GEMMA)
+provides a standard linear mixed model resolver with application in
+genome-wide association studies (GWAS).")
+ (license license:gpl3)))
+
+gemma-git
=====================================
src/gemma.cpp
=====================================
@@ -2550,8 +2550,8 @@ void GEMMA::BatchRun(PARAM &cPar) {
// LMM or mvLMM or Eigen-Decomposition
if (cPar.a_mode == M_LMM1 || cPar.a_mode == M_LMM2 || cPar.a_mode == M_LMM3 ||
- cPar.a_mode == M_LMM4 || cPar.a_mode == M_LMM5 ||
- cPar.a_mode == M_EIGEN) { // Fit LMM or mvLMM or eigen
+ cPar.a_mode == M_LMM4 || cPar.a_mode == M_LMM5 || cPar.a_mode == M_LMM9 ||
+ cPar.a_mode == M_EIGEN ) { // Fit LMM or mvLMM or eigen
gsl_matrix *Y = gsl_matrix_safe_alloc(cPar.ni_test, cPar.n_ph);
enforce_msg(Y, "allocate Y"); // just to be sure
gsl_matrix *W = gsl_matrix_safe_alloc(Y->size1, cPar.n_cvt);
@@ -2565,6 +2565,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
gsl_vector *eval = gsl_vector_calloc(Y->size1);
gsl_vector *env = gsl_vector_safe_alloc(Y->size1);
gsl_vector *weight = gsl_vector_safe_alloc(Y->size1);
+ debug_msg("Started on LMM");
assert_issue(is_issue(26), UtY->data[0] == 0.0);
// set covariates matrix W and phenotype matrix Y
@@ -2578,6 +2579,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
if (!(cPar.file_kin).empty()) {
ReadFile_kin(cPar.file_kin, cPar.indicator_idv, cPar.mapID2num,
cPar.k_mode, cPar.error, G);
+ debug_msg("Read K/GRM file");
if (cPar.error == true) {
cout << "error! fail to read kinship/relatedness file. " << endl;
return;
@@ -2609,7 +2611,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
}
}
- // eigen-decomposition and calculate trace_G
+ // eigen-decomposition and calculate trace_G - main track
cout << "Start Eigen-Decomposition..." << endl;
time_start = clock();
@@ -2663,7 +2665,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
if (cPar.a_mode == M_EIGEN) {
cPar.WriteMatrix(U, "eigenU");
cPar.WriteVector(eval, "eigenD");
- } else if (!cPar.file_gene.empty()) {
+ } else if (!cPar.file_gene.empty()) { // Run with gene file
// calculate UtW and Uty
CalcUtX(U, W, UtW);
CalcUtX(U, Y, UtY);
@@ -2682,6 +2684,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
cLmm.WriteFiles();
cLmm.CopyToParam(cPar);
} else {
+ debug_msg("Main LMM track");
// calculate UtW and Uty
CalcUtX(U, W, UtW);
CalcUtX(U, Y, UtY);
@@ -2736,6 +2739,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
CalcPve(eval, UtW, &UtY_col.vector, cPar.l_remle_null, cPar.trace_G,
cPar.pve_null, cPar.pve_se_null);
+ debug_msg("main print summary");
cPar.PrintSummary();
// calculate and output residuals
@@ -2771,15 +2775,16 @@ void GEMMA::BatchRun(PARAM &cPar) {
gsl_vector_safe_free(u_hat);
gsl_vector_safe_free(e_hat);
gsl_vector_safe_free(y_hat);
- }
+ } // output residuals
}
// Fit LMM or mvLMM (w. LOCO)
- if (cPar.a_mode == 1 || cPar.a_mode == 2 || cPar.a_mode == 3 ||
- cPar.a_mode == 4) {
+ if (cPar.a_mode == M_LMM1 || cPar.a_mode == M_LMM2 || cPar.a_mode == M_LMM3 ||
+ cPar.a_mode == M_LMM4 || cPar.a_mode == M_LMM9) {
if (cPar.n_ph == 1) {
+ debug_msg("fit LMM (one phenotype)");
LMM cLmm;
- cLmm.CopyFromParam(cPar);
+ cLmm.CopyFromParam(cPar); // set parameters
// if (is_check_mode()) disable_segfpe(); // disable fast NaN checking for now
@@ -2811,8 +2816,9 @@ void GEMMA::BatchRun(PARAM &cPar) {
cLmm.WriteFiles();
cLmm.CopyToParam(cPar);
} else {
+ debug_msg("fit mvLMM (multiple phenotypes)");
MVLMM cMvlmm;
- cMvlmm.CopyFromParam(cPar);
+ cMvlmm.CopyFromParam(cPar); // set parameters
// if (is_check_mode()) disable_segfpe(); // disable fast NaN checking
@@ -3377,7 +3383,8 @@ void GEMMA::WriteLog(int argc, char **argv, PARAM &cPar) {
}
if (cPar.a_mode == 1 || cPar.a_mode == 2 || cPar.a_mode == 3 ||
- cPar.a_mode == 4 || cPar.a_mode == 5 || cPar.a_mode == 11 ||
+ cPar.a_mode == 4 || cPar.a_mode == 5 || cPar.a_mode == 9 ||
+ cPar.a_mode == 11 ||
cPar.a_mode == 12 || cPar.a_mode == 13) {
outfile << "## REMLE log-likelihood in the null model = "
<< cPar.logl_remle_H0 << endl;
@@ -3546,7 +3553,7 @@ void GEMMA::WriteLog(int argc, char **argv, PARAM &cPar) {
outfile << "## time on calculating UtX = " << cPar.time_UtX << " min "
<< endl;
}
- if ((cPar.a_mode >= 1 && cPar.a_mode <= 4) ||
+ if ((cPar.a_mode >= 1 && cPar.a_mode <= 4) || cPar.a_mode == 9 ||
(cPar.a_mode >= 51 && cPar.a_mode <= 54)) {
outfile << "## time on optimization = " << cPar.time_opt << " min "
<< endl;
=====================================
src/gemma.h
=====================================
@@ -41,8 +41,10 @@ using namespace std;
// gw: 72
enum M_MODE { M_LMM1=1, M_LMM2=2, M_LMM3=3, M_LMM4=4, M_LMM5=5,
- M_BSLMM5=15,
- M_KIN=21, M_KIN2=22, M_EIGEN=31 };
+ M_LMM9=9, // GeneNetwork mode
+ M_BSLMM5=15,
+ M_KIN=21, M_KIN2=22, M_EIGEN=31
+};
class GEMMA {
=====================================
src/lmm.cpp
=====================================
@@ -42,6 +42,7 @@
#include "gsl/gsl_vector.h"
#include "gzstream.h"
+#include "gemma.h"
#include "gemma_io.h"
#include "fastblas.h"
#include "lapack.h"
@@ -109,27 +110,28 @@ void LMM::WriteFiles() {
}
auto common_header = [&] () {
- if (a_mode != 2) {
+ if (a_mode != M_LMM2) {
outfile << "beta" << "\t";
outfile << "se" << "\t";
}
- if (!is_legacy_mode())
- outfile << "logl_H1" << "\t"; // we may make this an option
+ if (a_mode != M_LMM3 && a_mode != M_LMM9)
+ outfile << "logl_H1" << "\t";
switch(a_mode) {
- case 1:
+ case M_LMM1:
outfile << "l_remle" << "\t"
<< "p_wald" << endl;
break;
- case 2:
+ case M_LMM2:
+ case M_LMM9:
outfile << "l_mle" << "\t"
<< "p_lrt" << endl;
break;
- case 3:
+ case M_LMM3:
outfile << "p_score" << endl;
break;
- case 4:
+ case M_LMM4:
outfile << "l_remle" << "\t"
<< "l_mle" << "\t"
<< "p_wald" << "\t"
@@ -142,27 +144,28 @@ void LMM::WriteFiles() {
auto sumstats = [&] (SUMSTAT st) {
outfile << scientific << setprecision(6);
- if (a_mode != 2) {
+ if (a_mode != M_LMM2) {
outfile << st.beta << "\t";
outfile << st.se << "\t";
}
- if (!is_legacy_mode())
+ if (a_mode != M_LMM3 && a_mode != M_LMM9)
outfile << st.logl_H1 << "\t";
switch(a_mode) {
- case 1:
+ case M_LMM1:
outfile << st.lambda_remle << "\t"
<< st.p_wald << endl;
break;
- case 2:
+ case M_LMM2:
+ case M_LMM9:
outfile << st.lambda_mle << "\t"
<< st.p_lrt << endl;
break;
- case 3:
+ case M_LMM3:
outfile << st.p_score << endl;
break;
- case 4:
+ case M_LMM4:
outfile << st.lambda_remle << "\t"
<< st.lambda_mle << "\t"
<< st.p_wald << "\t"
@@ -173,6 +176,7 @@ void LMM::WriteFiles() {
};
+
if (!file_gene.empty()) {
outfile << "geneID" << "\t";
@@ -1114,7 +1118,7 @@ void LogRL_dev12(double l, void *params, double *dev1, double *dev2) {
return;
}
-void LMM::CalcRLWald(const double &l, const FUNC_PARAM ¶ms, double &beta,
+void LMM::CalcRLWald(const double l, const FUNC_PARAM ¶ms, double &beta,
double &se, double &p_wald) {
size_t n_cvt = params.n_cvt;
size_t n_index = (n_cvt + 2 + 1) * (n_cvt + 2) / 2;
@@ -1156,7 +1160,8 @@ void LMM::CalcRLWald(const double &l, const FUNC_PARAM ¶ms, double &beta,
return;
}
-void LMM::CalcRLScore(const double &l, const FUNC_PARAM ¶ms, double &beta,
+
+void LMM::CalcRLScore(const double l, const FUNC_PARAM ¶ms, double &beta,
double &se, double &p_score) {
size_t n_cvt = params.n_cvt;
size_t n_index = (n_cvt + 2 + 1) * (n_cvt + 2) / 2;
@@ -1197,7 +1202,6 @@ void LMM::CalcRLScore(const double &l, const FUNC_PARAM ¶ms, double &beta,
gsl_matrix_free(Pab);
gsl_vector_safe_free(Hi_eval);
gsl_vector_safe_free(v_temp);
- return;
}
void CalcUab(const gsl_matrix *UtW, const gsl_vector *Uty, gsl_matrix *Uab) {
@@ -1418,7 +1422,7 @@ void LMM::AnalyzeGene(const gsl_matrix *U, const gsl_vector *eval,
CalcUab(UtW, Uty, Uab);
FUNC_PARAM param0 = {false, ni_test, n_cvt, eval, Uab, ab, 0};
- if (a_mode == 2 || a_mode == 3 || a_mode == 4) {
+ if (a_mode == M_LMM2 || a_mode == M_LMM3 || a_mode == M_LMM4) {
CalcLambda('L', param0, l_min, l_max, n_region, l_H0, logl_H0);
}
@@ -1427,16 +1431,16 @@ void LMM::AnalyzeGene(const gsl_matrix *U, const gsl_vector *eval,
FUNC_PARAM param1 = {false, ni_test, n_cvt, eval, Uab, ab, 0};
// 3 is before 1.
- if (a_mode == 3 || a_mode == 4) {
+ if (a_mode == M_LMM3 || a_mode == M_LMM4) {
CalcRLScore(l_H0, param1, beta, se, p_score);
}
- if (a_mode == 1 || a_mode == 4) {
+ if (a_mode == M_LMM1 || a_mode == M_LMM4) {
CalcLambda('R', param1, l_min, l_max, n_region, lambda_remle, logl_H1);
CalcRLWald(lambda_remle, param1, beta, se, p_wald);
}
- if (a_mode == 2 || a_mode == 4) {
+ if (a_mode == M_LMM2 || a_mode == M_LMM4) {
CalcLambda('L', param1, l_min, l_max, n_region, lambda_mle, logl_H1);
p_lrt = gsl_cdf_chisq_Q(2.0 * (logl_H1 - logl_H0), 1);
}
@@ -1526,17 +1530,17 @@ void LMM::Analyze(std::function< SnpNameValues(size_t) >& fetch_snp,
double logl_H1 = 0.0;
// 3 is before 1.
- if (a_mode == 3 || a_mode == 4) {
+ if (a_mode == M_LMM3 || a_mode == M_LMM4 || a_mode == M_LMM9 ) {
CalcRLScore(l_mle_null, param1, beta, se, p_score);
}
- if (a_mode == 1 || a_mode == 4) {
+ if (a_mode == M_LMM1 || a_mode == M_LMM4) {
// for univariate a_mode is 1
CalcLambda('R', param1, l_min, l_max, n_region, lambda_remle, logl_H1);
CalcRLWald(lambda_remle, param1, beta, se, p_wald);
}
- if (a_mode == 2 || a_mode == 4) {
+ if (a_mode == M_LMM2 || a_mode == M_LMM9 || a_mode == M_LMM4) {
CalcLambda('L', param1, l_min, l_max, n_region, lambda_mle, logl_H1);
p_lrt = gsl_cdf_chisq_Q(2.0 * (logl_H1 - logl_mle_H0), 1);
}
@@ -1848,17 +1852,17 @@ void LMM::AnalyzePlink(const gsl_matrix *U, const gsl_vector *eval,
FUNC_PARAM param1 = {false, ni_test, n_cvt, eval, Uab, ab, 0};
// 3 is before 1, for beta.
- if (a_mode == 3 || a_mode == 4) {
+ if (a_mode == M_LMM3 || a_mode == M_LMM4) {
CalcRLScore(l_mle_null, param1, beta, se, p_score);
}
- if (a_mode == 1 || a_mode == 4) {
+ if (a_mode == M_LMM1 || a_mode == M_LMM4) {
CalcLambda('R', param1, l_min, l_max, n_region, lambda_remle,
logl_H1);
CalcRLWald(lambda_remle, param1, beta, se, p_wald);
}
- if (a_mode == 2 || a_mode == 4) {
+ if (a_mode == M_LMM2 || a_mode == M_LMM4) {
CalcLambda('L', param1, l_min, l_max, n_region, lambda_mle, logl_H1);
p_lrt = gsl_cdf_chisq_Q(2.0 * (logl_H1 - logl_mle_H0), 1);
}
=====================================
src/lmm.h
=====================================
@@ -118,9 +118,9 @@ public:
const gsl_vector *env);
void WriteFiles();
- void CalcRLWald(const double &lambda, const FUNC_PARAM ¶ms, double &beta,
+ void CalcRLWald(const double lambda, const FUNC_PARAM ¶ms, double &beta,
double &se, double &p_wald);
- void CalcRLScore(const double &l, const FUNC_PARAM ¶ms, double &beta,
+ void CalcRLScore(const double l, const FUNC_PARAM ¶ms, double &beta,
double &se, double &p_score);
};
=====================================
src/param.cpp
=====================================
@@ -107,7 +107,8 @@ PARAM::PARAM(void)
time_UtZ(0.0), time_opt(0.0), time_Omega(0.0) {}
PARAM::~PARAM() {
- gsl_rng_free(gsl_r);
+ if (gsl_r)
+ gsl_rng_free(gsl_r);
}
// Read files: obtain ns_total, ng_total, ns_test, ni_test.
@@ -510,6 +511,7 @@ void PARAM::CheckParam(void) {
error = true;
}
if (a_mode != 1 && a_mode != 2 && a_mode != 3 && a_mode != 4 && a_mode != 5 &&
+ a_mode != M_LMM9 &&
a_mode != 11 && a_mode != 12 && a_mode != 13 && a_mode != 14 &&
a_mode != 15 && a_mode != 21 && a_mode != 22 && a_mode != 25 &&
a_mode != 26 && a_mode != 27 && a_mode != 28 && a_mode != 31 &&
@@ -649,7 +651,7 @@ void PARAM::CheckParam(void) {
// Only LMM option (and one prediction option) can deal with
// multiple phenotypes and no gene expression files.
if (n_ph > 1 && a_mode != 1 && a_mode != 2 && a_mode != 3 && a_mode != 4 &&
- a_mode != 43) {
+ a_mode != 9 && a_mode != 43) {
cout << "error! the current analysis mode " << a_mode
<< " can not deal with multiple phenotypes." << endl;
error = true;
@@ -919,7 +921,7 @@ void PARAM::CheckParam(void) {
enforce_fexists(file_anno, "open file");
if (!loco.empty()) {
- enforce_msg((a_mode >= 1 && a_mode <= 4) || a_mode == 21 || a_mode == 22,
+ enforce_msg((a_mode >= 1 && a_mode <= 4) || a_mode == 9 || a_mode == 21 || a_mode == 22,
"LOCO only works with LMM and K");
// enforce_msg(file_bfile.empty(), "LOCO does not work with PLink (yet)");
enforce_msg(file_gxe.empty(), "LOCO does not support GXE (yet)");
@@ -947,7 +949,7 @@ void PARAM::CheckParam(void) {
}
if ((a_mode == 1 || a_mode == 2 || a_mode == 3 || a_mode == 4 ||
- a_mode == 5 || a_mode == 31) &&
+ a_mode == 5 || a_mode == 9 || a_mode == 31) &&
(file_kin.empty() && (file_ku.empty() || file_kd.empty()))) {
cout << "error! missing relatedness file. " << endl;
error = true;
=====================================
src/param.h
=====================================
@@ -213,7 +213,7 @@ public:
size_t n_mh; // # MH steps in each iter.
double geo_mean; // Mean of geometric dist.
long int randseed; // holds -seed parameter
- gsl_rng *gsl_r; // Track the randomizer
+ gsl_rng *gsl_r = NULL; // Track the randomizer
double trace_G;
HYPBSLMM cHyp_initial;
=====================================
src/version.h
=====================================
@@ -1,5 +1,5 @@
// version.h generated by GEMMA scripts/gen_version_info.sh
-#define GEMMA_VERSION "0.98.3"
-#define GEMMA_DATE "2020-11-28"
-#define GEMMA_YEAR "2020"
-#define GEMMA_PROFILE ""
+#define GEMMA_VERSION "0.98.4"
+#define GEMMA_DATE "2021-01-29"
+#define GEMMA_YEAR "2021"
+#define GEMMA_PROFILE "/gnu/store/8mkllydvkgfy6ydlrymrx8wj0dy1x6lm-profile"
=====================================
test/dev_test_suite.sh
=====================================
@@ -82,6 +82,23 @@ testBXDLMMLikelihoodRatio() {
assertEquals "3088458213" `perl -nle 'foreach $x (split(/\s+/,$_)) { $sum += sprintf("%.2f",(substr($x,,0,6))) } END { printf "%.0f",$sum }' $outfn`
}
+testBXDLMM9LikelihoodRatio() {
+ # Test for GeneNetwork output
+ outn=BXD_LMM_LR
+ $gemma $gemmaopts -g ../example/BXD_geno.txt.gz \
+ -p ../example/BXD_pheno.txt \
+ -c ../example/BXD_covariates2.txt \
+ -a ../example/BXD_snps.txt \
+ -k ./output/BXD.cXX.txt \
+ -lmm 9 -no-check -maf 0.1 \
+ -o $outn
+ assertEquals 0 $?
+
+ outfn=output/$outn.assoc.txt
+ assertEquals "80498" `wc -w < $outfn`
+ assertEquals "3088496565" `perl -nle 'foreach $x (split(/\s+/,$_)) { $sum += sprintf("%.2f",(substr($x,,0,6))) } END { printf "%.0f",$sum }' $outfn`
+}
+
testCenteredRelatednessMatrixissue188() {
outn=issue188
rm -f output/$outn.*
View it on GitLab: https://salsa.debian.org/med-team/gemma/-/compare/acf8bc6b6fa4a1c66c03107240052ed7a4298f04...8cf6fd2b71c9510451401a4767f86ad116bfdc82
--
View it on GitLab: https://salsa.debian.org/med-team/gemma/-/compare/acf8bc6b6fa4a1c66c03107240052ed7a4298f04...8cf6fd2b71c9510451401a4767f86ad116bfdc82
You're receiving this email because of your account on salsa.debian.org.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://alioth-lists.debian.net/pipermail/debian-med-commit/attachments/20210202/f5295d65/attachment-0001.html>
More information about the debian-med-commit
mailing list