denovolyzeR intro

Introduction

This package provides functions to analyse de novo genetic variants using the statistical framework described in Samocha et al (2014) Nature Genetics 10.1038/ng.3050. This vignette demonstrates the usage of the denovolyzeR package to recapitulate the analyses described in this paper:

1. What is the overall burden of de novo variation in the study cohort: are there more de novos than expected? As well as looking at the total burden of de novos, results are returned for different classes of variant, such as loss-of-function (lof), missense (mis) and synonymous (syn).
2. Do de novo variants cluster in specific genes: are there more genes containing multiple de novos than expected?
3. Are there any individual genes that contain more de novos than expected?

If using the package, please cite Ware et al (2015) Curr Protoc Hum Genet. 10.1002/0471142905.hg0725s87.

Overall de novo burden

First we want to know whether there are more de novos than expected, using the denovolyzeByClass() function. These variants were obtained by sequencing 1,078 cases, so we use nsamples=1078.

denovolyzeByClass(genes=autismDeNovos$gene,
                  classes=autismDeNovos$class,
                  nsamples=1078)

##   class observed expected enrichment   pValue
## 1   syn      254    302.3      0.840 0.998000
## 2   mis      655    679.0      0.965 0.826000
## 3   lof      131     94.3      1.390 0.000199
## 4  prot      786    773.2      1.020 0.328000
## 5   all     1040   1075.5      0.967 0.864000

The total number of de novos is almost exactly as our model predicts. However, we see a statistically significant excess of LOF variants in this population.

Genes containing multiple de novos

Next, we look to see if the total number of genes that contain more than one de novo is greater than expected, using the denovolyzeMultiHits() function.

denovolyzeMultiHits(genes=autismDeNovos$gene,
                    classes=autismDeNovos$class,
                    nsamples=1078)

##   class obs expMean expMax pValue nVars
## 1   syn   3     3.2      9   0.66   254
## 2   mis  31    20.5     30   0.00   655
## 3   lof   5     0.9      4   0.00   131
## 4  prot  43    27.7     43   0.01   786
## 5   all  66    47.9     65   0.00  1040

obs = the number of genes in our dataset with >1 de novo variant
expMean = the expected number of genes containing >1 de novo: an average obtained by permutation
expMax = the maximum number of genes containing >1 de novo in nperms permutations (default nperms=100)
pValue = an empirical p value
nVars = the total number of de novo variants in each class
Note that the number of observed genes with >1 protein-altering variant does not equal the number of genes with >1 lof + number of genes with >1 missense, as genes containing 1 lof + 1 missense will only be counted as “multihits” in the combined analysis.

Here it looks like there may be an excess of genes with >1 lof variant, >1 missense, and >1 protein-altering variant. We will want to increase the number if permutations here to get a handle on our level of significance.

denovolyzeMultiHits(genes=autismDeNovos$gene,
                    classes=autismDeNovos$class,
                    nsamples=1078,
                    nperms=1000)

##   class obs expMean expMax pValue nVars
## 1   syn   3     3.3     10  0.646   254
## 2   mis  31    20.1     33  0.011   655
## 3   lof   5     1.1      6  0.006   131
## 4  prot  43    28.2     45  0.002   786
## 5   all  66    47.5     73  0.006  1040

There is another important option here. The expected number of genes containing >1 hit is obtained by permutation: Given n de novo variants, how many genes contain >1 de novo? There are two options for selecting n: by default it is derived from your data: e.g. in the example above autismDeNovos contains 131 lof variants, so this is the number used in the permutation. This is controlled by the default parameter nVars="actual"

sum(autismDeNovos$class %in% c("frameshift","non","splice"))

## [1] 131

This is a conservative approach, addressing the question: “given the number of variants in our dataset, do we see more genes with >1 variant than expected?”

An alternative approach simply asks whether there are more genes with >1 variant than our de novo model predicts. This is accessed by setting nVars="expected".

denovolyzeMultiHits(genes=autismDeNovos$gene,
                    classes=autismDeNovos$class,
                    nsamples=1078,
                    nperms=1000,
                    nVars="expected")

##   class obs expMean expMax pValue      nVars
## 1   syn   3     4.5     12  0.848  302.27646
## 2   mis  31    21.3     35  0.020  678.98816
## 3   lof   5     0.5      4  0.000   94.26139
## 4  prot  43    27.1     44  0.002  773.24955
## 5   all  66    50.6     73  0.011 1075.52601

Do any individual genes contain more de novos than expected

We see 5 genes containing >1 de novo lof variant. This is more than expected, but are any of these genes individually significant? We can denovolyzeByGene() to find out.

By default this function compares the number of LOF variants against expectation for each gene, and then the total number of protein-altering variants (LOF + missense). It can also be configured to return other classes if relevant.

head(
denovolyzeByGene(genes=autismDeNovos$gene,
                 classes=autismDeNovos$class,
                 nsamples=1078)
  )

##         gene lof_observed lof_expected lof_pValue prot_observed
## 239   DYRK1A            3            0   2.69e-08             3
## 706    SCN2A            3            0   1.83e-06             5
## 141     CHD8            3            0   7.19e-07             4
## 674     RFX8            0            0   1.00e+00             2
## 796 SUV420H1            1            0   6.37e-03             3
## 618     POGZ            2            0   1.23e-04             2
##     prot_expected prot_pValue
## 239           0.1    2.77e-05
## 706           0.1    3.70e-07
## 141           0.2    2.44e-05
## 674           0.0    2.34e-05
## 796           0.1    3.17e-05
## 618           0.1    5.07e-03

Several genes meet statistical significance after correcting for multiple testing. Default options apply two tests across 19618 genes, so a Bonferroni corrected p-value threshold at \(\alpha\) = 0.05 would be \(1.3\times 10^{-6}\).

Geneset analysis

The analyses presented so far have been exome-wide. It may be appropriate to restrict analyses to a geneset of interest - for example, it may be relevant to examine the burden of de novo variation in a pathway of interest, or initial variant detection may have been restricted to a set of candidate genes (rather than whole exome sequencing). All of the above funtions can be targeted to a subset of genes using the includeGenes argument.

The package includes as an example a list of 837 genes that interact with the fragile X mental retardation protein (FMRP). Is this geneset enriched for de novos, and recurrent de novos, in our autism trios?

nrow(fmrpGenes); head(fmrpGenes)

## [1] 837

##            ensgID          enstID hgncID hgncSymbol geneName
## 1 ENSG00000142599 ENST00000337907   9965       RERE     RERE
## 2 ENSG00000149527 ENST00000449969  29037      PLCH2    PLCH2
## 3 ENSG00000078369 ENST00000378609   4396       GNB1     GNB1
## 4 ENSG00000157933 ENST00000378536  10896        SKI      SKI
## 5 ENSG00000171735 ENST00000303635  18806     CAMTA1   CAMTA1
## 6 ENSG00000188157 ENST00000379370    329       AGRN     AGRN

denovolyzeByClass(genes=autismDeNovos$gene,
                  classes=autismDeNovos$class,
                  nsamples=1078,
                  includeGenes=fmrpGenes$geneName)

##   class observed expected enrichment   pValue
## 1   syn       28     33.6      0.835 8.53e-01
## 2   mis       83     74.4      1.110 1.74e-01
## 3   lof       32      9.1      3.500 3.18e-09
## 4  prot      115     83.6      1.380 6.47e-04
## 5   all      143    117.1      1.220 1.13e-02

denovolyzeMultiHits(genes=autismDeNovos$gene,
                    classes=autismDeNovos$class,
                    nsamples=1078,
                    nperms=1000,
                    includeGenes=fmrpGenes$geneName)

## Warning in denovolyzeMultiHits(genes = autismDeNovos$gene, classes =
## autismDeNovos$class, : De novo list includes 897 genes not specified for
## inclusion. These will not be analysed.

##   class obs expMean expMax pValue nVars
## 1   syn   0     0.6      4  1.000    28
## 2   mis   9     5.1     11  0.061    83
## 3   lof   3     1.0      6  0.061    32
## 4  prot  15     9.4     18  0.028   115
## 5   all  20    14.3     25  0.045   143