HeritSeq

https://github.com/kechrislab/heritseq

Repository

README.Rmd

---
output: pdf_document
---

## Overview
The **HeritSeq** package provides heritability score analyses under linear mixed models (LMM) or generalized linear mixed models (GLMM) for count data motivated by high-throughput sequencing. It is applicable to counts with biological replicates. This package includes functions to:

- compute heritability score under LMM for a normalized/transformed dataset, and under negative binomial mixed models (NBMM) or compound Poisson mixed models (CPMM) for data without transformations.
- test presence of heritability under LMM, NBMM or CPMM.
- generate confidence intervals of estimated heritability score. 
- simulate synthetic high-throughput sequencing datasets using NBMM or CPMM. 

See [1] for model details and performance comparisons. 

[1]: Rudra, Pratyaydipta, W. Jenny Shi, Brian Vestal, Pamela H. Russell, Aaron Odell, Robin D. Dowell, Richard A. Radcliffe, Laura M. Saba, and Katerina Kechris. *Model based heritability scores for high-throughput sequencing data.* BMC bioinformatics 18, no. 1 (2017): 143. [pdf]( https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-017-1539-6 )

## Example dataset
The **HeritSeq** package includes an example high-throughput sequencing dataset 
called *simData*. This dataset was generated based on a recombinant inbred mice panel miRNA sequencing counts. It is an 100 by 175 matrix, containing 100 features and 175 samples. The total number of strains is 59, and the strain labels are recorded 
by the variable *strains*. 

## Installation
The package requires **R** version >=3.2.3.

### Dependencies
The installation of dependencies only needs to be done once, if at all.
Instructions for installing the dependencies:
```{r, eval = FALSE}
install.packages("lme4", repos="http://cran.r-project.org")
install.packages("cplm", repos="http://cran.r-project.org")
install.packages("pbapply", repos="http://cran.r-project.org")
install.packages("statmod", repos="http://cran.r-project.org")
install.packages("MASS", repos="http://cran.r-project.org")

install.packages("glmmADMB", 
                repos=c("http://glmmadmb.r-forge.r-project.org/repos",
                        getOption("repos")),
                type="source")

source("https://bioconductor.org/biocLite.R")
biocLite("DESeq2")
biocLite("SummarizedExperiment")
```

For some systems an [alternative installation]( http://glmmadmb.r-forge.r-project.org/ ) might be needed for the glmmADMB package.  

Since R for OS X Maverick (and latter versions) was compiled using gfortran-4.8, Mac users might receive a “-lgfortran” error when installing the cplm package. This problem and its solution are stated [here]( http://thecoatlessprofessor.com/programming/rcpp-rcpparmadillo-and-os-x-mavericks-lgfortran-and-lquadmath-error/).


### Direct installation from CRAN
```{r, eval = FALSE}
install.packages("HeritSeq", repos="http://cran.r-project.org")
```


### Tarball installation
Download the [package tarball](https://github.com/KechrisLab/HeritSeq/blob/master/package/HeritSeq_1.0.1.tar.gz) and install the package:

```{r, eval = FALSE}
install.packages("HeritSeq_1.0.1.tar.gz",repos=NULL,type="source")
library("HeritSeq")
```
```{r, echo = FALSE}
setwd("~/Documents/heritseq/")
load("data/heritseq_example.RData")
```

## Estimate heritability scores
We will use *simData* to illustrate the procedure of heritability estimation 
under different models. 

Before fitting any model, make sure that the input dataset has been adjusted for 
library sizes and batch bias. The dataset *simData* is post such adjustment, we 
therefore omit the process. 

### NBMM
Under NBMM, an observed number of reads aligned to feature/gene $g$, $Y_{gsr}$, follows a negative binomial distribution with mean $\mu_{gs}$ and variance $\mu_{gs}+\phi_{g} \mu_{gs}^2$, where $\phi_{g}$ is the dispersion parameter for feature/gene $g$, shared across strains. The generalized linear model uses a $\log$-link: 
$\log(\mu_{gs}) = \alpha_{g}+ b_{gs}, \;\;b_{gs}\sim N(0, \sigma^{2}_{g}).$

The corresponding heritability score, aka Variance Partition Coefficient (VPC),
is $\frac{e^{\sigma^2_{g}} - 1}{e^{\sigma^2_{g}} - 1 + \phi_{g} e^{\sigma^2_{g}} + e^{-\alpha_g-\sigma^2_{g}/2} }.$

Compute VPC for all features using NBMM:
```{r, eval = FALSE}
library(glmmADMB)
result.nb <- fit.NB(CountMatrix = simData, Strains = strains, test = FALSE)
vpc.nb <- computeVPC.NB(para = result.nb[[1]])
```
The function *fit.NB( )* returns a list with two objects. The first object is a $G \times 3$ matrix indicating the fitted parameters for each feature, where $G$ is the total number of features/genes. The columns are ordered by the fitted parameters $\alpha_{g}, \sigma^{2}_{g}, \phi_{g}$. The second object provides the p-value for testing the presence of heritability if test = TRUE; it returns NULL otherwise (default).

The function *computeVPC.NB( )* takes in the list of NBMM parameters and outputs the corresponding VPC values. 

### CPMM 
For a CP random variable $Y_{gsr}$ with mean $\mu_{gs}$, its variance can be expressed as $\phi_{g}\mu_{gs}^{p_{g}}$, for some tweedie parameter $1
      

Owner

GitHub Events

Committers

Last synced: almost 3 years ago

Issues and Pull Requests

Last synced: about 1 year ago

Packages