13_cytokineResistancePhospho.Rmd

---
title: "Cytokine Resistance Phosphodata"
author: "Sara Gosline"
date: "4/9/2021"
output: html_document
---


```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE)
##process cytokine data
library(amlresistancenetworks)
library(dplyr)
library(ggplot2)
```

## Phospho data kinase enrichment

For this analysis will load in the phosphoproteomic data collected in the Molm13 cell lines. We use the stoichiometry corrected proteomic data to focus explicitly on the changes in kinase activity.  

```{r data load,echo=FALSE,warning=FALSE}

###load all the data

phosData<-querySynapseTable('syn22986341')%>%subset(!is.nan(LogRatio))%>%
  mutate(Gene=unlist(Gene))%>%
  mutate(site=unlist(site))


clinvars<-phosData%>%
  dplyr::select(Sample='sample',CellType,TimePoint,Treatment,Batch)%>%
  distinct()

##what are we doing again?
earlySummary<-phosData%>%
    subset(Batch=='Experiment 1')%>%
  dplyr::select(sample,CellType,TimePoint,Treatment,Batch)%>%
  distinct()%>%
  mutate(conditionName=stringr::str_c(CellType,TimePoint,Treatment,sep='_'))

lateSummary<-phosData%>%
  subset(TimePoint==0)%>%
    subset(Batch=='Experiment 2')%>%
  dplyr::select(sample,CellType,TimePoint,Treatment,Batch)%>%
  distinct()%>%
  mutate(conditionName=stringr::str_c(CellType,TimePoint,Treatment,sep='_'))


earlyPhosMat<-phosData%>%
  subset(Batch=='Experiment 1')%>%
  dplyr::select(sample,site,LogRatio)%>%
  tidyr::pivot_wider(values_from=LogRatio,names_from=sample,
                     values_fn=list(LogRatio=mean,na.rm=T))%>%#updated 9/6,values_fill=list(LogRatio=0.0))%>%
  tibble::remove_rownames()%>%
  tibble::column_to_rownames('site')



latePhosMat<-phosData%>%
  subset(TimePoint==0)%>%
  subset(Batch=='Experiment 2')%>%
  dplyr::select(sample,site,LogRatio)%>%
  tidyr::pivot_wider(values_from=LogRatio,names_from=sample,
                     values_fn=list(LogRatio=mean,na.rm=T))%>%#updated 9/6,values_fill=list(LogRatio=0.0))%>%
  tibble::remove_rownames()%>%
  tibble::column_to_rownames('site')
```

## Now we can define the processing functions needed


We can plot the data and the kinase activity as needed in heatmap or PCA form.
```{r processing functions,echo=TRUE,warning=FALSE}
##
#' @param dat.table
# 2- Create PCA plots of global and phospho for the late molm13 only, and the late molm 13 + parental

#' @param dat.table
plotAllData <- function(dat.table, prefix='',expand=0.01, alpha=0.1, ...) {
  library(ggfortify)
  met <- dat.table %>%
    dplyr::select(sample, CellType, TimePoint, Treatment) %>%
    distinct()
    
  mat <- dat.table %>% dplyr::select(Gene,LogRatio,sample) %>%
    distinct() %>%
    mutate(LogRatio=as.numeric(LogRatio)) %>%
    tidyr::pivot_wider(names_from='sample', values_from='LogRatio',
                       values_fn=list(LogRatio=function(x) mean(x,na.rm=T)),
                       values_fill=list(LogRatio=0)) %>%
    tibble::remove_rownames() %>%
    tibble::column_to_rownames('Gene')
  mat <- mat[complete.cases(mat),]
  x <- prcomp(t(mat))$x
  x <- as.data.frame(x)
  x$sample <- rownames(x)
  ggdata <- inner_join(x, met, by="sample")
  
  library(ggforce)
  ggplot(ggdata, aes_string(x="PC1", y="PC2", ...)) +
    geom_point(size=2.5) +
    ggforce::geom_mark_ellipse(color = NA, alpha=alpha, expand=expand)

}


##plot kinase activity
plotKinDat<-function(kindat,sig.kin=NULL,prefix='all'){
  library(pheatmap)
  
  ##create matrix of kinase scores
  if(!is.null(sig.kin)){
    kindat<-subset(kindat,Kinase%in%sig.kin)#$Kinase.Gene)
    #kinmat<-sig.kin%>%mutate(score='Yes')%>%
    #  tidyr::pivot_wider(names_from=Condition,values_from=score,values_fill=list(score='No'))%>%
    #  tibble::column_to_rownames('Kinase.Gene')
  
  }
  
  mat <-kindat%>%
    ungroup()%>%
    tidyr::pivot_wider(-c(meanNKINscore,numSubstr),
                                              values_from=meanLFC,
                                                names_from=Sample,
                                                values_fn=list(meanLFC=mean), values_fill=0.0)%>%
    tibble::remove_rownames()%>%
    tibble::column_to_rownames('Kinase')
  
  #if(!is.null(sig.kin))
  #  kinAts=kinmat
  #else
    kinAts<-kindat%>%
      ungroup()%>%
      dplyr::select(Kinase,numSubstr)%>%
      distinct()%>%
      group_by(Kinase)%>%
      summarize(substrates=mean(numSubstr))%>%
      tibble::remove_rownames()%>%
      tibble::column_to_rownames('Kinase')
  
  sampAts<-phosData%>%
    dplyr::select(sample,TimePoint,Treatment)%>%
    distinct()%>%
    tibble::remove_rownames()%>%
    tibble::column_to_rownames('sample')
  
  
  sampAts$TimePoint=as.factor(sampAts$TimePoint)
  
  pal<-nationalparkcolors::park_palette('GrandTeton')
  sampCols= lapply(names(sampAts),function(x) {
    p<-pal[1:length(unique(sampAts[[x]]))]
    names(p)<-unique(sampAts[[x]])
    return(p)})

  names(sampCols)=names(sampAts)
  #print(sampCols)
  #vars=names(sort(apply(mat,1,var,na.rm=T),decreasing=T))
  pheatmap(mat,cellwidth = 8,cellheight=8,clustering_distance_cols = 'correlation',
          clustering_distance_rows = 'correlation', 
          annotation_colors=sampCols,annotation_col=sampAts) 

  pheatmap(mat,cellwidth = 8,cellheight=8,clustering_distance_cols = 'correlation',
          clustering_distance_rows = 'correlation',
          annotation_col=sampAts,annotation_colors = sampCols,
          file=paste0(prefix,'cytokineKinaseHeatmap.pdf'),height=14,width=8) 

  
  pheatmap(mat,cellwidth = 8,cellheight=8,cluster_cols=FALSE,
          clustering_distance_rows = 'correlation',
          annotation_col=sampAts,annotation_colors = sampCols,
          file=paste0(prefix,'cytokineKinaseHeatmapUnclustered.pdf'),height=14,width=8) 
  }


####Show kinase activity
kindat<-phosData%>%
  subset(Batch=='Experiment 1')%>%
  dplyr::rename(Sample='sample', LogFoldChange='LogRatio')%>%
  mapPhosphoToKinase()


##plot PCA
phosData%>%
  subset(Batch=='Experiment 1')%>%
  subset(CellType=='MOLM-13')%>%
  plotAllData(.,colour='Treatment',fill='Treatment')+
  ggtitle("Phosphoproteomics Experiment 1")

```

## KSEA enrichment analysis

Next we define functions to carry out various KSEA analyses. 
```{r ksea, echo=FALSE,message=FALSE}
#' plot all the KSEA 
#' @param condList
#' @return data frame
doAllKSEAplots<-function(condList,pdat=phosData){
  
  gene.to.site<-dplyr::select(pdat,Gene,site,Peptide)%>%distinct()%>%
    dplyr::mutate(residue=stringr::str_replace(site,paste0(Gene,'-'),''))%>%
    dplyr::mutate(residue=stringr::str_replace_all(residue,"([STY])", ";\\1"))%>%
    dplyr::mutate(residue=stringr::str_replace(residue,"^;", ""))%>%
    dplyr::mutate(residue=stringr::str_replace_all(residue,"([sty])", ""))
  
  full.df<-purrr::map_df(names(condList),.f=function(clName){ 
    print(clName)
    condList[[clName]]%>%
      tibble::rownames_to_column('site')%>%
      left_join(gene.to.site)%>%
      dplyr::select(Gene,Peptide,residue,value='logFC',p_adj='adj.P.Val')%>%
      amlresistancenetworks::computeKSEA(.,prefix=clName,0.05,linkedSubs=1,suffix='pdf')%>%
      mutate(Condition=clName)%>%
      as.data.frame()
  })
  return(full.df)
  
}


#' build networks from data frame
#' @param data.res
#' @param gene.col
#' @param weight.col
#' @param condition.col
#' @return network list?
runNetworksFromDF<-function(data,gene.col='Kinase.Gene',
                              weight.col='aveSubstrateLog2FC',
                              pvalue.col='P.Value',
                              condition.col='Condition',extra.col=c('Substrate.Gene','Source','log2FC'),
                              signif=0.01,beta=5){
  conditions=unique(dplyr::select(data,condition.col))
  for(con in conditions$condition){
    dat = data%>%
      dplyr::select(cond=condition.col,value=weight.col,Gene=gene.col,p.value=pvalue.col)%>%
      subset(p.value<signif)%>%
      subset(cond==con)
    vals<-dat$value
    names(vals)<- stringr::str_replace_all(dat$Gene,"([sty])", "")
    print(length(vals))
    #print(con)
    #print(vals)
    amlresistancenetworks::computePhosphoNetwork(phos.vals=vals,beta=beta,fname=paste0('phosphoNetworkfor',con,'_beta',beta))
  }
   
}
```

### Single agent time course comparisons

We want to identify changes in each condition at 5 min and 60 min compare to zero. From here we can see the changes in predicted kinase activity upon treatment. Specifically these comparisons include:
1. 5 min Trametinib vs. control
2. 60 min Trametinib vs. control
3. 5 min MCP-1 vs. control
4. 60 min MCP-1 vs. control
5. 5 min MCP-1+Trametinib vs. control
6. 60 min MCP-1+Trametinib vs. control

```{r early time course,warning=FALSE,echo=FALSE}

phosMat <- earlyPhosMat
summary <- earlySummary
#phosMat<-fullPhosMat
#phosData<-uncorrectedPhosData
m13Phos<-list(Tram_5min=limmaTwoFactorDEAnalysis(phosMat,    
                         filter(summary,conditionName=='MOLM-13_0_none')$sample,
                        filter(summary,conditionName=='MOLM-13_5_Trametinib')$sample),
            Tram_60min=limmaTwoFactorDEAnalysis(phosMat,
                         filter(summary,conditionName=='MOLM-13_0_none')$sample,
                        filter(summary,conditionName=='MOLM-13_60_Trametinib')$sample),
              MCP1_5min=limmaTwoFactorDEAnalysis(phosMat,                                   
                         filter(summary,conditionName=='MOLM-13_0_none')$sample,
                         filter(summary,conditionName=='MOLM-13_5_MCP-1')$sample),
                MCP1_60min=limmaTwoFactorDEAnalysis(phosMat,
                        filter(summary,conditionName=='MOLM-13_0_none')$sample,                                     
                        filter(summary,conditionName=='MOLM-13_60_MCP-1')$sample),
                MCP1_tram_5min=limmaTwoFactorDEAnalysis(phosMat,
                        filter(summary,conditionName=='MOLM-13_0_none')$sample,
                        filter(summary,conditionName=='MOLM-13_5_Trametinib+MCP-1')$sample),
                MCP1_tram_60min=limmaTwoFactorDEAnalysis(phosMat,
                        filter(summary,conditionName=='MOLM-13_0_none')$sample,     
                        filter(summary,conditionName=='MOLM-13_60_Trametinib+MCP-1')$sample))
 
  p3<-doAllKSEAplots(m13Phos,phosData)
  
  sig.kin<-p3%>%
    subset(p.value<0.05)%>%
    select(Kinase.Gene,Condition)%>%distinct()
 # print(sig.kin)
  
  plotKinDat(subset(kindat,
                    Sample%in%subset(summary,CellType=='MOLM-13')$sample),
             sig.kin$Kinase.Gene,'exp1_singleTimeCourse')
  
  combPhos<-list(comb_5min=limmaTwoFactorDEAnalysis(phosMat,
                          filter(summary,conditionName=='MOLM-13_5_Trametinib+MCP-1')$sample,
                         filter(summary,conditionName%in%c('MOLM-13_5_Trametinib','MOLM-13_5_MCP-1'))$sample),
                 comb_60min=limmaTwoFactorDEAnalysis(phosMat,
                          filter(summary,conditionName=='MOLM-13_60_Trametinib+MCP-1')$sample,
                          filter(summary,conditionName%in%c('MOLM-13_60_Trametinib','MOLM-13_60_MCP-1'))$sample))
    p4<-doAllKSEAplots(combPhos,phosData)

  combDf<-do.call(rbind,lapply(names(combPhos),function(x) mutate(combPhos[[x]],condition=x)))
  
#  tramMCP=runNetworksFromDF(combDf,gene.col='featureID',
#                              weight.col='logFC',pvalue.col='P.Value',beta=5,
#                              condition.col='condition',extra.col=c('AveExpr','t','P.Value'))
#    tramMCP=runNetworksFromDF(combDf,gene.col='featureID',
#                              weight.col='logFC',pvalue.col='P.Value',beta=10,
#                              condition.col='condition',extra.col=c('AveExpr','t','P.Value'))
#  #lateHeatmap<-kseaZscoreHeatmap(list(p3,p5),'earlyTreatmentKSEAzscoreHeatmap.pdf')
  sig.kin<-p4%>%
    subset(p.value<0.05)%>%
    select(Kinase.Gene,Condition)%>%distinct()
 # print(sig.kin)
  
  plotKinDat(subset(kindat,Sample%in%subset(summary,CellType=='MOLM-13')$sample),
             sig.kin$Kinase.Gene,'exp1_comboTimeCourse')

```

### Resistant cells time course comparisons

We want to see how the resistant cells compared compare to resistant controls at both 5 minutes and 60 minutes.

```{r collect comparisons,echo=FALSE}
#####now do various comparisons


phosMat <- earlyPhosMat
summary <-earlySummary
latePhos<-list(resist_vs_parental=limmaTwoFactorDEAnalysis(phosMat,
                        filter(summary,conditionName=='MOLM-13_0_none')$sample,
                        filter(summary,conditionName=='MOLM-13 Tr Resistant_0_none')$sample),
               resist_vs_mcp1_5min=limmaTwoFactorDEAnalysis(phosMat,
                        filter(summary,conditionName=='MOLM-13 Tr Resistant_5_MCP-1')$sample,
                        filter(summary,conditionName=='MOLM-13 Tr Resistant_0_none')$sample),
              resist_vs_mcp1_60min=limmaTwoFactorDEAnalysis(phosMat,
                        filter(summary,conditionName=='MOLM-13 Tr Resistant_60_MCP-1')$sample,
                        filter(summary,conditionName=='MOLM-13 Tr Resistant_0_none')$sample))
  
  ph3<-doAllKSEAplots(latePhos)
             
#  earlyLatePh<-plotConditionsInFlow(earlyLatePhos,title='Phosphoproteomics in late',0.05)
#  ggsave('earlyLatePhos.png',earlyLatePh,width=11,height=6)
  
  latePhresdf<-do.call(rbind,lapply(names(latePhos),function(x) data.frame(latePhos[[x]],Condition=x)))
   # lateNets<-runNetworksFromDF(ph4)
  
   # tramMCP=runNetworksFromDF(latePhresdf,gene.col='featureID',
  #                            weight.col='logFC',pvalue.col='adj.P.Val',beta=5,
  #                            condition.col='Condition',extra.col=c('AveExpr','t','P.Value'))
  #      tramMCP=runNetworksFromDF(latePhresdf,gene.col='featureID',
  #                            weight.col='logFC',pvalue.col='adj.P.Val',beta=10,
  #                            condition.col='Condition',extra.col=c('AveExpr','t','P.Value'))
  
  
  sig.kin<-ph3%>%
    subset(p.value<0.05)%>%
    select(Kinase.Gene,Condition)%>%distinct()

  plotKinDat(subset(kindat,
                    Sample%in%subset(summary,CellType%in%c('MOLM-13 Tr Resistant','Late MOLM-13'))$sample),
             sig.kin$Kinase.Gene,'resistantCellLines')

```
These represent the late resistance samples


## Plot selected kinases

We have a list of kinases from Rucha that she wants to see. 

```{r kinase activity}

  kins=c('CDK1', 'CDK2', 'CDK4', 'CDK6', 'MAPK1', 'MAPK3', 'MAPK9', 'SRC', 'MTOR', 'JAK2', 'JAK3','TYK2','HCK','LCK')

 p<- kindat%>%dplyr::rename(sample='Sample')%>%left_join(rbind(earlySummary,lateSummary))%>%subset(Kinase%in%kins)%>%subset(!Treatment%in%c('Treatment Withdrawn',NA))%>%ggplot(aes(x=TimePoint,y=meanLFC,color=Treatment,shape=CellType))+geom_jitter()+facet_grid(CellType~Kinase)  
 #print(p)
 
tws<-subset(earlySummary,Treatment%in%c('Trametinib Withdrawn',NA))$sample
 
 kindat%>%
      subset(!Sample%in%tws)%>%
   plotKinDat(.,kins,prefix='selectedExp1')
 #ggsave('selectedKinaseActivity.pdf',p)
```
<!-- ### Network analysis -->
<!-- To do the network analysis we focus on the substrates only, from the un-corrected data. Skipping this now. -->

<!-- ```{r network, echo=FALSE} -->
<!-- uncorrectedPhosData<-querySynapseTable('syn24389738')%>%subset(!is.nan(LogRatio))%>% -->
<!--   mutate(Gene=unlist(Gene))%>% -->
<!--   mutate(site=unlist(site)) -->

<!-- fullPhosMat<-uncorrectedPhosData%>%dplyr::select(sample,site,LogRatio)%>% -->
<!--   tidyr::pivot_wider(values_from=LogRatio,names_from=sample, -->
<!--                      values_fn=list(LogRatio=mean))%>%#,values_fill=list(LogRatio=0.0))%>% -->
<!--   tibble::remove_rownames()%>% -->
<!--   tibble::column_to_rownames('site') -->

<!--  uncorrectedKinDat<-mapPhosphoToKinase(dplyr::rename(uncorrectedPhosData,Sample='sample', LogFoldChange='LogRatio')) -->

<!--  combPhos<-list(comb_5min=limmaTwoFactorDEAnalysis(fullPhosMat, -->
<!--                           filter(summary,conditionName=='MOLM-13_5_Trametinib+MCP-1')$sample, -->
<!--                           filter(summary,conditionName%in%c('MOLM-13_5_Trametinib','MOLM-13_5_MCP-1'))$sample), -->
<!--                  comb_60min=limmaTwoFactorDEAnalysis(fullPhosMat, -->
<!--                           filter(summary,conditionName=='MOLM-13_60_Trametinib+MCP-1')$sample, -->
<!--                           filter(summary,conditionName%in%c('MOLM-13_60_Trametinib','MOLM-13_60_MCP-1'))$sample)) -->
<!-- ##now we can load this into a network -->
<!-- ``` -->