---

output: github_document

---

```{r, include = FALSE}

knitr::opts_chunk$set(

  collapse = TRUE,

  comment = "#>"

)

```

# GPROB <img src="man/figures/gprob.svg" width="181px" align="right" />

Multiple diseases can present with similar initial symptoms, making it

difficult to clinically differentiate between these conditions. GPROB uses

patients' genetic information to help prioritize a diagnosis. This genetic

diagnostic tool can be applied to any situation with phenotypically similar

diseases with different underlying genetics.

<!-- badges: start -->

[![R build status](https://github.com/immunogenomics/GPROB/workflows/R-CMD-check/badge.svg)](https://github.com/immunogenomics/GPROB/actions)

<!-- badges: end -->

## Citation

Please cite:

- Knevel, R. et al. [Using genetics to prioritize diagnoses for rheumatology

  outpatients with inflammatory arthritis.][1] Sci. Transl. Med. 12, (2020)

[1]: http://dx.doi.org/10.1126/scitranslmed.aay1548

## License

Please see the [LICENSE] file for details. [Contact us] for other licensing

options.

[LICENSE]: LICENSE

[Contact us]: mailto:soumya@broadinstitute.org

## Installation

Install and load the GPROB R package.

```{r, eval = FALSE}

devtools::install_github("immunogenomics/GPROB")

library(GPROB)

```

## Synopsis

GPROB estimates the probability that each individual has a given phenotype.

We need three inputs:

- Population prevalences of the phenotypes of interest.

- Odds ratios for SNP associations with the phenotypes.

- SNP genotypes (0, 1, 2) for each individual.

### Example

Let's use a small example with artificial data to learn how to use GPROB.

Suppose we have 10 patients, and we know of 7 single nucleotide polymorphisms

(SNPs) associated with rheumatoid arthritis (RA) or systemic lupus

erythematosus (SLE).

#### Prevalence

First, we should find out the prevalence of RA and SLE in the population that

is representative of our patients.

```{r}

prevalence <- c("RA" = 0.001, "SLE" = 0.001)

```

#### Odds Ratios

Next, we need to obtain the odds ratios (ORs) from published genome-wide

association studies (GWAS). We should be careful to note which alleles are

associated with the phenotype to compute the risk in the correct direction.

```{r}

or <- read.delim(

  sep = "",

  row.names = 1,

  text = "

snp  RA SLE

SNP1 1.0 0.4

SNP2 1.0 0.9

SNP3 1.0 1.3

SNP4 0.4 1.6

SNP5 0.9 1.0

SNP6 1.3 1.0

SNP7 1.6 1.0

")

or <- as.matrix(or)

```

#### Genotypes

Finally, we need the genotype data for each of our 10 patients. Here, the data

is coded in the form (0, 1, 2) to indicate the number of copies of the risk

allele.

```{r}

geno <- read.delim(

  sep = "",

  row.names = 1,

  text = "

id SNP1 SNP2 SNP3 SNP4 SNP5 SNP6

 1    0    1    0    2    1    0

 2    0    0    1    0    2    2

 3    1    0    1    1    0    2

 4    1    1    0    2    0    0

 5    0    1    1    1    1    0

 6    0    0    1    3    0    2

 7    2    2    2    2    2    2

 8    1    2    0    2    1    1

 9    0    2    1   NA    1    2

10    1    0    2    2    2    0

")

geno <- as.matrix(geno)

```

#### Dealing with missing or invalid data

Before we run the `GPROB()` function, we need to deal with invalid and missing

data.

We remove individuals who have `NA` for any SNP:

```{r}

ix <- apply(geno, 1, function(x) !any(is.na(x)))

geno <- geno[ix,]

```

We remove individuals who have invalid allele counts:

```{r}

ix <- apply(geno, 1, function(x) !any(x < 0 | x > 2))

geno <- geno[ix,]

```

And we make sure that we use the same SNPs in the `or` and `geno` matrices:

```{r}

or <- or[colnames(geno),]

```

#### Run GPROB

Then we can run the GPROB function to estimate probabilities:

```{r}

library(GPROB)

res <- GPROB(prevalence, or, geno)

res

```

In this example, we might interpret the numbers as follows:

- Individual 2 has RA with probability 0.003, given individual genetic risk

  factors, disease prevalence, and the number of patients used in genetic risk

  score calculations.

- Individual 2 has RA with probability 0.77, conditional on the additional

  assumption that individual 2 has either RA or SLE.

## Calculations, step by step

Let's go through each step of GPROB to understand how how it works.

The genetic risk score <i>S<sub>ki</sub></i> of individual *i* for disease *k* is defined as:

<p align="center">

<img src="https://latex.codecogs.com/svg.latex?\Large&space;S_{ki}=\sum_{j}{\beta_{kj}x_{ij}}"/>

</p>

where:

- <i>x<sub>ij</sub></i> is the number of risk alleles of SNP *j* in individual *i*

- <i>β<sub>kj</sub></i> is the log odds ratio for SNP *j* reported in a genome-wide

  association study (GWAS) for disease *k*

<table><tr><td>

<b>Note:</b> We might want to consider shrinking the risk by some factor (e.g.

0.5) to correct for possible overestimation of the effect sizes due to

publication bias. In other words, consider running <code>geno <- 0.5 *

geno</code>.

</td></tr></table>

```{r}

risk <- geno %*% log(or)

risk

```

The known prevalence <i>V<sub>k</sub></i> of each disease in the general population:

```{r}

prevalence

```

We can calculate the population level probability <i>P<sub>ki</sub></i> that each individual

has the disease.

<p align="center">

<img src="https://latex.codecogs.com/svg.latex?\Large&space;P_{ki}=\frac{1}{1+\exp{(S_{ki}-\alpha_k)}}"/>

</p>

We find <i>α<sub>k</sub></i> for each disease *k* by minimizing

<i>(P&#773;<sub>k</sub> - V<sub>k</sub>)<sup>2</sup></i>. This ensures that the

mean probability <i>P&#773;<sub>k</sub></i> across individuals is equal to the

known prevalence <i>V<sub>k</sub></i> of the disease in the population.

```{r}

# @param alpha A constant that we choose manually.

# @param risk A vector of risk scores for individuals.

# @returns A vector of probabilities for each individual.

prob <- function(alpha, risk) {

  1 / (

    1 + exp(alpha - risk)

  )

}

alpha <- sapply(seq(ncol(risk)), function(i) {

  o <- optimize(

    f        = function(alpha, risk, prevalence) {

      ( mean(prob(alpha, risk)) - prevalence ) ^ 2

    },

    interval = c(-100, 100),

    risk = risk[,i],

    prevalence = prevalence[i]

  )

  o$minimum

})

alpha

```

Now that we have computed alpha, we can compute the population-level

probabilities of disease for each individual.

```{r}

# population-level disease probability

p <- sapply(seq_along(alpha), function(i) prob(alpha[i], risk[,i]))

p

```

Next we assume that each individual has one of the diseases:

<p align="center">

<img src="https://latex.codecogs.com/svg.latex?\Large&space;\text{Pr}(Y_k=1|(\textstyle\sum_k{Y_k})=1)"/>

</p>

Then, we calculate the conditional probability <i>C<sub>ki</sub></i> of each

disease *k*: 

<p align="center">

<img src="https://latex.codecogs.com/svg.latex?\Large&space;C_{ki}=\frac{P_{ki}}{\sum_k{P_{ki}}}"/>

</p>

```{r}

# patient-level disease probability

cp <- p / rowSums(p)

cp

```