Alex Anuta - Breast Cancer Cell Analysis
Working from the COVID.rmd model, fill in the blank chunks with code to analize the dataset of your choice. You will need to think about what you are copying and pasting and what needs to be updated to make the code relevant to your chosen dataset. Double check all syntax and be precise!
###First, let’s 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
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)
head(BC_mat1)
## log_MCF7 log_MDA231 log_MDA468 log_SKBR3
## [1,] -0.43 1.90 -0.74 0
## [2,] -0.96 -0.35 -0.71 0
## [3,] 0.27 0.29 -0.04 0
## [4,] 0.88 0.78 -0.05 0
## [5,] 0.79 1.19 0.92 0
## [6,] -0.79 -0.90 -1.28 0## 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 of interest
BC_data<- BC_data %>% mutate(neglog_MCF7=-log10(pvalue_MCF7_vs_MCF10A))
## Use ggplot to plot the log ratio a column (cell line or treatment) of interest against the -log p-value of that same column of interest
volcano_plot<-BC_data %>% ggplot(aes(x=log_MCF7,y=neglog_MCF7,description=Gene_Symbol))+
geom_point(alpha=0.7,color="royalblue")
volcano_plot
## BETTER VOLCANO PLOT
## Define the significant ones so they can be colored
BC_data<-BC_data %>% mutate(significance=ifelse((log_MCF7>2 & neglog_MCF7>2.99),"UP", ifelse((log_MCF7<c(-2) & neglog_MCF7>2.99),"DOWN","NOT SIG")))
## Some standard colors
plain_cols3<-c("red","gray","blue")
## volcano plot as before with some added graph features
better_volcano_plot<-BC_data %>% ggplot(aes(x=log_MCF7,y=neglog_MCF7,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 of MCF7 line compared to control (MCF10A)",y="-log(MCF7)")
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
BC_long<-pivot_longer(BC_data, cols = c(MCF10A_1:SKBR3_2), names_to = 'variable')%>% select(-c(pvalue_MCF7_vs_MCF10A:significance))%>%select(-Description)
# Make a pivot longer table of the data for initial columns. Exclude the p-values and log columns
head(BC_long)
## # A tibble: 6 × 4
## Gene_Symbol Peptides variable value
## <chr> <dbl> <chr> <dbl>
## 1 NES 7 MCF10A_1 9.54
## 2 NES 7 MCF10A_2 4.58
## 3 NES 7 MCF7_1 5.07
## 4 NES 7 MCF7_2 5.42
## 5 NES 7 MDA231_1 25.4
## 6 NES 7 MDA231_2 27.4
## 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=="GFRA1" | Gene_Symbol=="HLA-C" | Gene_Symbol=="SOX9")
head(Examples_Down)
## # A tibble: 6 × 4
## Gene_Symbol Peptides variable value
## <chr> <dbl> <chr> <dbl>
## 1 GFRA1 1 MCF10A_1 0.001
## 2 GFRA1 1 MCF10A_2 2.58
## 3 GFRA1 1 MCF7_1 46.1
## 4 GFRA1 1 MCF7_2 44.0
## 5 GFRA1 1 MDA231_1 1.56
## 6 GFRA1 1 MDA231_2 1.9
## make barplots facetted by Gene Symbol (if copying COVID example - change x=order to x = variable)
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 upregulated ones
Examples_Up<-BC_long %>% filter(Gene_Symbol=="GFRA1" | Gene_Symbol=="DST")
head(Examples_Up)
## # A tibble: 6 × 4
## Gene_Symbol Peptides variable value
## <chr> <dbl> <chr> <dbl>
## 1 GFRA1 1 MCF10A_1 0.001
## 2 GFRA1 1 MCF10A_2 2.58
## 3 GFRA1 1 MCF7_1 46.1
## 4 GFRA1 1 MCF7_2 44.0
## 5 GFRA1 1 MDA231_1 1.56
## 6 GFRA1 1 MDA231_2 1.9
Example_plot_up<-Examples_Up %>% ggplot(aes(x=variable,y=value))+
geom_bar(stat="identity",fill="royalblue")+
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
BONUS CODE- for figuring out what the gene symbols are
#your original table contains descriptions of each gene. Write some code to create yourself a table to filter for your gene symbols of interest and select for the gene symbol, description, and peptides. Do this for both the significantly upregulated and downregulated points of interest.
#Downregulated
downregulated <- BC_data %>%
filter(significance == 'DOWN') %>%
select(Gene_Symbol, Description, Peptides)
downregulated
## # A tibble: 12 × 3
## Gene_Symbol Description Pepti…¹
## <chr> <chr> <dbl>
## 1 FMNL1 Isoform 2 of Formin-like protein 1 OS=Homo sapiens GN=FM… 1
## 2 HLA-B HLA class I histocompatibility antigen, B-52 alpha chain… 1
## 3 CTSZ Cathepsin Z OS=Homo sapiens GN=CTSZ PE=1 SV=1 8
## 4 PLAU Isoform 2 of Urokinase-type plasminogen activator OS=Hom… 2
## 5 HIST1H1A Histone H1.1 OS=Homo sapiens GN=HIST1H1A PE=1 SV=3 1
## 6 TMBIM1 Transmembrane BAX inhibitor motif-containing protein 1 O… 1
## 7 ADCK3 Isoform 3 of Chaperone activity of bc1 complex-like, mit… 2
## 8 METTL7A Methyltransferase-like protein 7A OS=Homo sapiens GN=MET… 3
## 9 EXOG Nuclease EXOG, mitochondrial OS=Homo sapiens GN=EXOG PE=… 3
## 10 PKP2 Isoform 1 of Plakophilin-2 OS=Homo sapiens GN=PKP2 7
## 11 HLA-A HLA class I histocompatibility antigen, A-1 alpha chain … 1
## 12 APOA1 Apolipoprotein A-I OS=Homo sapiens GN=APOA1 PE=1 SV=1 2
## # … with abbreviated variable name ¹Peptides
#Upregulated
upregulated <- BC_data %>%
filter(significance == 'UP') %>%
select(Gene_Symbol, Description, Peptides)
upregulated
## # A tibble: 57 × 3
## Gene_Symbol Description Pepti…¹
## <chr> <chr> <dbl>
## 1 MUC5B Mucin-5B OS=Homo sapiens GN=MUC5B PE=1 SV=3 2
## 2 AGR3 Anterior gradient protein 3 homolog OS=Homo sapiens GN=A… 5
## 3 USP32 Ubiquitin carboxyl-terminal hydrolase 32 OS=Homo sapiens… 2
## 4 AARS Alanyl-tRNA synthetase, cytoplasmic OS=Homo sapiens GN=A… 23
## 5 EPB41L5 Isoform 2 of Band 4.1-like protein 5 OS=Homo sapiens GN=… 5
## 6 FBP1 Fructose-1,6-bisphosphatase 1 OS=Homo sapiens GN=FBP1 PE… 8
## 7 ISOC1 Isochorismatase domain-containing protein 1 OS=Homo sapi… 6
## 8 ELP3 Isoform 2 of Elongator complex protein 3 OS=Homo sapiens… 2
## 9 SYTL2 Isoform 8 of Synaptotagmin-like protein 2 OS=Homo sapien… 1
## 10 EML1 Isoform 3 of Echinoderm microtubule-associated protein-l… 2
## # … with 47 more rows, and abbreviated variable name ¹Peptides