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)
COLORS (feel free to mess around and find your own favorites if you like)
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
breast_cancer<-read_csv(file="breast_cancer_cells.csv")
## MCF10A cell line is non-tumorigenic (therefore control)
## make some new columns that store the log ratios
breast_cancer<-breast_cancer %>% 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)))%>%arrange(-log_MCF7)
## make a matrix of just the data
BC_mat1<-breast_cancer %>% select(MCF10A_1:SKBR3_2) %>% as.matrix() %>% round(.,2)
## make a heatmap from the data (NOTE: this is set to cluster the rows. You can remove that if you like)
pheatmap(BC_mat1, color=pats_cols,cellwidth=30,cellheight=.06,cluster_cols=FALSE,cluster_rows=TRUE,legend=TRUE,fontsize = 7,scale="row")
PLOTTING A VOLCANO
##BASIC VOLCANO PLOT
## Add a column that stores the negative log of the pvalue between 2 of the treatments or cell lines
breast_cancer<-breast_cancer %>% mutate(neglog_MCF=-log10(pvalue_MCF7_vs_MCF10A))
## Use ggplot to plot the log ratio at 24 hours against the -log p-value at 24 hours
volcano_plot<-breast_cancer %>% ggplot(aes(x=log_MCF7,y=neglog_MCF,description=Gene_Symbol))+geom_point(alpha=0.7,color="blue")
volcano_plot
BETTER VOLCANO PLOT
## Define the significant ones so they can be colored
breast_cancer<-breast_cancer %>% mutate(significance=ifelse((log_MCF7>2& neglog_MCF>2.99),"UP",ifelse((log_MCF7<c(-2)& neglog_MCF>2.99),"DOWN","NOT SIG")))
## Define some standard colors
plain_cols3<-c("red","gray","blue")
## volcano plot as before with some added things
optional_volcano_plot<-breast_cancer%>% ggplot(aes(x=log_MCF7,y=neglog_MCF,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 MCF7 compared to control",y="-log(p-value)")
## use ggplotly to see hover over the points to see the gene names
optional_volcano_plot<-breast_cancer%>% ggplot(aes(x=log_MCF7,y=neglog_MCF,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 MCF7 compared to control",y="-log(p-value)")
optional_volcano_plot
ggplotly(optional_volcano_plot)
EXAMPLES OF A COUPLE PROTEINS
## based on the gene symbols from plotly, select a few proteins from the table. Select just the data and gene symbols then melt
Examples_Down<-breast_cancer %>% filter(Gene_Symbol=="CTSZ" | Gene_Symbol=="APOA1" | Gene_Symbol=="HLA-A" | Gene_Symbol=="ADCK3") %>%
select(Gene_Symbol,MCF10A_1:SKBR3_2) %>%
melt()
## make barplots facetted by Gene Symbol
Example_plot_down<-Examples_Down %>% ggplot(aes(x=variable,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="sample",y="relative intensity")
Example_plot_down
## Same process for the those that are show opposite effects (decreased in a cell line or)
Examples_Up<-breast_cancer %>% filter(Gene_Symbol=="C17orf28" | Gene_Symbol=="VAV2" | Gene_Symbol=="EPPK1" ) %>%
select(Gene_Symbol,MCF10A_1:SKBR3_2) %>%
melt()
Example_plot_up<-Examples_Up %>% ggplot(aes(x=variable,y=value))+
geom_bar(stat="identity",fill="blue")+
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="sample",y="relative intensity")
Example_plot_up
Downregulated<-breast_cancer %>% filter(Gene_Symbol=="CTSZ" | Gene_Symbol=="APOA1" | Gene_Symbol=="HLA-A" | Gene_Symbol=="ADCK3") %>% select(Gene_Symbol, Description, Peptides)
Upregulated<-breast_cancer %>% filter(Gene_Symbol=="C17orf28" | Gene_Symbol=="VAV2" | Gene_Symbol=="EPPK1" )%>% select(Gene_Symbol, Description, Peptides)
Downregulated
## # A tibble: 7 x 3
## Gene_Symbol Description Peptides
## <chr> <chr> <dbl>
## 1 HLA-A HLA class I histocompatibility antigen, A-30 alpha chain… 1
## 2 HLA-A HLA class I histocompatibility antigen, A-23 alpha chain… 3
## 3 HLA-A HLA class I histocompatibility antigen, A-2 alpha chain … 21
## 4 ADCK3 Isoform 3 of Chaperone activity of bc1 complex-like, mit… 2
## 5 CTSZ Cathepsin Z OS=Homo sapiens GN=CTSZ PE=1 SV=1 8
## 6 APOA1 Apolipoprotein A-I OS=Homo sapiens GN=APOA1 PE=1 SV=1 2
## 7 HLA-A HLA class I histocompatibility antigen, A-1 alpha chain … 1
Upregulated
## # A tibble: 4 x 3
## Gene_Symbol Description Peptides
## <chr> <chr> <dbl>
## 1 C17orf28 Isoform 2 of UPF0663 transmembrane protein C17orf28 OS=H… 3
## 2 EPPK1 Uncharacterized protein OS=Homo sapiens GN=EPPK1 PE=4 SV… 4
## 3 VAV2 Isoform 2 of Guanine nucleotide exchange factor VAV2 OS=… 4
## 4 EPPK1 Epiplakin OS=Homo sapiens GN=EPPK1 PE=1 SV=2 71