chrom_loss_analysis.Rmd

---
title: "Analysis of NPS/AD"
subtitle: 'Loss of chrY'
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 /sc/arion/projects/CommonMind/hoffman/NPS-AD/work/nps_ad/
ml python git
git pull
R

system("git pull")
rmarkdown::render('chrom_loss_analysis.Rmd');




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


ml git
cd ~/build2/dreamlet
git pull
R CMD INSTALL .


--->

Examine distribution of chromosome loss events per cell 
  Is the probability of losing chrA related to loss of chrB after correcting for depth
  show overall loss rate per chromosome vs depth
  Is count per cell Poisson distributed?
  Use t-statistic as a proxy?


Donor level
  ashr within each cell type
  loss with age
  relate to phenotype



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

```{r load.packages, cache=FALSE}
library(tidyverse)
library(ggplot2)
library(cowplot)
library(lme4)
library(lmerTest)
library(SingleCellExperiment)
library(RNOmni)
```




```{r load.data}
# SingleCellExperiment 
outfolder = "/sc/arion/projects/CommonMind/hoffman/NPS-AD/work/nps_ad/freeze1_results/"
datafile = paste0(outfolder, "sceCombine_f1.RDS")
sceCombine = readRDS( datafile )

# collapse GABAergic neurons
sceCombine$celltype8 = sceCombine$celltype %>%
    as.character %>% 
    recode( 
    "GABAergic SST interneuron" = 'GABAergic neuron',
    "GABAergic VIP interneuron" = 'GABAergic neuron',
    "GABAergic PVALB interneuron" = 'GABAergic neuron') %>% 
    as.factor

# OUTRIDER results
outfolder = '/sc/arion/projects/CommonMind/hoffman/NPS-AD/work/nps_ad/chrom_loss_results/'
file = paste0(outfolder, "outliers.RDS")
res_combine = readRDS(file)

# adjust p-value due to merging
res_combine$padjust = p.adjust(res_combine$pValue, "fdr")
res_combine$aberrant = res_combine$padjust < 0.05
```

```{r analysis}
# get unique information for each Channel
df_uniq = with(colData(sceCombine), unique(data.frame(Channel, SubID, dx, Sex, Age)))

# allow grouping by Channel, geneID, and celltype8
# use data.table for faster merge
df_data = data.frame(sampleID = colnames(sceCombine), colData(sceCombine))

# total counts per chromosome
df_chr_counts = res_combine %>%
                group_by(geneID) %>%
                summarize(totalCounts = sum(rawcounts))

# combine outlier and sample data
df_res = left_join(res_combine, df_data, by="sampleID") %>%
  as_tibble %>% 
  mutate(Channel = gsub('-.*$','', sampleID),
         se = l2fc / zScore)  

# Transformation scores 
df_res$zScore.loss = -pmin(df_res$zScore, 0)
df_res$pValue.loss = pnorm(df_res$zScore, lower.tail=TRUE)
df_res$zScore.gain = pmax(df_res$zScore, 0)
df_res$pValue.gain = pnorm(df_res$zScore, lower.tail=FALSE)

# group and summarize
df_counts_aggr = df_res %>% 
  group_by(Channel, geneID, celltype8) %>% 
  summarize(nAbb = sum(aberrant),
            nObs = length(aberrant),
            nAbb.up = sum(aberrant*(l2fc>1)),
            nAbb.down = sum(aberrant*(l2fc<1)),
            abber.rate.up = sum(aberrant*(l2fc>1)) / length(aberrant),
            abber.rate.down = sum(aberrant*(l2fc<1)) / length(aberrant) ) %>% 
  left_join(df_uniq, by="Channel") %>%
  left_join(df_chr_counts, by="geneID")
```


# Correlation in gain/loss scores across chromosomes
```{r dependence, fig.height=20, fig.width=12}
library(broom)

chrs = c(1:22, "MT", "X", "Y")
chrs = factor(chrs, chrs)

grid = expand.grid( chrA = chrs, 
                    chrB = chrs)

df_dependence = lapply( 1:nrow(grid), function(i){
  message(i)

  chrA = grid$chrA[i]
  chrB = grid$chrB[i]

  a_join = df_res %>%
    select(geneID, sampleID, zScore, zScore.loss, zScore.gain) %>% 
    filter(., geneID == chrA)

  b_join = df_res %>%
    select(geneID, sampleID, zScore, zScore.loss, zScore.gain) %>% 
    filter(., geneID == chrB)

  ab = left_join( a_join, b_join, by='sampleID') 

  res1 = ab %>%
    group_modify(~ tidy(cor.test(.x$zScore.loss.x, .x$zScore.loss.y, 
      method="spearman"))) %>%
    mutate(chrA = chrA, chrB = chrB, Direction="loss")

  res2 = ab %>%
    group_modify(~ tidy(cor.test(.x$zScore.gain.x, .x$zScore.gain.y, 
      method="spearman"))) %>%
    mutate(chrA = chrA, chrB = chrB, Direction="gain")

  res3 = ab %>%
    group_modify(~ tidy(cor.test(.x$zScore.gain.x, .x$zScore.loss.y, 
      method="spearman"))) %>%
    mutate(chrA = chrA, chrB = chrB, Direction="mixed")

    rbind(res1, res2, res3)
})
df_dependence = do.call(rbind, df_dependence)
```

```{r show.dep, height=15}
library(Matrix)
library(corrplot)

col = rev(COL2('RdBu', 200))

C.list =lapply( c("loss", "gain", "mixed"), function(x){

  C = df_dependence %>%
      filter(Direction == x)  %>%
      with(sparseMatrix(i=as.numeric(chrA), j=as.numeric(chrB), x=estimate))
  rownames(C) = levels(df_dependence$chrA)
  colnames(C) = levels(df_dependence$chrA)
  as.matrix(C)
  })
names(C.list) = c("loss", "gain", "mixed")

par(mfrow=c(3,1))
corrplot(C.list[['loss']], col=col, tl.col="black", title="zScore.loss")
corrplot(C.list[['gain']], col=col, tl.col="black", title="zScore.gain")
corrplot(C.list[['mixed']], col=col, tl.col="black", title="zScore.mixed")
```

# Negative binomial model of aberrent counts
```{r aberrant.nb}
# only examine AD vs control
grid = expand.grid(geneID = sort(unique(df_counts_aggr$geneID)),
            celltype8 = sort(unique(df_counts_aggr$celltype8)))

coefs = c('scale(Age)','dxAD')

# test outlier rate
###################

res_lm = lapply(1:nrow(grid), function(i){
  message("\r", i, "   ")

  idx = with(df_counts_aggr,  (dx %in% c("AD", "Control")) & (geneID==grid$geneID[i]) & (celltype8==grid$celltype8[i]))
  data = droplevels(df_counts_aggr[idx,])
  data$dx = factor(data$dx, c("Control", "AD"))

  if( sum(data$nAbb> 0) < 100){
    df = NULL
  }else{

    df1 = df2 = NULL

    if( sum(data$nAbb.down > 1) > 10){
      fit1 = glmer.nb(nAbb.down ~ offset(log(nObs)) + (1|SubID) + dx + scale(Age), data)
      df1 = data.frame( geneID = grid$geneID[i],
                  celltype8 = grid$celltype8[i],
                  Direction="loss",
                  coef = coefs,
                  coef(summary(fit1))[coefs,,drop=FALSE])
    }

    if( sum(data$nAbb.up > 1) > 10){
      fit2 = glmer.nb(nAbb.up ~ offset(log(nObs)) + (1|SubID) + dx + scale(Age), data)

      
      df2 = data.frame( geneID = grid$geneID[i],
                  celltype8 = grid$celltype8[i],
                  Direction="gain",
                  coef = coefs,
                  coef(summary(fit2))[coefs,,drop=FALSE])
    }

    df = rbind(df1, df2)
  }
  df
  })
res_lm = do.call(rbind, res_lm)

colnames(res_lm)[colnames(res_lm) == "Pr...z.."] = "pvalue"
res_lm$FDR = p.adjust(res_lm$pvalue, "fdr")
res_lm = res_lm[order(res_lm$FDR),]

res_lm$geneID = factor(res_lm$geneID, c(1:22, "X", "Y", "MT"))
```

```{r plot1}
ratio = with(res_lm, nlevels(celltype8) / nlevels(geneID))

zmax = log10(max(abs(res_lm[res_lm$coef=="dxAD",]$z.value)))

ggplot(res_lm[res_lm$coef=="dxAD",], aes(geneID, celltype8, fill=sign(z.value)*log10(abs(z.value)))) + geom_tile() + theme_classic() + theme(aspect.ratio=ratio, plot.title = element_text(hjust = 0.5)) + scale_fill_gradient2(name=expression(sign~log[10]~z), low = "blue", mid="white", high="red", limits=c(-zmax, zmax))  + facet_wrap(~ Direction, ncol=1) + ggtitle("dxAD")

ggplot(res_lm[res_lm$coef=="scale(Age)",], aes(geneID, celltype8, fill=sign(z.value)*log10(abs(z.value)))) + geom_tile() + theme_classic() + theme(aspect.ratio=ratio, plot.title = element_text(hjust = 0.5)) + scale_fill_gradient2(name=expression(sign~log[10]~z), low = "blue", mid="white", high="red", limits=c(-zmax, zmax))  + facet_wrap(~ Direction, ncol=1) + ggtitle("Age")
```



chr21 trisomy is Down syndrome
 APP gene is on chr21

# Using z-score directly
```{r z.score}
grid = expand.grid(geneID = sort(unique(df_res$geneID)),
            celltype8 = sort(unique(df_res$celltype8)))

coefs = c('dxAD', 'scale(Age)')

res_lmz = lapply(1:nrow(grid), function(i){
  message("\r", i, "   ")

  idx = with(df_res, (dx %in% c("AD", "Control")) & (geneID==grid$geneID[i]) & (celltype8==grid$celltype8[i]))

  data = droplevels(df_res[idx,])
  data$dx = factor(data$dx, c("Control", "AD"))

  # test chromsome loss using pmin
  fit1 = lmer( RankNorm(zScore.loss) ~ scale(log10(rawcounts+1)) + (1|batch) + (1|Channel) + (1|SubID) + dx + scale(Age), data = data)
  
  # test chromsome gain using pmax
  fit2 = lmer( RankNorm(zScore.gain) ~ scale(log10(rawcounts+1)) + (1|batch) + (1|Channel) + (1|SubID) + dx + scale(Age), data = data)

  df1 = data.frame(chrom = grid$geneID[i], celltype8 = grid$celltype8[i],
                    Direction="loss", coef = coefs,
                    coef(summary(fit1))[coefs,])

  df2 = data.frame(chrom = grid$geneID[i], celltype8 = grid$celltype8[i],
                    Direction="gain", coef = coefs,
                    coef(summary(fit2))[coefs,])

  rbind(df1, df2)
})
res_lmz = do.call(rbind, res_lmz)

colnames(res_lmz)[colnames(res_lmz) == "Pr...t.."] = "pvalue"
res_lmz$FDR = p.adjust(res_lmz$pvalue, "fdr")
res_lmz = res_lmz[order(res_lmz$FDR),]

res_lmz$geneID = factor(res_lmz$chrom, c(1:22, "X", "Y", "MT"))
```

```{r plot}
ratio = with(res_lmz, nlevels(celltype8) / nlevels(geneID))
zlim = max(abs(res_lmz$t.value))

ggplot(res_lmz[res_lmz$coef=="dxAD",], aes(geneID, celltype8, fill=t.value)) + geom_tile() + theme_classic() + theme(aspect.ratio=ratio, plot.title = element_text(hjust = 0.5)) + scale_fill_gradient2(low = "blue", mid="white", high="red", limits=c(-zlim, zlim)) + facet_wrap(~ Direction, ncol=1) + ggtitle("dxAd")

ggplot(res_lmz[res_lmz$coef=="scale(Age)",], aes(geneID, celltype8, fill=t.value)) + geom_tile() + theme_classic() + theme(aspect.ratio=ratio, plot.title = element_text(hjust = 0.5)) + scale_fill_gradient2(low = "blue", mid="white", high="red", limits=c(-zlim, zlim)) + facet_wrap(~ Direction, ncol=1) + ggtitle("Age")
```



```{r save}
outfile = paste0("/sc/arion/projects/CommonMind/hoffman/NPS-AD/work/nps_ad/sce_mosaic_aneuploidy_AD", gsub("-", "_", Sys.Date()), ".tsv")
write.table(df_res, file=outfile, row.names=FALSE, quote=FALSE, sep="\t")
R.utils::gzip(outfile, overwrite=TRUE)
```


# Predict AD
Association with AD and number of cells with chromosome *loss*.

```{r predict.AD}
# df_res = fread(outfile)
grid = expand.grid(geneID = sort(unique(df_res$geneID)),
            celltype8 = sort(unique(df_res$celltype8)))

res_pred = lapply(1:nrow(grid), function(i){
  message("\r", i, "   ")

  idx = with(df_res, (dx %in% c("AD", "Control")) & (geneID==grid$geneID[i]) & (celltype8==grid$celltype8[i]))

  data = droplevels(df_res[idx,])
  data$dx = factor(data$dx, c("Control", "AD"))

  form1 = dx ~ zScore.loss
  if(grid$geneID[i] == "Y"){ 
    form2 = dx ~ RankNorm(ad_pred_cells) + Age 
  }else{
    form2 = dx ~ RankNorm(ad_pred_cells) + Age + Sex
  }

  # predict at cell level
  fit = glm(form1, data, family="binomial")
  data$ad_pred = predict(fit, type="response")

  # collapse by subject
  data2 = data %>%
            group_by(SubID) %>%
            summarize(nAbb = sum(aberrant),
            nObs = length(aberrant),
            nAbb.up = sum(aberrant*(l2fc>1)),
            nAbb.down = sum(aberrant*(l2fc<1)),
            abber.rate.up = sum(aberrant*(l2fc>1)) / length(aberrant),
            abber.rate.down = sum(aberrant*(l2fc<1)) / length(aberrant),
            Age = Age[1],
            Sex = Sex[1],
            Institution = Institution[1],
            dx = dx[1],
            # ad_pred_cells = sum(ad_pred > 0.7) / length(ad_pred),
            ad_pred_cells = median(ad_pred) )

  data2 %>% group_by(dx) %>% summarize(a = mean(ad_pred_cells))

  # predict AD at Subject-level using number of AD cells
  # + abber.rate.down + abber.rate.up
  fit = glm(form2, data2, family="binomial")

  # training vs testing
  # testing = c(sample(which(data2$dx=="AD"), 50), 
  #              sample(which(data2$dx=="Control"), 50))

  # fit = glm(form2, data2[-testing,], family="binomial")
  # y_pred = predict(fit, data2[testing,])

  # lm(y_pred ~ data2$dx[testing])

  data.frame(chrom = grid$geneID[i], celltype8 = grid$celltype8[i],
    t(coef(summary(fit))['RankNorm(ad_pred_cells)',]))
})
res_pred = do.call(rbind, res_pred)
```

```{r plot.res_pred, cache=FALSE}
ratio = with(res_pred, nlevels(celltype8) / nlevels(chrom))
zlim = max(abs(res_pred$z.value))

res_pred$chrom = factor(res_pred$chrom, c(1:22, "X", "Y", "MT"))
colnames(res_pred)[colnames(res_pred) == "Pr...z.."] = "pvalue"
res_pred$fdr = p.adjust(res_pred$pvalue, "fdr")

ggplot(res_pred, aes(chrom, celltype8, fill=z.value)) + geom_tile() + theme_classic() + theme(aspect.ratio=ratio) + scale_fill_gradient2(low = "blue", mid="white", high="red", limits=c(-zlim, zlim)) + geom_text(aes(label=ifelse(fdr < 0.05, '*', '')))
```

```{r across.chrom}
library(PRROC)
library(glmnet)

data = df_res %>%
        filter(dx %in% c("AD", "Control", "SCZ", "BP", "NA")) %>%
        select(geneID, celltype8, SubID, Channel, zScore, zScore.loss, zScore.gain, Age, Sex, dx)  %>%
        # mutate(dx = factor(dx, c("Control", "AD"))) %>%
        group_by(SubID, geneID, celltype8) %>%
        summarize( zScore = mean(zScore),
                   zScore.loss = mean(zScore.loss),
                   zScore.gain = mean(zScore.gain),
                   SubID = SubID[1],
                   Sex = Sex[1],
                   Age = Age[1],
                   dx = dx[1])

df = data %>%
      pivot_wider( names_from = c(celltype8, geneID), values_from=c(zScore, zScore.loss, zScore.gain)) %>%
      droplevels


keep = grep("loss|gain", colnames(df))

X = as.data.frame(df[,c(2,3,keep)])
X$Sex = as.numeric(X$Sex)
X[is.na(X)] = 0
X = scale(as.matrix(X))

set.seed(1)
training = sample.int(nrow(X), 0.7*nrow(X))

# Training
keys = c("Sex", "Age", "Sex + Age", "Sex + Age + MA")
df_pred = lapply( keys, function(x){

  if( x == "Sex" | x == "Age"){
    idx = c(100, which(colnames(X) == x))
  }else if(x == "Sex + Age"){
    idx = which(colnames(X) %in% c("Sex", "Age"))
  }else if( x == 'Sex + Age + MA'){
    idx = 1:ncol(X)
  }

  if( x == 'Sex + Age + MA' ){
    pf = rep(1, length(idx))
    pf[1:2] = 0
    fitcv = cv.glmnet(X[training,idx], df$dx[training], family="multinomial", type.multinomial = "grouped", penalty.factor=pf)

    png("test.png")
    plot(fitcv)
    dev.off()

  }else{    
    fitcv = glmnet(X[training,idx,drop=FALSE], df$dx[training], family="multinomial", type.multinomial = "grouped", lambda=0)
  }

  # same coef image
  if( x == "Sex + Age + MA"){
    beta = coef(fitcv, s = "lambda.1se")

    df_beta = lapply( names(beta), function(x){
      value = as.matrix(beta[[x]][-c(1:3),])
      data.frame(coef=rownames(value), key=x, value = c( value))
      })
    df_beta = do.call(rbind, df_beta)

    df_beta = cbind(df_beta, df_beta$coef %>% 
                strsplit("_") %>% 
                do.call(rbind, .))
    colnames(df_beta)[-c(1:3)] = c("score", "celltype8", "chrom")

    df_beta$chrom = factor(df_beta$chrom, c(1:22, "X", "Y", "MT")) 

    ratio = with(df_beta, length(unique(celltype8)) / nlevels(chrom))

    fig = ggplot(df_beta, aes(chrom, celltype8, fill=value)) + 
              geom_tile() +
              theme_classic() +
              theme(aspect.ratio=ratio) +
              facet_wrap(~ key + score, ncol=2) +
              scale_fill_gradient2(low="blue", mid="white", high="red")
    ggsave("test.png", fig, width=16, height=12)
  }

  # Performance on testing set
  y_pred = predict(fitcv, newx = X[-training,idx], s = "lambda.1se", type="response")
  y_pred = as.data.frame(y_pred[,,])

  # plotCorrMatrix(cor(y_pred))

  # evaluate performance
  df_pred = lapply(levels(df$dx), function(status){
    message(status)
    y_pred_case = y_pred[df$dx[-training]==status, status]
    y_pred_not = y_pred[df$dx[-training]=='Control', status]

    y_pred_case = y_pred_case[!is.na(y_pred_case)]
    y_pred_not = y_pred_not[!is.na(y_pred_not)]

    pr <- pr.curve(scores.class0 = y_pred_case, 
                   scores.class1 = y_pred_not, 
                   rand.compute=TRUE, 
                   curve=TRUE)

    data.frame(Class = status,
                Set = x,
               AUPR = pr$auc.integral,
               AUPR.rand = pr$rand$auc.integral,
               n.class = length(y_pred_case))
  })
  do.call(rbind, df_pred)
})
df_pred = do.call(rbind, df_pred)

df_pred$Set = factor(df_pred$Set, keys)

ggplot(df_pred, aes(Set, AUPR, fill=Set)) + 
  facet_wrap(~ Class) +
  geom_bar(stat="identity") +
  theme_classic() +
  theme(aspect.ratio=1, axis.text.x = element_text(angle = 75, vjust = 1, hjust=1), legend.position="none") +
  scale_fill_brewer(palette="Set1") +
  geom_hline(aes(yintercept=AUPR.rand), linetype="dashed", color="grey50") +
  scale_y_continuous(limits=c(0, 1), expand=c(0,0)) +
  geom_text(aes(y=AUPR+.01, label=format(AUPR, digits=3)), size=2, vjust=0)
```

```{r temp, eval=FALSE}

# status = c("Control", "SCZ")
# df_sub = df[-training,]
# df_sub$y_pred = y_pred[,status[2],1]
# idx = df_sub$dx %in% status
# df_sub = droplevels(df_sub[idx,])
# df_sub$dx = factor(df_sub$dx, status)

# fit = glm(dx ~ Sex + Age, df_sub, family="binomial")
# calcVarPart(fit)

# fit = glm(dx ~ y_pred + Sex + Age, df_sub, family="binomial")
# summary(fit)
# calcVarPart(fit)
```

```{r test, eval=FALSE}
# Evaluate AUPR
################

# train 
fit.train = glm(dx ~  Age, df[training,], family="binomial")

# test
X_new = model.matrix(~ Sex + Age, df[-training,])
y_pred = predict(fit.train, newx=X_new)
y_pred_ad = y_pred[df$dx[-training]=="AD"]
y_pred_ctrl = y_pred[df$dx[-training]=="Control"]

y_pred_ad = y_pred_ad[!is.na(y_pred_ad)]
y_pred_ctrl = y_pred_ctrl[!is.na(y_pred_ctrl)]

pr <- pr.curve(scores.class0 = y_pred_ad, 
               scores.class1 = y_pred_ctrl, 
               rand.compute=TRUE, 
               curve=TRUE)

pr


# lasso
#####

idx = 1:ncol(X)
# train 
fit.train = cv.glmnet(X[training,idx], df$dx[training], family="binomial", dfmax=5 )

# test
y_pred = predict(fit.train, newx = X[-training,idx], s = "lambda.1se")
y_pred_ad = y_pred[df$dx[-training]=="AD"]
y_pred_ctrl = y_pred[df$dx[-training]=="Control"]

y_pred_ad = y_pred_ad[!is.na(y_pred_ad)]
y_pred_ctrl = y_pred_ctrl[!is.na(y_pred_ctrl)]

pr <- pr.curve(scores.class0 = y_pred_ad, 
               scores.class1 = y_pred_ctrl, 
               rand.compute=TRUE, 
               curve=TRUE)

pr$auc.integral
```


# Fraction of cells
```{r cell.frac}
library(qvalue)
get_pi1 = function(p){

  res <- tryCatch( pi0est(p), error = function(e) NULL)

  pi1 = NA
  if( ! is.null(res)){
    pi1 = 1 - res$pi0
  }
  pi1
}

 data2 = df_res %>%
            filter(celltype == 'Glutamatergic neuron', dx %in% c("AD", "Control")) %>%
            group_by(SubID, geneID, celltype) %>%
            summarize(nAbb = sum(aberrant),
              nObs = length(aberrant),
              nAbb.up = sum(aberrant*(l2fc>1)),
              nAbb.down = sum(aberrant*(l2fc<1)),
              abber.rate.up = sum(aberrant*(l2fc>1)) / length(aberrant),
              abber.rate.down = sum(aberrant*(l2fc<1)) / length(aberrant),
              mu.z.loss = mean(zScore.loss),
              mu.z.gain = mean(zScore.gain),
              pi.loss = get_pi1(2*pnorm(zScore.loss, lower.tail=FALSE)),
              pi.gain = get_pi1(2*pnorm(zScore.gain, lower.tail=FALSE)),
              Age = Age[1],
              Sex = Sex[1],
              Institution = Institution[1],
              dx = dx[1]) %>% 
           mutate( geneID = factor(geneID, c(1:22, "MT", "X", "Y")))
               

```            

```{r mu.z.loss, fig.height=12}
idx = with(data2, (dx %in% c("AD", "Control", "SCZ") ))

ggplot(data2[idx,], aes(dx, mu.z.loss, fill=Sex)) + geom_violin(color=NA) + geom_boxplot(width=0.1, position=position_dodge(width=0.9), outlier.size=1) + facet_wrap(~ geneID, ncol=4) + theme_classic() + theme(aspect.ratio=1) + scale_fill_manual(values=c("red", "blue"))
```

```{r mu.z.gain, fig.height=12}
ggplot(data2[idx,], aes(dx, mu.z.gain, fill=Sex)) + geom_violin(color=NA) + geom_boxplot(width=0.1, position=position_dodge(width=0.9), outlier.size=1) + facet_wrap(~ geneID, ncol=4) + theme_classic() + theme(aspect.ratio=1) + scale_fill_manual(values=c("red", "blue"))
```


```{r abber.rate.down, fig.height=12}
ggplot(data2[idx,], aes(dx, abber.rate.down, fill=Sex)) + geom_violin(color=NA) + geom_boxplot(width=0.1, position=position_dodge(width=0.9), outlier.size=1) + facet_wrap(~ geneID , ncol=4) + theme_classic() + theme(aspect.ratio=1) + scale_fill_manual(values=c("red", "blue"))  + scale_y_continuous(limits=c(0,1), expand=c(0,0))
```

```{r abber.rate.up, fig.height=12}
ggplot(data2[idx,], aes(dx, abber.rate.up, fill=Sex)) + geom_violin(color=NA) + geom_boxplot(width=0.1, position=position_dodge(width=0.9), outlier.size=1) + facet_wrap(~ geneID, ncol=4) + theme_classic() + theme(aspect.ratio=1) + scale_fill_manual(values=c("red", "blue"))  + scale_y_continuous(limits=c(0,1), expand=c(0,0))
```



```{r pi.loss, fig.height=12}
ggplot(data2[idx,], aes(dx, pi.loss, fill=Sex)) + geom_violin(color=NA) + geom_boxplot(width=0.1, position=position_dodge(width=0.9), outlier.size=1) + facet_wrap(~ geneID , ncol=4) + theme_classic() + theme(aspect.ratio=1) + scale_fill_manual(values=c("red", "blue"))  + scale_y_continuous(limits=c(0,1), expand=c(0,0))
```



```{r pi.gain, fig.height=12}
ggplot(data2[idx,], aes(dx, pi.gain, fill=Sex)) + geom_violin(color=NA) + geom_boxplot(width=0.1, position=position_dodge(width=0.9), outlier.size=1) + facet_wrap(~ geneID , ncol=4) + theme_classic() + theme(aspect.ratio=1) + scale_fill_manual(values=c("red", "blue")) + scale_y_continuous(limits=c(0,1), expand=c(0,0))
```




```{r lm}
data2$dx = factor(data2$dx, c("Control", "AD"))
fit = lm(mu.z.loss ~ dx + Age + Institution + Sex, data2 %>% filter(geneID == "21") )
summary(fit)

fit = glm(dx ~ mu.z.loss + Age + Institution + Sex, data2 %>% filter(geneID == "21"), family="binomial" )
summary(fit)

```


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













```{r pi1, eval=FALSE}
library(qvalue)
library(ashr)
library(lmerTest)
library(data.table)
library(tidyverse)

# saveRDS(res, "/sc/arion/scratch/hoffmg01/res.RDS")
# saveRDS(ods, "/sc/arion/scratch/hoffmg01/ods.RDS")
# res = readRDS("/sc/arion/scratch/hoffmg01/res.RDS")
# ods = readRDS("/sc/arion/scratch/hoffmg01/ods.RDS")

pi1 = function(beta, se){
  if( length(beta) < 10) return(NA)
  # fit = ash(beta, se)
  fit = ash(beta/se, rep(1, length(se)))
  1 - get_pi0(fit)
}

# get unique information for each Channel
df_uniq = with(colData(sceCombine), unique(data.frame(Channel, SubID, dx, Sex, Age)))

# allow grouping by Channel, geneID, and celltype
# use data.table for faster merge
df_data = data.table(data.frame(sampleID = colnames(sceCombine), colData(sceCombine)))
setkey(df_data, sampleID)
setkey(res, sampleID)

# total counts per chromosome
df_chr_counts = res %>%
                group_by(geneID) %>%
                summarize(totalCounts = sum(rawcounts))

df_res = left_join(res, df_data, by="sampleID") %>%
  as_tibble %>% 
  mutate(Channel = gsub('-.*$','', sampleID),
         se = l2fc / zScore)  


# Get outliers for aggregated across groups
df_counts = df_res %>% 
  group_by(Channel, geneID, celltype) %>% 
  summarize(nAbb = sum(aberrant),
            nObs = length(aberrant),
            nAbb.up = sum(aberrant*(l2fc>1)),
            nAbb.down = sum(aberrant*(l2fc>1)),
            abber.rate.up = sum(aberrant*(l2fc>1)) / length(aberrant),
            abber.rate.down = sum(aberrant*(l2fc<1)) / length(aberrant) )


df_counts = merge(df_counts, df_uniq, by="Channel")

library(lme4)
# fit = glmer.nb(nAbb ~ offset(log(nObs)) + (1|SubID) + (1|celltype), df_counts[df_counts$geneID==2,])

df_vp = lapply(1:22, function(x){
  fit = glmer(nAbb ~ offset(log(nObs)) + (1|SubID) + (1|celltype) , df_counts[df_counts$geneID==x,], family="poisson")

  data.frame(chrom = x, t(calcVarPart(fit)))
})
df_vp = do.call(rbind, df_vp)



df_counts_aggr = df_res %>% 
  group_by(Channel, geneID, celltype) %>% 
  summarize(nAbb = sum(aberrant),
            nObs = length(aberrant),
            nAbb.up = sum(aberrant*(l2fc>1)),
            nAbb.down = sum(aberrant*(l2fc>1)),
            abber.rate.up = sum(aberrant*(l2fc>1)) / length(aberrant),
            abber.rate.down = sum(aberrant*(l2fc<1)) / length(aberrant) ) %>% 
  left_join(df_uniq, by="Channel") %>%
  left_join(df_chr_counts, by="geneID")


fit = glmer(nAbb ~ offset(log(nObs)) + (1|Channel) + (1|SubID) + (1|geneID) + (1|celltype) + scale(log(totalCounts)) + dx, df_counts_aggr[df_counts_aggr$dx %in% c("AD", "Control"),], family="poisson")

calcVarPart(fit)
summary(fit)

df_sub = df_counts_aggr[df_counts_aggr$dx %in% c("AD", "Control"),]

grid = expand.grid(geneID = sort(unique(df_sub$geneID)),
            celltype = sort(unique(df_sub$celltype)))

# test outlier rate
###################
# note: outliers only called in Males
rm(res_lm)
res_lm = mclapply(1:nrow(grid), function(i){
  cat("\r", i, "   ")

  idx = with(df_sub, (geneID==grid$geneID[i]) & (celltype==grid$celltype[i]))
  data = droplevels(df_sub[idx,])
  data$dx = factor(data$dx, c("Control", "AD"))

  if( sum(data$nAbb> 0) < 10){
    df = NULL
  }else{

    fit = glmer(nAbb ~ offset(log(nObs)) + (1|SubID) + dx + scale(Age)*dx, data, family="poisson")

    coefs = c('scale(Age)','dxAD', 'dxAD:scale(Age)')
    df = data.frame( geneID = grid$geneID[i],
                celltype = grid$celltype[i],
                coef = coefs,
                coef(summary(fit))[coefs,,drop=FALSE])
  }
  df
  }, mc.cores=12)
res_lm = do.call(rbind, res_lm)

colnames(res_lm)[colnames(res_lm) == "Pr...z.."] = "pvalue"
res_lm$FDR = p.adjust(res_lm$pvalue, "fdr")
res_lm = res_lm[order(res_lm$FDR),]

res_lm$geneID = factor(res_lm$geneID, c(1:22, "X", "Y", "MT"))

ratio = with(res_lm, nlevels(celltype) / nlevels(geneID))

fig = ggplot(res_lm[res_lm$coef=="scale(Age)",], aes(geneID, celltype, fill=z.value)) + geom_tile() + theme_classic() + theme(aspect.ratio=ratio) + scale_fill_gradient2(low = "blue", mid="white", high="red")

ggsave("test.png", fig)



# Using z-score directly
#########################

chr21 trisomy is Down syndrome
APP gene is on chr21

df_sub = df_res[df_res$dx %in% c("AD", "Control"),]

grid = expand.grid(geneID = sort(unique(df_sub$geneID)),
            celltype = sort(unique(df_sub$celltype)))

res_lm = lapply(1:nrow(grid), function(i){
  cat("\r", i, "   ")

  idx = with(df_res, (geneID==grid$geneID[i]) & (celltype==grid$celltype[i]))

  data = droplevels(df_res[idx,])
  data$dx = factor(data$dx, c("Control", "AD"))

  # test chromsome loss using pmin
  fit1 = lmer( pmin(zScore, 0) ~ scale(log10(rawcounts+1)) + (1|batch) + (1|Channel) + (1|SubID) + dx + scale(Age), data = data)
  
  # test chromsome gain using pmax
  fit2 = lmer( pmax(zScore, 0) ~ scale(log10(rawcounts+1)) + (1|batch) + (1|Channel) + (1|SubID) + dx + scale(Age), data = data)

  coefs = c('dxAD', 'scale(Age)')
  df1 = data.frame(chrom = grid$geneID[i], celltype = grid$celltype[i],
                    Direction="loss", coef = coefs,
                    coef(summary(fit1))[coefs,])

  df2 = data.frame(chrom = grid$geneID[i], celltype = grid$celltype[i],
                    Direction="gain", coef = coefs,
                    coef(summary(fit2))[coefs,])

  rbind(df1, df2)
})
res_lm = do.call(rbind, res_lm)

colnames(res_lm)[colnames(res_lm) == "Pr...t.."] = "pvalue"
res_lm$FDR = p.adjust(res_lm$pvalue, "fdr")
res_lm = res_lm[order(res_lm$FDR),]

res_lm$geneID = factor(res_lm$chrom, c(1:22, "X", "Y", "MT"))

ratio = with(res_lm, nlevels(celltype) / nlevels(geneID))

fig = ggplot(res_lm[res_lm$coef=="scale(Age)",], aes(geneID, celltype, fill=z.value)) + geom_tile() + theme_classic() + theme(aspect.ratio=ratio) + scale_fill_gradient2(low = "blue", mid="white", high="red")

ggsave("test.png", fig)












# Using z-score directly
# lm() uses too much memory
# test AD-control directly
idx = with(df_res, (geneID=="21") & (dx %in% c('AD', 'Control')) & (celltype=='Glutamatergic neuron'))
df_sub = df_res[idx,]
df_sub = droplevels(df_sub)
df_sub$dx = factor(df_sub$dx, c("Control", "AD"))
fit = lmer( pmax(zScore, 0) ~ scale(rawcounts) + (1|Channel) + (1|SubID) + scale(Age) + dx, data = df_sub)

summary(fit)




df_ad = lapply( sort(unique(df_res$geneID)), function(chrom){
  df_sub = df_res[(df_res$geneID=="Y") &(df_res$dx %in% c('AD', 'Control')),]
  df_sub = droplevels(df_sub)

  # test chromsome loss using pmin
  fit1 = lmer( pmin(zScore, 0) ~ scale(log10(rawcounts+1)) + (1|celltype) + (1|Channel) + (1|SubID) + as.numeric(dx), data = df_sub)
  
  # test chromsome gain using pmax
  fit2 = lmer( pmax(zScore, 0) ~ scale(log10(rawcounts+1)) + (1|celltype) + (1|Channel) + (1|SubID) + as.numeric(dx), data = df_sub)

  rbind(data.frame(chrom = chrom, Direction="loss", t(coef(summary(fit1))['as.numeric(dx)',])),
        data.frame(chrom = chrom, Direction="gain", t(coef(summary(fit2))['as.numeric(dx)',])))
})
df_ad = do.call(rbind, df_ad)





# pi1 = 1 - qvalue(pValue)$pi0
df_pi1 = df_res %>% 
  group_by(Channel, geneID, celltype) %>% 
  summarize( pi1.ash = pi1(l2fc, se))


idx = with(df_res, which(Channel == "H1004_1" & geneID == "10"))

with(df_res[idx,], pi1(l2fc, se))
```



# need to exclude replicates




# or should I use a negative binomial model directly and the produce residauls.
```{r results}
head(results(ods))

nrow(results(ods,all=TRUE))

table(results(ods,all=TRUE)$padjust < 0.05)
```


fit = MGLMreg(cbind(y1, y2, y3, y4) ~ 1, data=ydata, dist='DM')

fit@coefficients / fit@SE