---
title: "GWAS solution overview"
---

# Principal component analysis

```{r}
library("MASS")
data(biopsy)
names(biopsy)
predictors <- biopsy[complete.cases(biopsy),2:10]
fit <- prcomp(predictors, scale=TRUE)
summary(fit)
fit
```

1.  The `rotation` element is the interesting part. 

2.  (and 3 + 4)

```{r}
plot(fit)
biplot(fit)
```

Only one real PC needed (from the elbow plot).  The biplot shows that one variable drives PC2 while all of them have essentially the same effect on PC1. Not surprising as they are all more and more critical indicators of tumor features.


5.  Fitting using a formula model

```{r}
fit2 <- prcomp(~ V1 + V2 + V3 + V3 + V5 + V6 + V7 + V8 + V9, 
              data=biopsy, na.action=na.omit)
print(fit2)
```


6.  

```{r}
cor(predictors)
```




# Principal component regression


We already computed the PCs above so we can just use those. The computed values can either be extracted as the `x` element or by using the `predict()` function.

1.  The outcome is binary so we will use a logistic regression model. Trying with the first four PCs.

    ```{r}
    y <- biopsy$class[complete.cases(biopsy)]
    pcr <- glm(y ~ fit$x[,1] + fit$x[,2] + fit$x[,3] + fit$x[,4],
               family=binomial)
    summary(pcr)
    ```
    

2.  PC 1 provides a lot of info about the outcome, and the coefficient is negative which means that at increase for the PC value will lower the odds of a malignant outcome. However, recall that the rotation matrix for the first principal component only had negative weights so when that is combined we get that an increase in any predictor will result in an increase in log odds of a malignant tumor.  
    We also see that even though the PCs are ranked from highest variation to lowest variation their p-values are not necessarily ranked in the same order.
    
3.  Primarily prediction as the components are difficult to interpret.

4.  
    ```{r}
    pcr2 <- glm(y ~ fit$x[,1] + fit$x[,2],
               family=binomial)
    summary(pcr2)
    ```
    
    We see that there is hardly any difference for the predictors. The reason why the predictors are not completely identical to the fit above where we had 4 predictors is because of the non-identity link function.
    
5.  We will need to see how the predictability improves with number of components.     



# GWAS analysis

First we get the data

```{r cache=TRUE}
options(timeout = 5*60)
load(url("https://www.biostatistics.dk/teaching/bioinformatics/data/gwasgt.rda"))
load(url("https://www.biostatistics.dk/teaching/bioinformatics/data/gwaspt.rda"))
```

1.  The number of columns in the `genotype` dataframe shows the number of SNPs available.

```{r}
dim(genotypes)
```

2.  If we had just a single SNP available we would start by doing a multiple linear regression analysis since the phenotype outcome (BMI) is quantitative. That would also allow us to adjust for the other variables that are available in the `phenotypes` dataframe.

3.  For example (just using the 1st column in `genotypes`) as my one SNP:

   ```{r}
   m1 <- lm(BMI ~ genotypes[,1] + gender + age, data=phenotypes) 
   summary(m1)
   ```
   
4.  In the following we are discarding `age`  and `gender` since the function we will be using only takes a single matrix of predictors. We could include their effects indirectly by subtracting their respective (estimated) effects from `BMI` and analysing the residuals as outcomes instead.

    ```{r warning=FALSE, message=FALSE}
    library("MESS")
    full <- mfastLmCpp(phenotypes$BMI, genotypes)
    ```
   
5.  For the Manhattan plot we use that the 3rd column in `full` is the t test statistic *and* we know there are 1324 rows.

    ```{r}
    logpvals <- -log10(2*pt(-abs(full[,3]), df=1324-2)) 
    plot(logpvals, col=(seq(32019)>11283)+1)
    ```
    
6.  Now running lasso

    ```{r warning=FALSE, message=FALSE}
    library("glmnet")
    l1 <- glmnet(genotypes, phenotypes$BMI)
    l12 <- cv.glmnet(genotypes, phenotypes$BMI) # Get the lambda parameter
    ```
    
    We can check if the predictors selected are the same for the two methods.
    
    ```{r}
    estcoefs <- coef(l1, s=l12$lambda.1se)
    pick <- rownames(estcoefs)[as.numeric(estcoefs) != 0]  # Remove the intercept
    pick
    ```
    
    The previous analysis produced this set
    
    ```{r}
    seq(NROW(full))[logpvals > -log10(0.05/NCOL(genotypes))]
    ```
    
7.  Same as above with `alpha` parameter set.
8.  For computing the genomic controls we find the median of the $p$-values

    ```{r}
    median(full[,3]^2)
    ```
    Since this number is klower than the default level of 0.456 there is no need for correcting for genomic controls.
    
    
# Population admixture

1.  Because hidden population admixture may act as a confounder between allele frequencies and the outcome and may mask the true relationship.

2.  Using the `onlinePCA` package for fast computation of the first principal components.

    ```{r, messages=FALSE, warnings=FALSE}
    library("onlinePCA")
    admixturepca <- batchpca(genotypes, q=1)
    ```
    
    Let's extract the effect of the PC 
    
    ```{r}
    am1 <- lm(phenotypes$BMI ~ admixturepca$vectors)
    coef(am1)
    summary(am1)
    adjustedBMI <- phenotypes$BMI-coef(am1)[2]*admixturepca$vectors  # Subtracting the effect
    full2 <- mfastLmCpp(adjustedBMI, genotypes)
    logpvals <- -log10(2*pt(-abs(full2[,3]), df=1324-3)) 
    plot(logpvals, col=(seq(32019)>11283)+1)
    ```


3.  If at all possible we should use SNPs that are not our primary SNPs of interest.




# Hertability

1.  The minor allele frequencies are

    ```{r}
    MAF <- colMeans(genotypes)/2 # Minor allele frequencies
    MAF[MAF>.5] <- 1-MAF[MAF>.5]
    ```
    
2.  Compute the SDs for each column

    ```{r}
    sds <- sqrt(2*MAF*(1-MAF))
    ```

3.  Running the code

    ```{r warnings=FALSE, messages=FALSE, cache=TRUE}
    W <- t((t(genotypes) - 2*MAF)/sqrt(2*MAF*(1-MAF)))
    y <- phenotypes$BMI
    N <- length(y)
    A <- tcrossprod(W) / ncol(genotypes)
    library("coxme")
    id <- 1:N
    res <- lmekin(y ~ 1 + (1|id), varlist=list(id=A))
    res
    ```


4.  We are computing the SNP heritability. The `A` matrix contains the correlation structure derived from the SNPs.

5.  The SNP heritability is 

    ```{r}
    0.2276165 / (0.2276165 + 1.114715^2)
    ```
    
    
    
