So, lets load RMarkdown

library(rmarkdown)

set messages to FALSE on everything (prevents certain boring things from being shown in the results)

knitr::opts_chunk$set(echo = TRUE, message=FALSE,warning=FALSE,collapse = TRUE)

PACKAGES

library(reshape2)
library(ggplot2)
library(dplyr)
library(plotly)
library(viridis)
library(data.table)
library(pheatmap)
library(tidyverse)
library(ggthemes)
library(clipr)
library(tidyr)
library(Rcpp)

COLORS

mycolors<-c(viridis(15))
felix_cols<-mycolors[c(5,2)]
felix_4cols<-mycolors[c(15,10,8,2)]
plain_cols1<-c("blue","gray")
plain_cols2<-c("red","gray")

pats_cols<-colorRampPalette(c("#FDE725FF", "white","#440154FF"))(21)
leos_cols<-colorRampPalette(c("white","blue"))(10)

LOAD DATA AND MAKE AN OVERALL HEATMAP

## load the dataset (cells treated with BREAST CANCER)
BC_data<-read_csv(file="breast_cancer_cells.csv")

## make some new columns that store the log ratios  
BC_data<-BC_data %>% mutate(
  log_MCF7=log2((MCF7_1+MCF7_2)/(MCF10A_1+MCF10A_2)), 
  log_MDA231=log2((MDA231_1+MDA231_2)/(MCF10A_1+MCF10A_2)), 
  log_MDA468=log2((MDA468_1+MDA468_2)/(MCF10A_1+MCF10A_2)),
  log_SKBR3=log2((SKBR3_1+SKBR3_2)/(MCF10A_1+MCF10A_2))) 

## make a matrix of just the data
BC_mat1<-BC_data %>% select(log_MCF7:log_SKBR3) %>% as.matrix() %>% round(.,2)
## make a heatmap from the data
pheatmap(BC_mat1, color=pats_cols,cellwidth=30,cellheight=.03,cluster_cols=FALSE,cluster_rows=FALSE,
                    legend=TRUE,fontsize = 7,scale="column")

PLOTTING A VOLCANO

##BASIC VOLCANO PLOT

## Add a column that stores the negative log of the pvalue
BC_data<-BC_data %>% mutate(neglog_SKBR3=-log10(pvalue_SKBR3_vs_MCF10A))

## Use ggplot to plot the log ratio at 24 hours against the -log p-value at 24 hours
volcano_plot<-BC_data %>% ggplot(aes(x=log_SKBR3,y=neglog_SKBR3,description=Gene_Symbol))+
  geom_point(alpha=0.7,color="blue")

volcano_plot


  
## BETTER VOLCANO PLOT

## Define the siginificant ones so they can be colored
BC_data<-BC_data %>% mutate(significance=ifelse((log_SKBR3>2 & neglog_SKBR3>2.99),"UP", ifelse((log_SKBR3<c(-2) & neglog_SKBR3>2.99),"DOWN","NOT SIG")))

## Some standard colors
plain_cols3<-c("red","gray","blue")

## volcano plot as before with some added things
better_volcano_plot<-BC_data %>% ggplot(aes(x=log_SKBR3,y=neglog_SKBR3,description=Gene_Symbol,color=significance))+
  geom_point(alpha=0.7)+
  scale_color_manual(values=plain_cols3)+
  xlim(-6,6)+
  theme_bw()+
  theme(axis.text = element_text(colour = "black",size=14))+
  theme(text = element_text(size=14))+
  labs(x="log ratio at 24 hrs COVID-19 infection compared to control",y="-log(p-value)")
  

  

better_volcano_plot

## use ggplotly to see hover over the points to see the gene names
ggplotly(better_volcano_plot)

Barplots of significant points of interest

## EXAMPLES OF A COUPLE PROTEINS or GENES

# Make a pivot longer table of the C19 data for Control_2h_1:Virus_24h_2
BC_long<-pivot_longer(BC_data, cols = c(MCF10A_1:SKBR3_2), names_to = 'variable')%>% select(-c(pvalue_MCF7_vs_MCF10A:significance))%>%select(-Description)

## Because ggplot will automatically plot alphabetically on the x-axis, we need to make sure that our timepoints will be in the proper order by adding a new column using the factor() function

BC_long$order <- as.character(BC_long$variable)
BC_long$order <- factor(BC_long$order,levels=c('MCF10A_1','MCF10A_2','MCF7_1','MCF7_2','MDA231_1','MDA231_2','MDA468_1','MDA468_2','SKBR3_1','SKBR3_2'))



## based on the gene symbols from plotly, select a few proteins from the table. Select just the data and gene symbols from the pivot longer table you just made 
Examples_Down<-BC_long %>% filter(Gene_Symbol=="APOA1" | Gene_Symbol=="HLA-A" | Gene_Symbol=="MYO1B" | Gene_Symbol=="HMGN5")
  
## make barplots facetted by Gene Symbol
Examples_Down %>% 
  ggplot(aes(x=order,y=value))+ 
  geom_bar(stat="identity",fill="red")+
  facet_wrap(~Gene_Symbol)+
  theme_bw()+
  theme(axis.text = element_text(colour = "black",size=10))+
  theme(text = element_text(size=14))+
  theme(axis.text.x = element_text(angle=45, hjust=1))+
  labs(x="cell line",y="relative intensity")

## Same process for the upregulated ones
Examples_Down<-BC_long %>% filter(Gene_Symbol=="GCLC" | Gene_Symbol=="FNBP1L" | Gene_Symbol=="DENND4C" | Gene_Symbol=="KRT23" | Gene_Symbol=="TDP2")

Examples_Down %>% 
  ggplot(aes(x=order,y=value))+ 
  geom_bar(stat="identity",fill="orange")+
  facet_wrap(~Gene_Symbol)+
  theme_bw()+
  theme(axis.text = element_text(colour = "black",size=10))+
  theme(text = element_text(size=14))+
  theme(axis.text.x = element_text(angle=45, hjust=1))+
  labs(x="cell line",y="relative intensity")