genelist_generation_Code.R

library(tidyverse)
library(dplyr)
library(PCAtools)
library(factoextra)
library(FactoMineR)
library(grid)

setwd("/Users/pdhingra/OneDrive - Palatin Technologies, Inc/Documents/PROJECTS/DiabeticRetinopathy/proteomics/")

na_to_median <- function(x){
  x[is.na(x)] <- median(x, na.rm = TRUE)
  return(x)
}

# 1. Data loading & preprocessing ---------------------------------------------------------
sample_info <- read.table('sample_description.txt',header=TRUE) 
protein_data <- read.table("rat_retina_protein_normalized.txt",header = TRUE,check.names = FALSE)
protein_data<-protein_data[,-c(4,9,19,21,31,37,40)]

protein_data<-as.data.frame(protein_data)
rownames(protein_data)<-protein_data$ProteinName
proteins <- protein_data %>%
  select(-c(ProteinID, ProteinName)) %>%
  as.matrix()%>%t()


head(proteins)
dim(proteins)
#remove rows with less than 35 entries
#remove rows with less than 80% values meaning proteins with values for less than 35 samples
narows<-(apply(proteins,2, function(x) sum(is.na(x))))

#proteins with more than 8 NA values
toremove<-names(which(narows>8))
dim(proteins)
proteins<-proteins[,!(colnames(proteins)%in% toremove)]
dim(proteins)
proteins <- apply(proteins, 2, na_to_median) #  replace missing values
proteins <- proteins[, apply(proteins, 2, var) > 0 & apply(proteins, 2, function(y) !all(is.na(y)))] #  remove empty proteins



ggplot(data = as.data.frame(res_lda$x[, 1:2]),
       aes(x=LD1, y=LD2, colour= as.factor(pca_data$Treatment))) +
  geom_point(size = 1.3) + 
  theme_minimal() +
  scale_color_manual(values = brewer.pal(8, "Dark2")) +
  theme(axis.ticks.x = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title = element_blank(),
        legend.position = "")

#run PCA and use gene loadings to find genes which are <-0.1 and >0.1 
#select these genes from top 37 PCs 
#use these genes to plot heatmap
#see if they explain any variation
sample_info<-read.table("sample_description.txt",header=T)
proteins<-as.data.frame(t(proteins))

meta<-sample_info[match(colnames(proteins),sample_info$sample_name),]
rownames(meta)<-meta$sample_name
p<-PCAtools::pca(proteins,metadata=meta,removeVar = 0.1)
eigencorplot(p,metavars = c('sample_name','compound','Response'))
eigencorplot(p,metavars = c('sample_name','compound','Response'),components =getComponents(p, seq_len(31)) )

PC<-as.data.frame(p$loadings[,c(2,3,8,9,18,19,24)])
PC<-as.data.frame(p$loadings[,1:31])
#find top genes for each PC with coefficient >or <-0.1
genesel<-function(PC)
{
  index=which(PC < -0.1 | PC >0.1) 
  return(index)
}

genelist<-data.frame(genes=NULL)
for(i in 1:31)
{  
  final<-NULL
  final<-genesel(PC[,i])
  tmp<-data.frame(genes=rownames(PC[final,]))
  genelist<-rbind(genelist,tmp)
}  
genelist<-unique(genelist)

# perform wilcox test tp see differentialy expressed proteins between treatment groups
#read each selprotein and divide sample into categories, find significance

selprotein<-subset(proteins,rownames(proteins)%in%genelist$genes)
#read sample description
sampledesc<-read.delim("sample_description.txt",header=T,sep="\t")
sampledesc
dim(selprotein)
sink("protein_logfile.txt")
functest<-function(s1,s2,name,prot1,prot2)
{
  result<-wilcox.test(as.numeric(s1),as.numeric(s2))
  means1<-log2(mean(as.numeric(s1))) #s1- STZ
  means2<-log2(mean(as.numeric(s2))) #s2- treatment
  df<-means1-means2 #negative mean low gene expression in treatment, posiitve means up regulation
  cat(as.character(prot1)," vs ",
      as.character(prot2),"\t",
      as.character(name),"\t",
      result$p.value,"\t",
      result$statistic,"\t",
      df,"\n")
  
  
  
}



for(i in 1:nrow(selprotein))
{
  tmp<-NULL
  tmp<-as.data.frame(t(selprotein[i,]))
  # tmp<-as.data.frame(t(new_prodata[i,]))
  colnames(tmp)<-c("normprot")
  tmp$sample<-rownames(tmp)
  #tmp<-tmp[-c(1:2),]
  head(tmp)
  tmp<-merge(tmp,sampledesc,by.x=c("sample"),by.y=c("sample_name"))
  #protname=protname
  # protname=new_prodata$ProteinName[i]
  protname=rownames(selprotein[i,])
  
  #perform wilcox test
  ######################
  # STZ vs PL9654_0.05 mg/kg
  ######################
  rm(s1)
  rm(s2)
  s1<-tmp[tmp$compound=="STZ",]$normprot
  s2<-tmp[tmp$compound=="PL9654_0.05",]$normprot
  functest(s1,s2,protname,"STZ", "PL9654_0.05")
  
  ######################
  # STZ vs PL9654_0.1 mg/kg
  ######################
  rm(s1)
  rm(s2)
  s1<-tmp[tmp$compound=="STZ",]$normprot
  s2<-tmp[tmp$compound=="PL9654_0.1",]$normprot
  functest(s1,s2,protname,"STZ", "PL9654_0.1")
  
  ######################
  # STZ vs PL9654_0.5 mg/kg
  ######################
  rm(s1)
  rm(s2)
  s1<-tmp[tmp$compound=="STZ",]$normprot
  s2<-tmp[tmp$compound=="PL9654_0.5",]$normprot
  functest(s1,s2,protname,"STZ","PL9654_0.5")
  
  
  ######################
  # STZ vs healthy
  ######################
  rm(s1)
  rm(s2)
  s1<-tmp[tmp$compound=="STZ",]$normprot
  s2<-tmp[tmp$compound=="Healthy",]$normprot
  functest(s1,s2,protname,"STZ","Healthy")
  
}
sink()
closeAllConnections()

# remove PL8177 from sample_info
num=which(sample_info$compound=="PL8177_1")
sample_info<-sample_info[-num,]
sample_info[sample_info$Response=="R",]$Response="High Ranked"
sample_info[sample_info$Response=="NR",]$Response="Low Ranked"
#make complex heatmap
sub_protein<-subset(proteins,rownames(proteins)%in%genelist$genes)
dim(sub_protein)
sub_protein<-sub_protein[,match(sample_info[order(sample_info$Treatment,sample_info$Response),]$sample_name,colnames(sub_protein))]
metaheatmap<-sample_info[match(colnames(sub_protein),sample_info$sample_name),]

#make heatmap for cluster proteins
ann<-data.frame(metaheatmap$Treatment,metaheatmap$Response)
colnames(ann)<-c("Treatment","Response")

colours<-list('Treatment'=c('Healthy, naive'="#02BA26","STZ, 0.05mg/kg PL9654"="#AD07E3",
                            "STZ, 0.1mg/kg PL9654"="light pink","STZ, 0.5mg/kg PL9654"="magenta",
                            "STZ, vehicle"="black"),
              'Response'=c('STZ'='black','Healthy'='#02BA26','Low Ranked'="#ECB5FD",'High Ranked'="#CF49F9"))

colAnn<-ComplexHeatmap::HeatmapAnnotation(df=ann,which='col',col=colours)
ComplexHeatmap::Heatmap(sub_protein,top_annotation = colAnn,cluster_columns = FALSE,row_names_gp = gpar(fontsize = 8),
                        row_km=2,row_km_repeats = 100)


#boxplots
#make box plot for LEG2, FGF2, GFAP

sample_info <- read.table('sample_description.txt',header=TRUE) 
protein_data <- read.table("rat_retina_protein_normalized.txt",header = TRUE,check.names = FALSE)
protein_data<-protein_data[,-c(4,9,19,21,31,37,40)]

protein_data<-as.data.frame(protein_data)
rownames(protein_data)<-protein_data$ProteinName
protein_data <- protein_data %>%
  select(-c(ProteinID, ProteinName)) %>%
  as.matrix()%>%t()

#LEG3
leg3prot<-data.frame(protein_data[,"LEG3_RAT"])
colnames(leg3prot)<-"protein"
leg3prot$protein_name<-"LEG3_PROT"
leg3prot$sample_name<-rownames(leg3prot)

leg3prot<-merge(leg3prot,sample_info,by="sample_name")


ggplot(data=leg3prot,aes(x=protein,y=Treatment,fill=Treatment,group=Treatment))+geom_boxplot()+
  scale_fill_manual(values=c("#02BA26","#AD07E3","light pink","magenta","black"))+
  theme(panel.background = element_blank(),
        axis.title = element_text(size=14),
        axis.text=element_text(size=14,colour = "black"),legend.position = "none",
        axis.line = element_line(color="black"))+xlab("Normalized expression of LGALS3 protein")

#STAT3
Stat3prot<-data.frame(protein_data[,"STAT3_RAT"])
colnames(Stat3prot)<-"protein"
Stat3prot$protein_name<-"Stat3_PROT"
Stat3prot$sample_name<-rownames(Stat3prot)

Stat3prot<-merge(Stat3prot,sample_info,by="sample_name")


ggplot(data=Stat3prot,aes(x=protein,y=Treatment,fill=Treatment,group=Treatment))+geom_boxplot()+
  scale_fill_manual(values=c("#02BA26","#AD07E3","light pink","magenta","black"))+
  theme(panel.background = element_blank(),
        axis.title = element_text(size=14),
        axis.text=element_text(size=14,colour = "black"),legend.position = "none",
        axis.line = element_line(color="black"))+xlab("Normalized expression of LGALS3 protein")



#GFAP

Gfapprot<-data.frame(protein_data[,"GFAP_RAT"])
colnames(Gfapprot)<-"protein"
Gfapprot$protein_name<-"GFAP_PROT"
Gfapprot$sample_name<-rownames(Gfapprot)

Gfapprot<-merge(Gfapprot,sample_info,by="sample_name")


ggplot(data=Gfapprot,aes(x=protein,y=Treatment,fill=Treatment,group=Treatment))+geom_boxplot()+
  scale_fill_manual(values=c("#02BA26","#AD07E3","light pink","magenta","black"))+
  theme(panel.background = element_blank(),
        axis.title = element_text(size=18),
        axis.text=element_text(size=18,colour = "black"),legend.position = "none",
        axis.line = element_line(color="black"))+xlab("Normalized expression of GFAP protein")

#FGF2
Fgf2prot<-data.frame(protein_data[,"FGF2_RAT"])
colnames(Fgf2prot)<-"protein"
Fgf2prot$protein_name<-"FGF2_PROT"
Fgf2prot$sample_name<-rownames(Fgf2prot)

Fgf2prot<-merge(Fgf2prot,sample_info,by="sample_name")

ggplot(data=Fgf2prot,aes(x=protein,y=Treatment,fill=Treatment,group=Treatment))+geom_boxplot()+
  scale_fill_manual(values=c("#02BA26","#AD07E3","light pink","magenta","black"))+
  theme(panel.background = element_blank(),
        axis.title = element_text(size=14),
        axis.text=element_text(size=14,colour = "black"),legend.position = "none",
        axis.line = element_line(color="black"))+xlab("Normalized expression of FGF2 protein")

#B2M
B2Mprot<-data.frame(protein_data[,"B2MG_RAT"])
colnames(B2Mprot)<-"protein"
B2Mprot$protein_name<-"B2M_PROT"
B2Mprot$sample_name<-rownames(B2Mprot)

B2Mprot<-merge(B2Mprot,sample_info,by="sample_name")

ggplot(data=B2Mprot,aes(x=protein,y=Treatment,fill=Treatment,group=Treatment))+geom_boxplot()+
  scale_fill_manual(values=c("#02BA26","#AD07E3","light pink","magenta","black"))+
  theme(panel.background = element_blank(),
        axis.title = element_text(size=18),
        axis.text=element_text(size=18,colour = "black"),legend.position = "none",
        axis.line = element_line(color="black"))+xlab("Normalized expression of B2M protein")

#make final plot
leg3prot$LEG3_RAT<-NULL
leg3prot$TYPE<-"LGALS3"


Gfapprot$GFAO_RAT<-NULL
B2Mprot$B2M_RAT<-NULL
Fgf2prot$Fgf2_RAT<-NULL

Gfapprot$TYPE<-"GFAP"
B2Mprot$TYPE<-"B2M"
Fgf2prot$TYPE<-"FGF2"

final<-rbind(leg3prot,Fgf2prot,B2Mprot)


final<-rbind(Gfapprot,Fgf2prot,B2Mprot)
final$Treatment<-factor(final$Treatment,levels=c("Healthy, naive","STZ, vehicle","STZ, 0.05mg/kg PL9654","STZ, 0.1mg/kg PL9654","STZ, 0.5mg/kg PL9654"))
ggplot(data=final,aes(x=as.numeric(protein),y=Treatment,fill=Treatment,group=Treatment))+geom_boxplot()+
  scale_fill_manual(values=c("#02BA26","black","#AD07E3","light pink","magenta"))+
  theme(panel.background = element_blank(),
        axis.title = element_text(size=14),
        strip.text.x = element_text(size = 14),
        axis.text=element_text(size=14,colour = "black"),legend.position = "none",
        axis.line = element_line(color="black"))+facet_wrap(.~TYPE)+xlab("log normalized protein expression")


#Function for boxplots
boxplotprotein<-function(proname)
{
  prot<-data.frame(t(protein_data[proname,]))
  prot$sample_name<- colnames(protein_data)
  prot<-data.frame(prot[3:45,])
  colnames(prot)<-c("protein","sample_name")
  prot<-merge(prot,sample_info,by="sample_name")
  prot$protein<-as.numeric(prot$protein)
  
  p<-ggplot(data=prot,aes(x=protein,y=Treatment,fill=Treatment,group=Treatment))+geom_boxplot()+
    scale_fill_manual(values=c("blue","#CC79A7","#E7B800","#5D3A9B","purple","red"))+
    theme(panel.background = element_blank(),
          axis.title = element_text(size=14),
          axis.text=element_text(size=14,colour = "black"),legend.position = "none",
          axis.line = element_line(color="black"))+
    xlab(as.character(proname))
  return(p)
}