Paper1_code.Rmd

---
title: "R Notebook"
output: html_notebook
---

This is an [R Markdown](http://rmarkdown.rstudio.com) Notebook. When you execute code within the notebook, the results appear beneath the code. 

Try executing this chunk by clicking the *Run* button within the chunk or by placing your cursor inside it and pressing *Ctrl+Shift+Enter*. 

```{r}
#Clean up workspace
rm(list=ls(all=TRUE)) ##Clean workspace
```

```{r}
ActiveDir = dirname(rstudioapi::getActiveDocumentContext()$path)

#From: https://gist.github.com/stevenworthington/3178163
ipak <- function(pkg){
  new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
  if (length(new.pkg)) 
    install.packages(new.pkg, dependencies = TRUE)
  sapply(pkg, require, character.only = TRUE)
}

# Load or install packages from list
packages <- c("heatmaply","matrixStats","sunburstR","d3r","htmlwidgets","dplyr",
              "cluster","intergraph","visNetwork","igraph","ComplexHeatmap","seriation","circlize",
              "stringr")
ipak(packages)
```
```{r}
knitr::opts_knit$set(root.dir = ActiveDir)
#setwd(ActiveDir)
```


Read in all data
```{r}
RNAseq <- as.data.frame(read.csv("./Paper1_RNAseq.csv"))
ssGSEA <- as.data.frame(read.csv("./Paper1_Pathway.csv", row.names=1))
Drug <- as.data.frame(read.csv("./Paper1_Drug.csv", row.names=1))
```


```{r}
Build_from = which(str_detect(colnames(Drug),"^PIM")=="TRUE")
##Clean inputs 
dat = as.data.frame(t(Drug[which(Drug$Type == "Meta"),Build_from]))

dat = data.frame(lapply(dat, function(v) {
  return(toupper(v))
}))

df = dat %>% count(VB_Subtype)
colnames(df)[ncol(df)] <- "size"

                            # RCB = c("pCR"="palegreen","RCB-I"="palegreen4","RCB-II"="red","RCB-III"="red4"),
                             #VB_subtype=c("BL1"="darkolivegreen","BL2"="darkorange","LAR"="darkorchid","M"="deepskyblue","IM"="darkslateblue","UNS"="beige")


tree <- d3_nest(df, value_cols = "size")
sb1 <- sunburst(tree, width="100%", height=400,count = T, sortFunction = htmlwidgets::JS(
    "
    function(a,b) {
      // sort by count descending
      //   unlike the other example using data.name, value is at the top level of the object
      return b.value - a.value
    }
    "    
  )
  
)
sb1
htmltools::save_html(sb1, file = "./PDX_demographics.html")
```


Do clustering on ssGSEA
```{r}
#Select data rows and columns
Data_Columns = 2:ncol(ssGSEA)
Data_Rows = which(str_detect(row.names(ssGSEA),"^HALLMARK_")=="TRUE")

Meta_Columns = 2:5

Select_meta = c("VB_Subtype","RCB_status")
Meta_Rows= which(rownames(ssGSEA) %in% Select_meta)

#Build input dataframe from selection
Input = Drug[Data_Rows,Data_Columns]

Input = ssGSEA[Data_Rows,Data_Columns]
Input[,] <- apply(Input[,],2,function(x) as.numeric(as.character(x)))

Col_annotations = as.data.frame(t(ssGSEA[Meta_Rows,Data_Columns]))

#Calculate row StdDev
StdDev = c()
for(i in 1:nrow(Input)){
  StdDev[i] = rowSds(as.matrix(Input[i,]))
}

#Publication figure
set.seed(1536)
Anno_C = HeatmapAnnotation(RCB = Col_annotations$RCB_status,
                           VB_subtype=Col_annotations$VB_Subtype,col=list(
                             RCB = c("pCR"="palegreen","RCB-I"="palegreen4",
                                     "RCB-II"="red","RCB-III"="red4"),
                             VB_subtype=c("BL1"="darkolivegreen","BL2"="darkorange",
                                          "LAR"="darkorchid","M"="deepskyblue",
                                          "IM"="darkslateblue","UNS"="beige")
                             )
                           )
Anno_R = rowAnnotation(StDev=StdDev)

o1 = seriate(dist(Input), method = "GW_complete")
o2 = seriate(dist(t(Input)), method = "GW_complete")

#Input = t(scale(t(Input),scale=TRUE,center=TRUE))
#pdf("Hallmarks_20200524.pdf", width = 12, height = 8)
ht = Heatmap(Input,
        clustering_distance_columns = function(m) dist(m, method = "euclidean"),
        clustering_method_columns = "complete",
        clustering_distance_rows = function(m) dist(m, method = "euclidean"),
        clustering_method_rows = "complete",
        cluster_rows = as.dendrogram(o1[[1]]),
        cluster_columns = as.dendrogram(o2[[1]]),
        show_row_names = TRUE,
        top_annotation = Anno_C,
        left_annotation = Anno_R,
        col = colorRamp2(c(min(as.matrix(Input)),max(as.matrix(Input))),c("blue","orange")),
        row_names_max_width = max_text_width(rownames(ssGSEA),gp = gpar(fontsize = 12)))
draw(ht,heatmap_legend_side = "left",annotation_legend_side = "left")
#dev.off()
```


```{r}
Data_Columns = 6:ncol(Drug)
Data_Rows = which(Drug$Type == "Exp")

Meta_Columns = 2:5

Select_meta = c("VB_Subtype")
Meta_Rows= which(rownames(Drug) %in% Select_meta)

#Meta_Rows = which(RNAseq$Type == "Meta")


Input = Drug[Data_Rows,Data_Columns]

Input[,] <- apply(Input[,],2,function(x) as.numeric(as.character(x)))

StdDev = c()
for(i in 1:nrow(Input)){
  StdDev[i] = rowSds(as.matrix(Input[i,]))
}
top_N = which(StdDev>=min(tail(sort(StdDev),1000)))

Input = Input[top_N,]
Drug_temp = as.data.frame(t(Input))

Row_annotations = Drug[which(rownames(Drug) %in% rownames(Input)),Meta_Columns]
Col_annotations = t(Drug[Meta_Rows,Data_Columns])


heatmaply(Input,
          hclust_method = "complete",
          dist_method = "euclidean",
          seriate = "GW", 
          colors = cool_warm(25),
          k_col = NA,
          k_row = NA,
          titleX = FALSE, 
          titleY = FALSE,
          #plot_method = "plotly",
          RowSideColors = Row_annotations,
          col_side_colors = Col_annotations,
          showticklabels=c(TRUE,FALSE))

```


Clustering for the top X most variant drugs
```{r}
Data_Columns = 6:ncol(Drug)
Data_Rows = which(Drug$Type == "Exp")

Meta_Columns = 2:5

Select_meta = c("VB_Subtype","RCB_status")
Meta_Rows= which(rownames(Drug) %in% Select_meta)

#Meta_Rows = which(RNAseq$Type == "Meta")


Input = Drug[Data_Rows,Data_Columns]

Input[,] <- apply(Input[,],2,function(x) as.numeric(as.character(x)))

StdDev = c()
for(i in 1:nrow(Input)){
  StdDev[i] = rowSds(as.matrix(Input[i,]))
}
top_N = which(StdDev>=min(tail(sort(StdDev),1000)))

Input = Input[top_N,]
Drug_temp = as.data.frame(t(Input))

Row_annotations = Drug[which(rownames(Drug) %in% rownames(Input)),Meta_Columns]
Col_annotations = as.data.frame(t(Drug[Meta_Rows,Data_Columns]))

#Publication figure
set.seed(1536)
Anno_C = HeatmapAnnotation(RCB = Col_annotations$RCB_status,
                           VB_subtype=Col_annotations$VB_Subtype,col=list(
                             RCB = c("pCR"="palegreen","RCB-I"="palegreen4","RCB-II"="red","RCB-III"="red4"),
                             VB_subtype=c("BL1"="darkolivegreen","BL2"="darkorange","LAR"="darkorchid","M"="deepskyblue","IM"="darkslateblue","UNS"="beige")
                             )
                           )
Anno_R = HeatmapAnnotation(df=Row_annotations)

o1 = seriate(dist(Input), method = "GW_complete")
o2 = seriate(dist(t(Input)), method = "GW_complete")

#pdf("Fig_3B.pdf", width = 12, height = 8)
Heatmap(Input,
        clustering_distance_columns = function(m) dist(m, method = "euclidean"),
        clustering_method_columns = "complete",
        clustering_distance_rows = function(m) dist(m, method = "euclidean"),
        clustering_method_rows = "complete",
        cluster_rows = as.dendrogram(o1[[1]]),
        cluster_columns = as.dendrogram(o2[[1]]),
        show_row_names = FALSE,
        top_annotation = Anno_C,
        col = colorRamp2(c(0,1),c("blue","orange")))
#dev.off
```


Do log2(TPM to drug correlations)
```{r}
Data_Columns = 6:ncol(Drug)
Data_Rows = which(Drug$Type == "Exp")
Meta_Columns = 2:5
Select_meta = c("VB_Subtype","RCB_status")
Meta_Rows= which(rownames(Drug) %in% Select_meta) 
Input = Drug[Data_Rows,Data_Columns]
Drug_temp = as.data.frame(t(Input))

Data_Columns = 3:ncol(RNAseq)
Data_Rows = which(RNAseq$Type == "Exp")
Meta_Columns = 1
Select_meta = c("VB_Subtype","RCB_status")
Meta_Rows= which(rownames(RNAseq) %in% Select_meta)
Input = RNAseq[Data_Rows,Data_Columns]

names = as.vector(RNAseq[which(RNAseq$Sample=="PIM_ID"),Data_Columns])
names[,] <- apply(names[,],2,function(x) as.character(x))
colnames(Input) = names

rownames(Input) = RNAseq$Sample[Data_Rows]
Input[,] <- apply(Input[,],2,function(x) log2(as.numeric(as.character(x))))
Input[mapply(is.infinite, Input)] <- 0

Pathway_temp = as.data.frame(t(Input))

Spear_R = data.frame()
Spear_P = data.frame()

for(i in 1:ncol(Pathway_temp))
{
  for(j in 1:ncol(Drug_temp))
    {
      Active_Gene_name = colnames(Pathway_temp)[i]
      Active_Gene = as.data.frame(Pathway_temp[,i])
      row.names(Active_Gene) = names
      
      Active_Drug_name = colnames(Drug_temp)[j]
      Active_Drug = as.data.frame(Drug_temp[,j])
      
      Cor_matrix = merge(Active_Gene,Active_Drug,by="row.names")
      Cor_matrix[,c(2,3)] <- apply(Cor_matrix[,c(2,3)],2,function(x) as.numeric(as.character(x)))
      
      test = cor.test(Cor_matrix[,2],Cor_matrix[,3],method="pearson",verbose=FALSE)
      Spear_R[Active_Drug_name,Active_Gene_name] = test$estimate
      Spear_P[Active_Drug_name,Active_Gene_name] = test$p.value
    }
}
```

```{r}
Data_Columns = 6:ncol(Drug)
Data_Rows = which(Drug$Type == "Exp")
Meta_Columns = 2:5

Select_meta = c("VB_Subtype")
Meta_Rows= which(rownames(Drug) %in% Select_meta) 
Input = Drug[Data_Rows,Data_Columns]
Drug_temp = as.data.frame(t(Input))

Data_Columns = 2:ncol(ssGSEA)
Data_Rows = which(ssGSEA$Type == "Exp")
Input = ssGSEA[Data_Rows,Data_Columns]

# names = as.vector(colnames(ssGSEA[,Data_Columns]))
# Input[,] <- apply(Input[,],2,function(x) log2(as.numeric(as.character(x))))
Input[mapply(is.infinite, Input)] <- 0
Input[mapply(is.na, Input)] <- 0

Pathway_temp = as.data.frame(t(Input))

Spear_R = data.frame()
Spear_P = data.frame()

for(i in 1:4425) #ncol())
{
  for(j in 1:637) #ncol())
    {
      Active_Gene_name = colnames(Pathway_temp)[i]
      Active_Gene = as.data.frame(Pathway_temp[,i])
      row.names(Active_Gene) = names
      
      Active_Drug_name = colnames(Drug_temp)[j]
      Active_Drug = as.data.frame(Drug_temp[,j])
      
      Cor_matrix = merge(Active_Gene,Active_Drug,by="row.names")
      Cor_matrix[,c(2,3)] <- apply(Cor_matrix[,c(2,3)],2,function(x) as.numeric(as.character(x)))
      
      test = cor.test(Cor_matrix[,2],Cor_matrix[,3],method="pearson",verbose=FALSE)
      Spear_R[Active_Drug_name,Active_Gene_name] = test$estimate
      Spear_P[Active_Drug_name,Active_Gene_name] = test$p.value
    }
}

temp = cbind(Drug[which(Drug$Type == "Exp"),1:5],Spear_P)
write.csv(temp,"./Pathway_to_Drug_P.csv",row.names=TRUE)

temp = cbind(Drug[which(Drug$Type == "Exp"),1:5],Spear_R)
write.csv(temp,"./Pathway_to_Drug_R.csv",row.names=TRUE)
```

"Manual" exploration between drug and ssGSEA pathways analysis
```{r}
Search_drug = "PRIMA-1MET (Broad_2)"
Search_pathway = "MDM2"

Drug_temp = t(Drug[Search_drug,])
Path_temp = t(RNAseq[which(RNAseq$Sample == Search_pathway),3:ncol(RNAseq)])

Path_temp[] <- apply(Path_temp[],1,function(x) log2(as.numeric(as.character(x))))
Path_temp[mapply(is.infinite, Path_temp)] <- 0

names = as.vector(RNAseq[which(RNAseq$Sample=="PIM_ID"),3:ncol(RNAseq)])
names[,] <- apply(names[,],2,function(x) as.character(x))
row.names(Path_temp) = names

df_temp = merge(Drug_temp,Path_temp,by="row.names")
df_temp = df_temp[which(str_detect(df_temp$Row.names,"PIM")=="TRUE"),]
df_temp[,c(2,3)] <- apply(df_temp[,c(2,3)],2,function(x) as.numeric(as.character(x)))
colnames(df_temp) = c("Name","Drug","Pathway")

ggplot(df_temp,aes(x=Drug, y=Pathway))+geom_point()+geom_smooth(method=lm)#+xlim(-1,1)+ylim(-1,1)+xlab(Search_drug)+ylab(Search_pathway)
cor.test(df_temp$Drug,df_temp$Pathway,method="pearson")
```


```{r}
Search_drug = "PRIMA-1MET (Broad_2)"
Search_pathway = "REACTOME_GLUTATHIONE_CONJUGATION"

Drug_temp = t(Drug[Search_drug,])
Path_temp = t(ssGSEA[Search_pathway,2:ncol(ssGSEA)])

#Path_temp[] <- apply(Path_temp[],1,function(x) #log2(as.numeric(as.character(x))))
#Path_temp[mapply(is.infinite, Path_temp)] <- 0

#names = as.vector(RNAseq[which(RNAseq$Sample=="PIM_ID"),3:ncol(RNAseq)])
#names[,] <- apply(names[,],2,function(x) as.character(x))
#row.names(Path_temp) = names

df_temp = merge(Drug_temp,Path_temp,by="row.names")
df_temp = df_temp[which(str_detect(df_temp$Row.names,"PIM")=="TRUE"),]
df_temp[,c(2,3)] <- apply(df_temp[,c(2,3)],2,function(x) as.numeric(as.character(x)))
colnames(df_temp) = c("Name","Drug","Pathway")

ggplot(df_temp,aes(x=Drug, y=Pathway))+geom_point()+geom_smooth(method=lm)+ xlab(Search_drug)+ylab(Search_pathway)
cor.test(df_temp$Drug,df_temp$Pathway,method="pearson")
```

```{r}
Data_Columns = 6:ncol(Drug)
Data_Rows = which(Drug$Target %in% c("EGFR","Survivin","CDK4","HSP70","HSP90"))

Meta_Columns = 2:5

Select_meta = c("VB_Subtype","RCB_status")
Meta_Rows= which(rownames(Drug) %in% Select_meta)

Input = Drug[Data_Rows,Data_Columns]

Input[,] <- apply(Input[,],2,function(x) as.numeric(as.character(x)))

Row_annotations = Drug[which(rownames(Drug) %in% rownames(Input)),Meta_Columns]
Col_annotations = t(Drug[Meta_Rows,Data_Columns])

#Publication figure
set.seed(1536)
Anno_C = HeatmapAnnotation(df=Col_annotations)
Anno_R = HeatmapAnnotation(df=Row_annotations)

o1 = seriate(dist(Input), method = "GW_complete")
o2 = seriate(dist(t(Input)), method = "GW_complete")

Heatmap(Input,
        clustering_distance_columns = function(m) dist(m, method = "euclidean"),
        clustering_method_columns = "complete",
        clustering_distance_rows = function(m) dist(m, method = "euclidean"),
        clustering_method_rows = "complete",
        cluster_rows = as.dendrogram(o1[[1]]),
        cluster_columns = as.dendrogram(o2[[1]]),
        show_row_names = TRUE,
        top_annotation = Anno_C,
        col = colorRamp2(c(0,1),c("blue","orange")),
        row_names_max_width = max_text_width(rownames(Input),gp = gpar(fontsize = 8)))
```


Generate heatmap for target mitotic spindle figure
```{r}
Data_Columns = 6:ncol(Drug)
Data_Rows = which(Drug$Target %in% c("KSP","Tubulin/micotubules","CENP1","AURK","Eg5"))

Meta_Columns = 2:5

Select_meta = c("VB_Subtype","RCB_status")
Meta_Rows= which(rownames(Drug) %in% Select_meta)

Input = Drug[Data_Rows,Data_Columns]

Input[,] <- apply(Input[,],2,function(x) as.numeric(as.character(x)))

Row_annotations = as.data.frame(Drug[which(rownames(Drug) %in% rownames(Input)),Meta_Columns])
Col_annotations = as.data.frame(t(Drug[Meta_Rows,Data_Columns]))

#Publication figure
set.seed(1536)
Anno_C = HeatmapAnnotation(RCB = Col_annotations$RCB_status,
                           VB_subtype=Col_annotations$VB_Subtype,col=list(
                             RCB = c("pCR"="palegreen","RCB-I"="palegreen4","RCB-II"="red","RCB-III"="red4"),
                             VB_subtype=c("BL1"="darkolivegreen","BL2"="darkorange","LAR"="darkorchid","M"="deepskyblue","IM"="darkslateblue","UNS"="beige")
                             )
                           )
Anno_R = rowAnnotation(Target=Row_annotations$Target,col=list(Target=c("AURK"="navy","Eg5"="lightsalmon","KSP"="green","Tubulin/micotubules"="hotpink")))

o1 = seriate(dist(Input), method = "GW_complete")
o2 = seriate(dist(t(Input)), method = "GW_complete")

pdf("Mitotic_inhibitor_panel.pdf", width = 12, height = 8)
Heatmap(Input,
        clustering_distance_columns = function(m) dist(m, method = "euclidean"),
        clustering_method_columns = "complete",
        clustering_distance_rows = function(m) dist(m, method = "euclidean"),
        clustering_method_rows = "complete",
        cluster_rows = as.dendrogram(o1[[1]]),
        cluster_columns = as.dendrogram(o2[[1]]),
        show_row_names = TRUE,
        top_annotation = Anno_C,
        left_annotation = Anno_R,
        col = colorRamp2(c(0,1),c("blue","orange")),
        row_names_max_width = max_text_width(rownames(ssGSEA),gp = gpar(fontsize = 8)))
dev.off()
```




```{r}
Data_Columns = 6:ncol(Drug)
Data_Rows = which(Drug$Target %in% c("p53"))

Meta_Columns = 3

Select_meta = c("p53_status","VB_Subtype","RCB_status")
Meta_Rows= which(rownames(Drug) %in% Select_meta)

Input = Drug[Data_Rows,Data_Columns]
Input[,] <- apply(Input[,],2,function(x) as.numeric(as.character(x)))

Row_annotations = Drug[which(rownames(Drug) %in% rownames(Input)),Meta_Columns]
Col_annotations = as.data.frame(t(Drug[Meta_Rows,Data_Columns]))

#Publication figure
set.seed(1536)
Anno_C = HeatmapAnnotation(RCB = Col_annotations$RCB_status,
                           VB_subtype=Col_annotations$VB_Subtype,
                           p53_status = Col_annotations$p53_status,
                           col=list(
                             RCB = c("pCR"="palegreen","RCB-I"="palegreen4","RCB-II"="red","RCB-III"="red4"),
                             VB_subtype=c("BL1"="darkolivegreen","BL2"="darkorange","LAR"="darkorchid","M"="deepskyblue","IM"="darkslateblue","UNS"="beige"),
                             p53_status=c("Wt"="azure3","Nonsyn_SNV"="royalblue3","Stopgain"="palegreen2")
                             )
                           )
Anno_R = rowAnnotation(df=Row_annotations)

o1 = seriate(dist(Input), method = "GW_complete")
o2 = seriate(dist(t(Input)), method = "GW_complete")

pdf("P53_heatmap.pdf", width = 12, height = 8)
Heatmap(Input,
        clustering_distance_columns = function(m) dist(m, method = "euclidean"),
        clustering_method_columns = "complete",
        clustering_distance_rows = function(m) dist(m, method = "euclidean"),
        clustering_method_rows = "complete",
        cluster_rows = as.dendrogram(o1[[1]]),
        cluster_columns = TRUE,
        show_row_names = TRUE,
        top_annotation = Anno_C,
        #left_annotation = Anno_R,
        col = colorRamp2(c(0,0.4),c("blue","orange")),
        use_raster = FALSE,
        row_names_max_width = max_text_width(rownames(Drug),gp = gpar(fontsize = 12)))
dev.off()
```

Drugs that are different between subtype
```{r}
Data_Columns = which(t(Drug["VB_Subtype",]) %in% c("BL1","BL2","LAR","M"))
Data_Rows = which(Drug$Type == "Exp")

Meta_Columns = 3

Select_meta = c("VB_Subtype","RCB_status")
Meta_Rows= which(rownames(Drug) %in% Select_meta)

Input = Drug[Data_Rows,Data_Columns]

AOV_P = c()
for(i in 1:nrow(Input))
{
  df.test = as.data.frame(t(rbind(Input[i,],Drug["VB_Subtype",Data_Columns])))
  colnames(df.test) = c("Drug_val","Factor")
  df.test = transform(df.test,Drug_val=as.numeric(as.character(Drug_val)))
  result = summary(aov(Drug_val~Factor,df.test))
  AOV_P[i] = result[[1]][["Pr(>F)"]][1]
}

Input = Input[which(AOV_P<0.05),]
Input[,] <- apply(Input[,],2,function(x) as.numeric(as.character(x)))

Row_annotations = as.data.frame(Drug[which(rownames(Drug) %in% rownames(Input)),Meta_Columns])
colnames(Row_annotations) = "Targets"
Col_annotations = as.data.frame(t(Drug[Meta_Rows,Data_Columns]))

#Publication figure
set.seed(1536)
Anno_C = HeatmapAnnotation(RCB = Col_annotations$RCB_status,
                           VB_subtype=Col_annotations$VB_Subtype,col=list(
                             RCB = c("pCR"="palegreen","RCB-I"="palegreen4","RCB-II"="red","RCB-III"="red4"),
                             VB_subtype=c("BL1"="darkolivegreen","BL2"="darkorange","LAR"="darkorchid","M"="deepskyblue","IM"="darkslateblue","UNS"="beige"),
                             )
                           )
Anno_R = rowAnnotation(df=Row_annotations)

o1 = seriate(dist(Input), method = "GW_complete")
o2 = seriate(dist(t(Input)), method = "GW_complete")

pdf("Subtype_selective.pdf", width = 12, height = 8)
Heatmap(Input,
        clustering_distance_columns = function(m) dist(m, method = "euclidean"),
        clustering_method_columns = "complete",
        clustering_distance_rows = function(m) dist(m, method = "euclidean"),
        clustering_method_rows = "complete",
        cluster_rows = as.dendrogram(o1[[1]]),
        cluster_columns = FALSE,
        show_row_names = TRUE,
        top_annotation = Anno_C,
        left_annotation = Anno_R,
        col = colorRamp2(c(0,1),c("blue","orange")),
        use_raster = FALSE,
        row_names_max_width = max_text_width(rownames(Drug),gp = gpar(fontsize = 12)))
dev.off()
```


```{r}
Data_Columns = which(t(Drug["VB_Subtype",]) %in% c("BL1","BL2","LAR","M"))
Data_Rows = which(Drug$Type == "Exp")

Meta_Columns = 3

Select_meta = c("VB_Subtype","RCB_status","p53_status")
Meta_Rows= which(rownames(Drug) %in% Select_meta)

Input = Drug[Data_Rows,Data_Columns]

AOV_P = c()
for(i in 1:nrow(Input))
{
  df.test = as.data.frame(t(rbind(Input[i,],Drug["p53_status",Data_Columns])))
  colnames(df.test) = c("Drug_val","Factor")
  df.test = transform(df.test,Drug_val=as.numeric(as.character(Drug_val)))
  result = summary(aov(Drug_val~Factor,df.test))
  AOV_P[i] = result[[1]][["Pr(>F)"]][1]
}

Input = Input[which(AOV_P<0.05),]
Input[,] <- apply(Input[,],2,function(x) as.numeric(as.character(x)))

Row_annotations = as.data.frame(Drug[which(rownames(Drug) %in% rownames(Input)),Meta_Columns])
colnames(Row_annotations) = "Targets"
Col_annotations = as.data.frame(t(Drug[Meta_Rows,Data_Columns]))

#Publication figure
set.seed(1536)
# Anno_C = HeatmapAnnotation(RCB = Col_annotations$RCB_status,
#                            VB_subtype=Col_annotations$VB_Subtype,col=list(
#                              RCB = c("pCR"="palegreen","RCB-I"="palegreen4","RCB-II"="red","RCB-III"="red4"),
#                              VB_subtype=c("BL1"="darkolivegreen","BL2"="darkorange","LAR"="darkorchid","M"="deepskyblue","IM"="darkslateblue","UNS"="beige"),
#                              )
#                            )
Anno_R = rowAnnotation(df=Row_annotations)

o1 = seriate(dist(Input), method = "GW_complete")
o2 = seriate(dist(t(Input)), method = "GW_complete")

pdf("p53_selective.pdf", width = 12, height = 8)
Heatmap(Input,
        clustering_distance_columns = function(m) dist(m, method = "euclidean"),
        clustering_method_columns = "complete",
        clustering_distance_rows = function(m) dist(m, method = "euclidean"),
        clustering_method_rows = "complete",
        cluster_rows = as.dendrogram(o1[[1]]),
        cluster_columns = FALSE,
        show_row_names = TRUE,
        top_annotation = Anno_C,
        left_annotation = Anno_R,
        col = colorRamp2(c(0,1),c("blue","orange")),
        use_raster = FALSE,
        row_names_max_width = max_text_width(rownames(Drug),gp = gpar(fontsize = 12)))
dev.off()
```


```{r}
Data_Columns = 6:14
Data_Rows = which(Drug$Type == "Exp")

Meta_Columns = 3

Select_meta = c("VB_Subtype","RCB_status")
Meta_Rows= which(rownames(Drug) %in% Select_meta)

Input = Drug[Data_Rows,Data_Columns]

AOV_P = c()
for(i in 1:nrow(Input))
{
  df.test = as.data.frame(t(rbind(Input[i,],Drug["VB_Subtype",Data_Columns])))
  colnames(df.test) = c("Drug_val","Factor")
  df.test = transform(df.test,Drug_val=as.numeric(as.character(Drug_val)))
  result = summary(aov(Drug_val~Factor,df.test))
  AOV_P[i] = result[[1]][["Pr(>F)"]][1]
}

Input = Input[which(AOV_P<0.05),]
Input[,] <- apply(Input[,],2,function(x) as.numeric(as.character(x)))

Row_annotations = Drug[which(rownames(Drug) %in% rownames(Input)),Meta_Columns]
Col_annotations = t(Drug[Meta_Rows,Data_Columns])

#Publication figure
set.seed(1536)
Anno_C = HeatmapAnnotation(df=Col_annotations)
Anno_R = rowAnnotation(df=Row_annotations)

o1 = seriate(dist(Input), method = "GW_complete")
o2 = seriate(dist(t(Input)), method = "GW_complete")

pdf("Fig_4B.pdf", width = 12, height = 8)
ht = Heatmap(Input,
        clustering_distance_columns = function(m) dist(m, method = "euclidean"),
        clustering_method_columns = "complete",
        clustering_distance_rows = function(m) dist(m, method = "euclidean"),
        clustering_method_rows = "complete",
        cluster_rows = as.dendrogram(o1[[1]]),
        cluster_columns = as.dendrogram(o2[[1]]),
        show_row_names = TRUE,
        top_annotation = Anno_C,
        left_annotation = Anno_R,
        col = colorRamp2(c(0,0.5),c("blue","orange")),
        use_raster = FALSE,
        row_names_max_width = max_text_width(rownames(Drug),gp = gpar(fontsize = 12)))

draw(ht,heatmap_legend_side = "left",annotation_legend_side = "left")
dev.off()
```

```{r}
summary(Row_annotations)
```


```{r}
save.image()
```


The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike *Knit*, *Preview* does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.