cca_networks.Rmd

---
title: "CCA network analysis of NPS/AD"
subtitle: 'CCA networks with decorrelate'
author: "Developed by [Gabriel Hoffman](http://gabrielhoffman.github.io/)"
date: "Run on `r Sys.time()`"
documentclass: article
output: 
  html_document:
  toc: true
  smart: false
vignette: >
  %\VignetteIndexEntry{Decorrelate}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
  %\usepackage[utf8]{inputenc}
---


<!--- 

cd /hpc/users/hoffmg01/work/nps_ad
ml python #udunits proj gdal/2.4.1 geos unixODBC giflib tiff
R
# rm -rf analysis_1_cache/

system("ml git; git pull")
rmarkdown::render("cca_networks.Rmd");


# https://hoffmg01.u.hpc.mssm.edu/nps_ad/





--->



```{r setup, include=FALSE}
knitr::opts_chunk$set(
  echo = FALSE,
  warning=FALSE,
  message=FALSE,
  error = FALSE,
  tidy = FALSE,
  dev = c("png", "pdf"),
  cache = TRUE,
  cache.lazy = FALSE)
```

```{r load.packages, cache=FALSE}
# Use cache=FALSE so that package are fully loaded each time
# This ensures that forks within mclapply() have these loaded
# Othewise, mclapply() not have access to these libraries and will fail 
#   unless the libraries are manually loaded within each fork
suppressPackageStartupMessages({
library(SingleCellExperiment)
library(dreamlet)
library(decorrelate)
library(variancePartition)
library(RNOmni)
library(Rfast)
library(decorate)
})
```




```{r load.data}
# read pseudobulk
pbObj = readRDS("/sc/arion/projects/CommonMind/hoffman/NPS-AD/work/pbObj_v2.RDS")

# read dreamlet fit
fit.dl = readRDS("/sc/arion/projects/CommonMind/hoffman/NPS-AD/work/fit.dl.RDS")

# compute residuals and average replicates
residList = lapply(fit.dl, function(fit){

  # get residuals
  res = residuals(fit)

  # extract metadata for these samples
  md = metadata(pbObj)[[1]]
  ids = intersect(colnames(res), md$sample_id)
  df = md[match(ids, md$sample_id),]
  df$group_id = droplevels(df$group_id)

  # Compute mean for each donor
  resAgg = lapply( levels(df$group_id), function(gid){

    idx = which(gid == df$group_id)

    rowMeans(res[,idx,drop=FALSE])
  })
  resAgg = do.call(cbind, resAgg)
  colnames(resAgg) = levels(df$group_id)

  resAgg
  })
names(residList) = names(fit.dl)
```



```{r outliers, eval=FALSE}
scoreList = lapply(residList, function(res){

  zScores = t(scale(t(res)))
})
saveRDS(scoreList, file="scoreList.RDS")

# scp sklar1:/sc/arion/work/hoffmg01/nps_ad/scoreList.RDS ~/Downloads
```
```{r eval.outliers, eval=FALSE}
library(gtools)
library(Rfast)
library(org.Hs.eg.db)

scoreList = readRDS("~/Downloads/scoreList.RDS")

# res = scoreList[[1]]
library(org.Hs.eg.db)
library(gtools)

geneExpr = lapply( names(scoreList), function(CT){
  cat(CT)

  geneExpr = scoreList[[CT]]

  # df = select(org.Hs.eg.db, keys = rownames(geneExpr), columns="CHR", keytype="ENSEMBL")
  # keepGenes = unique(df$SYMBOL[df$CHR %in% 1:22])
  # geneExpr = geneExpr[keepGenes,]

  rownames(geneExpr) = paste0(CT, '_', rownames(geneExpr))
  data.frame(Symbol = rownames(geneExpr), geneExpr)
  })
geneExpr = do.call(smartbind, geneExpr)
rownames(geneExpr) = geneExpr$Symbol
geneExpr = as.matrix(geneExpr[,-1])


# se per sample
apply(geneExpr, 2, sd, na.rm=TRUE)

# count outliers per sample
count = apply(geneExpr, 2, function(x) sum(x>0, na.rm=TRUE))
mean(count)
var(count)



g = which.max(rowMaxs(geneExpr, value=TRUE))
id = which.max(geneExpr[g,])
geneExpr[g,id]

rownames(geneExpr)[g]
hist(geneExpr[g,])


C = cor( t(geneExpr[g,,drop=FALSE]), t(geneExpr), use="pair")

hist(C[,-g], xlim=c(-1,1))

median(C[,-g]^2)

max(C[,-g]^2)
C[,order(C)[1:5]]

```





```{r cca.analysis.functions}
pairwise_cca = function(exprList, verbose=TRUE){

  lvls = names(exprList)

  idx_combn = combn(length(lvls), 2)

  df = apply(idx_combn, 2, function(x){

    if(verbose) cat("\r", x, '   ')

    CT1 = names(exprList)[x[1]]
    CT2 = names(exprList)[x[2]]

    X = exprList[[CT1]]
    Y = exprList[[CT2]]

    ids = intersect(colnames(X), colnames(Y))

    if( length(ids) < 15){
      return(NULL)
    }

    X_sub = scale(t(X[,ids,drop=FALSE]))
    Y_sub = scale(t(Y[,ids,drop=FALSE]))

    group = c(rep(0, nrow(X_sub)), rep(1, nrow(Y_sub)))
    keep = intersect(colnames(X_sub), colnames(Y_sub))
    Z = rbind(X_sub[,keep], Y_sub[,keep]) 

    fit = boxM_fast(Z, group)

    fit.cca = cca(X_sub,Y_sub)

    data.frame(CT1 = CT1, CT2 = CT2, N = nrow(X_sub), lambda.x = fit.cca$lambdas[1],
      lambda.y = fit.cca$lambdas[2], cramer.V = fit.cca$cramer.V, cor1 = fit.cca$cor[1], rho1 = fit.cca$rho.mod[1], ri.x=sum(fit.cca$x.ri), ri.y=sum(fit.cca$y.ri), n.x=ncol(X_sub), n.y=ncol(Y_sub), BoxM.p = fit$pval, BoxM.ld1 = fit$logdet[1], BoxM.ld2 = fit$logdet[2], n.z=ncol(Z)) 
  })

  if( !is.null(df) ){
    df = do.call(rbind, df)
  }

  df
}

dfToMatrix = function(df, V1, V2, value){

  lvls = unique(c(df[,V1], df[,V2]))
  C = diag(1, length(lvls))
  C[C!=1] = NA
  rownames(C) = lvls
  colnames(C) = lvls

  for(i in 1:nrow(df) ){
    C[df[,V1][i], df[,V2][i]] = df[[value]][i]
    C[df[,V2][i], df[,V1][i]] = df[[value]][i]
  }

  C
}
```

```{r cca.analysis}
df_cca = pairwise_cca( residList )

```

```{r}
with(df_cca, plot(rho, rho.emp))
```

```{r plots}

C = dfToMatrix(df_cca, "CT1", "CT2", "rho.emp")

plotCorrMatrix(C^2)
```

```{r zenith, message=FALSE}
library(zenith)
library(GSEABase) 
library(ggplot2)

gsGO = get_GeneOntology(c("BP", "MF"), "SYMBOL")

len = sapply(gsGO, function(x) length(geneIds(x)))
idx = which(len >=10 & len <1000)
gsGO = gsGO[idx]
length(gsGO)

exclude = c('Endothelial', 'Unknown', 'OPC', 'Mural',"GABAergic PVALB+", "GABAergic SST+", "GABAergic VIP+")
keep = ! names(residList) %in% exclude

res = lapply(gsGO, function(gs){

  i = match(setName(gs), names(gsGO))
  message("\r", i, '   ', appendLF=FALSE)

  res = lapply(residList[keep], function(x){
    idx = intersect(rownames(x), geneIds(gs))

    if( length(idx) < 10){
      x = NULL
    }else{
      x = x[idx,,drop=FALSE]
      # quantile normalize
      t(apply(x, 1, RankNorm))
    }
  })
  names(res) = names(residList[keep])
 
  df_cca = pairwise_cca( res, FALSE )

  if( ! is.null(df_cca) ){
    df = data.frame(Geneset = setName(gs), df_cca )
  }else{
    df = NULL
  }
  df
})
res = do.call(rbind, res)

saveRDS(res, "./res_cca.RDS")
saveRDS(residList, "./residList.RDS")

# scp sklar1:/sc/arion/work/hoffmg01/nps_ad/res_cca.RDS ~/Downloads
# scp sklar1:/sc/arion/work/hoffmg01/nps_ad/residList.RDS ~/Downloads

# res = readRDS("~/Downloads/res_cca.RDS")
# res = res[res$CT1 != "Unknown"& res$CT2 != "Unknown",]

# residList = readRDS("~/Downloads/residList.RDS")
```

```{r test,eval=FALSE}

library(mgcv)
library(data.table)

gsn = res$Geneset[1]


C = dfToMatrix(res[res$Geneset==gsn,], "CT1", "CT2", "cramer.V")

plotCorrMatrix(C^2)


with(res, plot(cramer.V, abs(cor1)))

with(res, plot(ri.x, ri.y))

ggplot(res, aes(ri.x, ri.y, color=log(n.x)-log(n.y))) + geom_point() + theme_classic()

i = with(res, which.max(abs(ri.x - ri.y)))


fit = gam(log(ri.x - ri.y+1e-4) ~ log(n.x)*log(n.y), data=res)

summary(fit)


res$chisq = scale(residuals(fit))^2

head(res[order(res$chisq, decreasing=TRUE),])

df2 = res[order(res$chisq, decreasing=TRUE),]

remove = c("Endothelial", "Oligodendrocyte", "Mural")

df3 = df2[!(df2$CT1 %in% remove) & !(df2$CT2 %in% remove),c(1,2,3,7, 10, 11, 14)]
df3 = data.table(df3)
df3[,id:=ifelse(CT1< CT2, paste(CT1, CT2), paste(CT2, CT1))]

# which cell pairs show hightest difference in redundancy
df3[,median(chisq), by="id"][order(V1),]





df = lapply( unique(res$Geneset), function(gsn){

  C = dfToMatrix(res[res$Geneset==gsn,], "CT1", "CT2", "cramer.V")
  dcmp = eigen(C)

  data.frame(Geneset = gsn, n = nrow(C), ev2 = dcmp$values[2] / sum(dcmp$values))
})
df = do.call(rbind, df)

ggplot(df, aes(n, ev2)) + geom_point() + geom_smooth(method="gam") + theme_classic()

df2 = df[df$n>7,]

fit = gam(ev2 ~ s(n,k=2), data=df2)

z = scale(residuals(fit))

head(df2[order(z, decreasing=TRUE),])




```



# vizualize results
gs = gsGO[[gsn]]

residSet = lapply(residList, function(X){
  ids = intersect(geneIds(gs), rownames(X))
  X[ids,]
})

CT1 = "Astrocyte"
CT2 = "Microglia"

CT1 = "Glutamatergic"
CT2 = "GABAergic"

X = t(residSet[[CT1]])
Y = t(residSet[[CT2]])
ids = intersect(rownames(X), rownames(Y))
fit = fastcca( X[ids,], Y[ids,] )

with(fit, plot(x.vars[,1], y.vars[,1]))
with(fit, cor(x.vars[,1], y.vars[,1]))



rs = lapply(fit.dl, function(fit){
  ids = intersect(geneIds(gs), rownames(fit))
  fit$Amean[ids]
})











```{r exit, cache=FALSE, eval=FALSE}
knitr::knit_exit()
```