15_exp2_cytokineResistancePhospho.Rmd

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


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



## 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}

fmean<-function(x){
  mean(x,na.rm=T)
}
###load all the data
sync<-synapseLogin()
phosData<-readr::read_csv(sync$get('syn51753813')$path)|>#querySynapseTable('syn22986341')%>%subset(!is.nan(LogRatio
  subset(site!='')|>
 # mutate(Gene=unlist(Gene))%>%
#  mutate(site=unlist(site))%>%
  subset(TimePoint!=120)

protData<-readr::read_csv(sync$get('syn51753816')$path)|>#querySynapseTable('syn22986326')%>%
  subset(!is.nan(LogRatio))%>%
  subset(Gene!='')|>
  #mutate(Gene=unlist(Gene))%>%
  subset(TimePoint!=120)


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

##what are we doing again?

protMat<-protData%>%
  # subset(Batch=='Experiment 1')%>%
    subset(TimePoint==0)%>%
  dplyr::select(sample,Gene,LogRatio)%>%
  distinct()|>
  tidyr::pivot_wider(values_from=LogRatio,names_from=sample,
                     values_fn=list(LogRatio=fmean))%>%
  as.data.frame()|>
  subset(Gene!='')|>
  tibble::remove_rownames()|>
  tibble::column_to_rownames('Gene')

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


earlyPhosMat<-phosData%>%
  subset(Batch=='Experiment 1')%>%
  dplyr::select(sample,site,LogRatio)%>%
  tidyr::pivot_wider(values_from=LogRatio,names_from=sample,
                     values_fn=list(LogRatio=fmean))%>%#updated 9/6,values_fill=list(LogRatio=0.0))%>%
    as.data.frame()|>
  subset(site!='')|>
  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=fmean))%>%#updated 9/6,values_fill=list(LogRatio=0.0))%>%
  as.data.frame()|>
  subset(site!='')|>
  tibble::remove_rownames()%>%
  tibble::column_to_rownames('site')

lateProtMat<-protData%>%
    subset(Batch=='Experiment 2')%>%
    subset(TimePoint==0)%>%
  dplyr::select(sample,Gene,LogRatio)%>%
  tidyr::pivot_wider(values_from=LogRatio,names_from=sample,
                     values_fn=list(LogRatio=fmean))%>%
    as.data.frame()|>
    subset(Gene!='')|>
  tibble::remove_rownames()|>
  tibble::column_to_rownames('Gene')
```

## Experiment 2 data


We can plot the data and the kinase activity as needed in heatmap or PCA form. Here we want to visualize only the samples measured in Experiment 2 (minute the 120 min time point)
```{r processing functions,echo=TRUE,warning=FALSE}

#' @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)
    #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,CellType)%>%
    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)
  #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,
          file=paste0(prefix,'cytokineKinaseHeatmap.pdf'),height=14,width=8) 
  
  res<- 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,cluster_cols=FALSE,
          clustering_distance_rows = 'correlation',
          annotation_col=sampAts,annotation_colors = sampCols,
          file=paste0(prefix,'cytokineKinaseHeatmapUnclustered.pdf'),height=14,width=8) 
  return(res)
}

####Show kinase activity
####Show kinase activity
kindat<-phosData%>%
  subset(Batch=='Experiment 2')%>%
  dplyr::rename(Sample='sample', LogFoldChange='LogRatio')%>%
  mapPhosphoToKinase()
##plot PCA
##plot PCA
phosData%>%
  subset(Batch=='Experiment 2')%>%
 # subset(CellType=='MOLM-13')%>%
  plotAllData(.,colour='Treatment',fill='Treatment')+
  ggtitle("Phosphoproteomics Experiment 2")


##plot PCA
protData%>%
  subset(Batch=='Experiment 2')%>%
 # subset(CellType=='MOLM-13')%>%
  plotAllData(.,colour='Treatment',fill='Treatment')+
  ggtitle("Proteomics Experiment 2")
```


## KSEA enrichment analysis

Next we define functions to carry out various KSEA analyses. 
```{r ksea, echo=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.1,linkedSubs=1,suffix='pdf')%>%
      mutate(Condition=clName)%>%
      as.data.frame()
  })
  return(full.df)
  
}

#' plot all the GO 
#' @param condList
#' @return data frame
doAllGOplots<-function(condList){
  
  full.df<-purrr::map_df(names(condList),.f=function(clName){ 
    condList[[clName]]%>%
      tibble::rownames_to_column('Gene')%>%
      dplyr::select(Gene,value='logFC')%>%
      amlresistancenetworks::plotOldGSEA(.,prefix=clName,0.1)%>%
      as.data.frame()
  })
  
}

#' 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))
  }
   
}
```

### Global proteomic comparisons

We want to identify global changes between the combination treatment 

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

#phosMat <- protMat
summary <- lateSummary
#phosMat<-fullPhosMat
#phosData<-uncorrectedPhosData
m13Prot<-list(Tram_vs_control=limmaTwoFactorDEAnalysis(protMat,
                         filter(summary,conditionName=='MOLM-13_0_none')$sample,    
                        filter(summary,conditionName=='Late MOLM-13_0_Trametinib')$sample),
            TramMCP1_vs_control=limmaTwoFactorDEAnalysis(protMat,
                         filter(summary,conditionName=='MOLM-13_0_none')$sample,
                        filter(summary,conditionName=='Late MOLM-13_0_Trametinib+MCP-1')$sample),
              TramMCP1_vs_Tram=limmaTwoFactorDEAnalysis(protMat,
                         filter(summary,conditionName=='Late MOLM-13_0_Trametinib')$sample,             
                         filter(summary,conditionName=='Late MOLM-13_0_Trametinib+MCP-1')$sample),
            Resist_vs_control=limmaTwoFactorDEAnalysis(protMat,
              filter(earlySummary,conditionName=='MOLM-13_0_none')$sample,
              filter(earlySummary,conditionName=='MOLM-13 Tr Resistant_0_none')$sample))
 
  p3<-doAllGOplots(m13Prot)
  

```


### Phosphoproteomic comparisons

We want to identify phosphoproteomic changes between the combination treatment 

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

phosMat <- latePhosMat
summary <- lateSummary
#phosMat<-fullPhosMat
#phosData<-uncorrectedPhosData
m13Phos<-list(Tram_vs_control=limmaTwoFactorDEAnalysis(phosMat,
                         filter(summary,conditionName=='MOLM-13_0_none')$sample,    
                        filter(summary,conditionName=='Late MOLM-13_0_Trametinib')$sample),
            TramMCP1_vs_control=limmaTwoFactorDEAnalysis(phosMat,
                         filter(summary,conditionName=='MOLM-13_0_none')$sample,
                        filter(summary,conditionName=='Late MOLM-13_0_Trametinib+MCP-1')$sample),
              TramMCP1_vs_Tram=limmaTwoFactorDEAnalysis(phosMat,
                         filter(summary,conditionName=='Late MOLM-13_0_Trametinib')$sample,             
                         filter(summary,conditionName=='Late MOLM-13_0_Trametinib+MCP-1')$sample))
 
  p3<-doAllKSEAplots(m13Phos,phosData)
  
  sig.kin<-p3%>%
    subset(p.value<0.1)%>%
    select(Kinase.Gene,Condition)%>%distinct()
 # print(sig.kin)
  
  res=plotKinDat(subset(kindat,Sample%in%summary$sample), sig.kin$Kinase.Gene,'exp2Results')
  res  
```
### Plot kinases of interest

Here we are plotting kinases of interest by expression level

```{r plot kinases}

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

#tws<-subset(earlySummary,Treatment%in%c('Trametinib Withdrawn',NA))$sample
 
 kindat%>%
 #     subset(!Sample%in%tws)%>%
   plotKinDat(.,kins,prefix='selectedExp2')
 
 p<- protData%>%
   subset(Gene%in%kins)%>%
   subset(Batch=='Experiment 2')%>%
   subset(!Treatment%in%c('Trametinib Withdrawn',NA))%>%
   ggplot(aes(x=Treatment,y=LogRatio,color=Treatment,shape=CellType))+geom_jitter()+facet_grid(CellType~Gene)+scale_x_discrete(guide = guide_axis(angle = 90))
 print(p)
 ggsave('selectedExp2KinaseExpression.pdf',p,width=10)
```

### Combined phospho and protein

combine into network analysisi

```{r network analysis, message=FALSE,echo=FALSE}

  combDf<-do.call(rbind,lapply(names(m13Prot),function(x) mutate(m13Prot[[x]],condition=x)))

#  tramMCP=runNetworksFromDF(combDf,gene.col='featureID',
#                              weight.col='logFC',pvalue.col='adj.P.Val',beta=5,
#                              condition.col='condition',extra.col=c('AveExpr','t','P.Value'))
```