#### 01_EAT_Proteomics_ForUnsupervisedAnalysis.R #####
## Access DB SQlite where Proteomics data.
## Append patient meta information
library(DBI)
library(RSQLite)
library(pheatmap)
library(plyr)
library(RColorBrewer)
library(DelayedMatrixStats)
library(randomcoloR)
library(cluster)
Deane_proteins <- read.csv("H:/Projects/COVID19/Proteomics/Files/131ProteinsofInterest_UniProtID.csv",
header = TRUE, sep = ",", stringsAsFactors = FALSE)
##### Establish a connection with DB #####
con <- dbConnect(SQLite(), dbname = "P:/All_20200428_COVID_plasma_multiomics/SQLite Database/Covid-19 Study DB.sqlite")
##### Pull data from DB #####
dbListTables(con)
# [1] "biomolecules" "deidentified_patient_metadata" "lipidomics_measurements" "lipidomics_runs"
# [5] "metabolomics_measurements" "metabolomics_runs" "metadata" "omes"
# [9] "patient_metadata" "patient_samples" "proteomics_measurements" "proteomics_runs"
# [13] "rawfiles" "sqlite_sequence"
# Extract tables from SQLite database. The square prakets provide the first row or column names of table.
dbReadTable(con, "biomolecules")[0,]
dbReadTable(con, "metadata")[1:10,]
dbReadTable(con, "metabolomics_measurements")[0,]
dbReadTable(con, "metabolomics_runs")[0,]
dbReadTable(con, "biomolecules")[0,]
dbReadTable(con, "deidentified_patient_metadata")[0,]
dbReadTable(con, "proteomics_measurements")[0,]
#[1] measurement_id replicate_id biomolecule_id raw_abundance normalized_abundance raw_ibaq normalized_ibaq
dbReadTable(con, "proteomics_runs")[0,]
#[1] replicate_id rawfile_id unique_identifier
##### Knit tables into a long dataframe #####
df <- dbGetQuery(con, "SELECT proteomics_runs.unique_identifier, proteomics_measurements.normalized_abundance, proteomics_measurements.biomolecule_id, batch
FROM proteomics_measurements
INNER JOIN proteomics_runs ON proteomics_runs.replicate_id = proteomics_measurements.replicate_id
INNER JOIN rawfiles ON rawfiles.rawfile_id = proteomics_runs.rawfile_id
INNER JOIN biomolecules on biomolecules.biomolecule_id = proteomics_measurements.biomolecule_id
WHERE rawfiles.keep = 1
AND ome_id = 1
AND biomolecules.keep = '1'
")
# Select Sample_label, Gender, ICU_1, Mech_Ventilation from deidentified_patient_metadata
deidentified_patient_meta <- dbGetQuery(con, "SELECT Sample_label, Gender, ICU_1, Mech_Ventilation
FROM deidentified_patient_metadata
")
proteins_meta <- dbGetQuery(con, "SELECT standardized_name, metadata.biomolecule_id, metadata_type, metadata_value
FROM metadata
INNER JOIN biomolecules ON biomolecules.biomolecule_id = metadata.biomolecule_id
WHERE biomolecules.omics_id = 1
AND biomolecules.keep = 1
")
dbDisconnect(con)
##### Data formatting #####
# Reshape data into wide format with first two columns containing meta data
# Column names is "normalized_abundance.X" where X is the biomolecule.id
df_wide <- reshape(df, timevar = "biomolecule_id", v.names = "normalized_abundance",
idvar = "unique_identifier", direction = "wide" )
#136 samples and 519 protein groups
# Note HC were removed
# Match rawfiles to match "20200514_ES_HC_afterBatch1" -> "20200514_ES_COON_HC_afterBatch1_single-shot"
# HC_rows <- grep("HC",df_wide$unique_identifier)
#
# for(i in 1:length(HC_rows)){
# row_df <- HC_rows[i]
# old <- df_wide$unique_identifier[row_df]
# new <- gsub('^(.{11})(.*)$', '\\1_COON\\2', old)
# new <- gsub('^(.{31})(.*)$', '\\1_single-shot\\2', new)
#
# df_wide$unique_identifier[row_df] <- new
# }
##### Meta formatting #####
## meta contains the relevant information provided in the rawfile name (aka unique_identifier)
# Subset meta data and break up string by "_"
meta <- as.data.frame(stringr::str_split_fixed(df_wide$unique_identifier,"_",7))
# Fix typo in Disease state. Update factors and eleves
meta$V4 <- sub("19","", meta$V4)
meta$V4 <- as.factor(sub("-","",meta$V4))
meta$V4 <- sub("control","pooled_plasma",meta$V4)
batch_row <- grep("single",meta$V5)
for(i in 1:length(batch_row)){
row_meta <- batch_row[i]
v5 <- meta$V5[row_meta]
v6 <- meta$V6[row_meta]
#switch selected cells
meta$V5[row_meta] <- v6
meta$V6[row_meta] <- v5
#add a number to the pooled samples
meta$V7[row_meta] <- as.integer(i)
}
meta$V5 <- sub("*.atch","", meta$V5)
meta$V5 <- as.factor(sub("after","", meta$V5))
# Pad single digit numbers with a zero ex. 2 -> '02'
samples <- vector()
for (i in 1:nrow(meta)) {
if(nchar(as.character(meta$V7[i]),type = "char")==1){
one_sample <- paste0("0",meta$V7[i])
samples <- append(samples,one_sample)
}else{
one_sample <- as.character(meta$V7[i])
samples <- append(samples,one_sample)
}
}
#overwrite sample_number padded with a zero if single digit
meta$V7 <- samples
meta$V7 <- sub("_.*","",meta$V7)
meta$V8 <- paste0(meta$V4,"_",meta$V7)
# The patient id's are used more than once.
# To make all id's unique, the make.unique function will add numbers to each duplicated entry
meta$V9 <- make.unique(as.character(meta$V8))
#add batch info column
meta <- meta[,c(1,4,5,7,8,9)]
#rename meta columns
colnames(meta) <- c("Rawfile.Date.Stamp","Disease.State", "Batch",
"Padded.Patient.Number", "Sample_Label",
"Sample_Label.unique")
##### Protein meta formattting #####
proteins_meta <- proteins_meta[which(proteins_meta$metadata_type == "gene_name"),]
# Isolate only the first protein in protein group for Majority.protein.IDs and Protein.IDs
majorityProtein <- lapply(strsplit(proteins_meta$standardized_name, ";"), '[', 1)
proteins_meta$majority.protein <- vapply(majorityProtein, paste, collapse = ", ", character(1L))
##### Prepare for heatmap ######
# Subset dataframe to include abunance values only
proteomics <- df_wide[,c(3:ncol(df_wide))]
# rename row names to Patient unique Id. (non-duplicated)
rownames(proteomics) <- meta$Sample_Label.unique
# rename columns to only the metabolite unique identifier
colnames(proteomics) <- sub(".*abundance.","",colnames(proteomics))
# Annotations for patients
patient_annotation <- meta[,c(2,3)]
patient_annotation$Disease.State <- as.factor(patient_annotation$Disease.State)
rownames(patient_annotation) <- meta$Sample_Label.unique
#patient_annotation$Date <- as.Date(patient_annotation$Date)
##### Plot Heatmap #####
my_colour = list(
Disease.State = c( `COVID` = "#D0D3CA", `NONCOVID` = "#A3A79D", `pooled_plasma` = "#191919"),
Batch = c(`1` = "#FDF0FE", `2` = "#E8C5E9", `3` = "#9FA0C3", `4` = "#8B687F", `5` = "#7B435B", `6` = "#5C3344", `7` = "#1C093D"))
scaleRYG <- colorRampPalette(c("#3C99B2","#E8CB2E","#EF2D00"), space = "rgb")(20)
pheatmap(t(proteomics),
color = scaleRYG,
annotation_colors = my_colour,
annotation_col = patient_annotation,
cluster_cols = F,
#scale = "column",
show_colnames = F,
show_rownames = F,
main = "Proteomics COVID-19 HeatMap")
#### Remove Controls ####
## Create a heatmap without Controls and append more patient metadata
##### Patient Meta formatting ####
# Remove controls from meta
#meta_patientONLY <- meta[-grep("HC",meta$Sample_Label.unique),]
meta_patientONLY <- meta[-grep("pooled",meta$Sample_Label.unique),]
# What patient meta information is missing?
PatientSample_Missing <- meta_patientONLY[which(!meta_patientONLY$Sample_Label %in% deidentified_patient_meta$Sample_label),]
if(length(which(!deidentified_patient_meta$Sample_label %in% meta_patientONLY$Sample_Label)) == 0){
Deident_patientSample_Missing <- NULL
}else{
Deident_patientSample_Missing <- deidentified_patient_meta[which(!deidentified_patient_meta$Sample_label %in% meta_patientONLY$Sample_Label),]
}
#Missing samples in Patient meta data
# Rawfile.Date.Stamp Disease.State Batch Padded.Patient.Number Sample_Label Sample_Label.unique
# 94 202005010 NONCOVID 5 26 NONCOVID_26 NONCOVID_26
# match sample labels in meta and deidentified_patient_meta == All but one match 129 samples == 128 samples
meta_patientONLY <- meta_patientONLY[which(meta_patientONLY$Sample_Label.unique %in% deidentified_patient_meta$Sample_label),]
# reorder deidentified_patient_meta data to match sample label in meta_patientONLY
#deidentified_patient_meta <- deidentified_patient_meta[-which(rownames(deidentified_patient_meta) %in% rownames(Deident_patientSample_Missing)),]
deidentified_patient_meta_ordered <- deidentified_patient_meta[order(match(deidentified_patient_meta$Sample_label,meta_patientONLY$Sample_Label)),]
# Merge deidentified_patient_meta with meta_patientONLY by sample label
meta_patientONLY_merged <- merge(meta_patientONLY,deidentified_patient_meta_ordered, by.x = "Sample_Label", by.y = "Sample_label" )
#write.csv(meta_patientONLY_merged, "P:/All_20200428_COVID_plasma_multiomics/Proteomics/EAT_unsupervised_analysis/Meta_Proteomics_PatientsONLY.csv")
##### Data without Controls #####
# Remove rows that are connected with Control Data (128 x 517)
proteomics_noControls <- proteomics[which(rownames(proteomics) %in% meta_patientONLY_merged$Sample_Label),]
# For metaboanalyst
Disease.state <- sub("_.*", "",rownames(proteomics_noControls))
metaboanalyst_df <- cbind(Disease.state,proteomics_noControls)
metaboanalyst_df <- t(metaboanalyst_df)
rownames(metaboanalyst_df) <- c("Disease.state",proteins_meta$majority.protein)
#write.csv(metaboanalyst_df, file = "H:/Projects/COVID19/Proteomics/Files/metaboanalyst_COVID_proteomics_onlyPatients_normalized.csv")
# Annotations for patients
patient_annotation <- meta_patientONLY_merged[,c(3,4,7,8,9)]
patient_annotation$Gender <- factor(patient_annotation$Gender, levels = c("F","M"))
# Substitute all 0 -> NO and 1 -> Yes in ICU_1
patient_annotation$ICU_1 <- as.character(patient_annotation$ICU_1)
patient_annotation$ICU_1 <- sub("0","No",patient_annotation$ICU_1)
patient_annotation$ICU_1 <- sub("1","Yes",patient_annotation$ICU_1)
patient_annotation$ICU_1 <- factor(patient_annotation$ICU_1, levels = c("No","Yes"))
# Substitute all 1 -> NO and 1 -> Yes in ICU_1
patient_annotation$Mech_Ventilation <- as.character(patient_annotation$Mech_Ventilation)
patient_annotation$Mech_Ventilation <- sub("0","No",patient_annotation$Mech_Ventilation)
patient_annotation$Mech_Ventilation <- sub("1","Yes",patient_annotation$Mech_Ventilation)
patient_annotation$Mech_Ventilation <- factor(patient_annotation$Mech_Ventilation, levels = c("No","Yes"))
rownames(patient_annotation) <- meta_patientONLY_merged$Sample_Label
my_colour = list(
Disease.State = c( `COVID` = "#D0D3CA", `NONCOVID` = "#A3A79D"),
Batch = c(`1` = "#FDF0FE", `2` = "#E8C5E9", `3` = "#9FA0C3", `4` = "#8B687F", `5` = "#7B435B", `6` = "#5C3344", `7` = "#1C093D"),
Gender = c(`F` = "#C0B9DD", `M` = "#234947"),
ICU_1 = c(`No` = "#D0D3CA", `Yes` = "#368BA4"),
Mech_Ventilation = c(`No` = "#D0D3CA", `Yes` = "#463F3A"))
#scaleRYG <- colorRampPalette(c("#3C99B2","#E8CB2E","#EF2D00"), space = "rgb")(20)
#scaleRYG <- colorRampPalette(c("#3C99B2","#ffffff","#EF2D00"), space = "rgb")(20)
out <- pheatmap(t(proteomics_noControls),
color = scaleRYG,
annotation_colors = my_colour,
annotation_col = patient_annotation,
#annotation_row = metabolite_annotation,
cluster_cols = T,
#scale = "column",
show_colnames = F,
show_rownames = F,
#scale = "row",
main = "Proteomics COVID-19 patients only HeatMap")
# Plot extracted gene-t0-cluster assignments
plot(out$tree_col)
#### Coagulation Cascade ####
coagulation_proteins_meta <- proteins_meta[which(proteins_meta$majority.protein %in% Deane_proteins$UniProt.ID),]
coagulation_proteomics <- proteomics_noControls[,which(colnames(proteomics_noControls) %in% coagulation_proteins_meta$biomolecule_id)]
colnames(coagulation_proteomics) <- coagulation_proteins_meta$majority.protein
coagulation_proteomics_COVID <- coagulation_proteomics[-grep("NONCOVID",rownames(coagulation_proteomics)),]
coagulation_proteomics_NONCOVID <- coagulation_proteomics[grep("NONCOVID",rownames(coagulation_proteomics)),]
out <- pheatmap(t(coagulation_proteomics),
color = scaleRYG,
annotation_colors = my_colour,
annotation_col = patient_annotation,
#annotation_row = metabolite_annotation,
cluster_cols = T,
#scale = "column",
show_colnames = F,
show_rownames = F,
#scale = "row",
main = "Coagulation Proteins of Interest\n COVID-19 patients only HeatMap")
plot(out$tree_row,
main = "Dendogram of Proteins in Subset - coagulation ")
plot(out$tree_col,
main = "Dendogram of Patients in subset - coagulation")
##### Exploratory Figures #####
proteomics_COVID <- t(proteomics_noControls)
rownames(proteomics_COVID) <- proteins_meta$majority.protein
proteomics_NONCOVID <- proteomics_COVID[,grep("NONCOVID",colnames(proteomics_COVID))]
proteomics_COVID <- proteomics_COVID[,-grep("NONCOVID",colnames(proteomics_COVID))]
# Boxplot of unknown proteomics
boxplot(proteomics_COVID,
main = "Boxplot of COVID proteomics",
ylab = "Log2(Intensity)")
par(mar = c(22.1, 4.1, 4.1, 4.1) # change the margins bottow, left, top, and right
#lwd = 2, # increase the line thickness
#cex.axis = 1.2 # increase default axis label size
)
# Boxplot of COVID proteomics
boxplot(proteomics_COVID,
main = "Boxplot of COVID proteomics",
ylab = "Log2(Intensity)",
las = 1,
xaxt = "n")
text(x = 1:length(colnames(proteomics_COVID)),
y = par("usr")[3] - 0.30, #subtract from y axis to push labels down
labels = colnames(proteomics_COVID),
xpd = NA, #print below axis
srt = 90,
adj = 1)
# > median(proteomics_COVID)
# [1] 26.30916
# > mean(proteomics_COVID)
# [1] 26.72138
par(mar = c(10, 4.1, 4.1, 4.1) # change the margins bottow, left, top, and right
#lwd = 2, # increase the line thickness
#cex.axis = 1.2 # increase default axis label size
)
# Boxplot of NONCOVID proteomics
boxplot(proteomics_NONCOVID,
main = "Boxplot of NONCOVID proteomics",
ylab = "Log2(Intensity)",
las = 1,
xaxt = "n")
text(x = 1:length(colnames(proteomics_NONCOVID)),
y = par("usr")[3] - 0.30, #subtract from y axis to push labels down
labels = colnames(proteomics_NONCOVID),
xpd = NA, #print below axis
srt = 90,
adj = 1)
# Histogram of Unknowns and Knowns
par(mar = c(6.1, 4.1, 4.1, 4.1) # change the margins bottow, left, top, and right
#lwd = 2, # increase the line thickness
#cex.axis = 1.2 # increase default axis label size
)
ICU <- patient_annotation[which(patient_annotation$ICU_1 == "Yes"),]
noICU <- patient_annotation[which(patient_annotation$ICU_1 == "No"),]
# NONCOVID ICU vs. noICU histograms
noICU_NONCOVID <- noICU[grep("NONCOVID",rownames(noICU)),]
ICU_NONCOVID <- ICU[grep("NONCOVID",rownames(ICU)),]
proteomics_NONCOVID_noICU <- proteomics_NONCOVID[,which(colnames(proteomics_NONCOVID) %in% rownames(noICU_NONCOVID))]
proteomics_NONCOVID_ICU <- proteomics_NONCOVID[,which(colnames(proteomics_NONCOVID) %in% rownames(ICU_NONCOVID))]
hist(proteomics_NONCOVID_ICU,
main = "Histogram of NONCOVID in ICU samples\n features Log2(Intensity)",
xlab = "Log2(Intensity)",
#xlim = c(0,30),
las = 1,
breaks = 20)
abline(v= median(proteomics_NONCOVID_ICU), col = "red", lwd = 2)
# median of NONCOVID patients in the ICU 26.00141
hist(proteomics_NONCOVID_noICU,
main = "Histogram of NONCOVID NOT in the ICU samples\n features Log2(Intensity)",
xlab = "Log2(Intensity)",
#xlim = c(0,30),
las = 1,
breaks = 20)
abline(v= median(proteomics_NONCOVID_noICU), col = "red", lwd = 2)
# median of NONCOVID patients in the ICU 25.95939
# COVID ICU vs noICU histograms
noICU_COVID <- noICU[-grep("NONCOVID",rownames(noICU)),]
ICU_COVID <- ICU[-grep("NONCOVID",rownames(ICU)),]
proteomics_COVID_noICU <- proteomics_COVID[,which(colnames(proteomics_COVID) %in% rownames(noICU_COVID))]
proteomics_COVID_ICU <- proteomics_COVID[,which(colnames(proteomics_COVID) %in% rownames(ICU_COVID))]
hist(proteomics_COVID_ICU,
main = "Histogram of COVID in ICU samples\n features Log2(Intensity)",
xlab = "Log2(Intensity)",
#xlim = c(0,30),
las = 1,
breaks = 20)
abline(v= median(proteomics_COVID_ICU), col = "red", lwd = 2)
# median of NONCOVID patients in the ICU 26.38699
hist(proteomics_COVID_noICU,
main = "Histogram of COVID NOT in the ICU samples\n features Log2(Intensity)",
xlab = "Log2(Intensity)",
#xlim = c(0,30),
las = 1,
breaks = 20)
abline(v= median(proteomics_COVID_noICU), col = "red", lwd = 2)
# median of NONCOVID patients in the ICU 26.22727
##### Fold Change to Median NONCOVID #####
# Calculate the median for each column = metabolite/feature
NONCOVID_Median.ProteinIntensity <- rowMedians(as.matrix(proteomics_NONCOVID))
# Foldchange function -> take every row and divide it by the median vector
fold_change.FUNC <- function(x) x-NONCOVID_Median.ProteinIntensity
df_fc <- apply(proteomics_COVID,2, fold_change.FUNC)
# New color palette patients
my_colour = list(
Gender = c(`F` = "#C0B9DD", `M` = "#234947"),
Batch = c(`1` = "#FDF0FE", `2` = "#E8C5E9", `3` = "#9FA0C3", `4` = "#8B687F", `5` = "#7B435B", `6` = "#5C3344", `7` = "#1C093D"),
ICU_1 = c(`No` = "#D0D3CA", `Yes` = "#368BA4"),
Mech_Ventilation = c(`No` = "#D0D3CA", `Yes` = "#463F3A"))
# Heatmap of the Fold Change calculated from the median in NONCOVID cohort
scaleRYG <- colorRampPalette(c("#3C99B2","#ffffff","#EF2D00"), space = "rgb")(20)
out <- pheatmap(df_fc,
color = scaleRYG,
annotation_colors = my_colour,
annotation_col = patient_annotation[-grep("NONCOVID",rownames(patient_annotation)),-c(1)],
#annotation_row = metabolite_annotation,
cluster_cols = T,
#scale = "column",
show_colnames = F,
show_rownames = F,
#scale = "row",
main = "Heatmap Fold Change COVID/NONCOVID(median)")
plot(out$tree_col,
main = "Fold Change COVID/NONCOVID(median)\n Patients")
#Export fold change data frame ordered by protein clustering
# Re-0rder original data (proteins) to match ordering in heatmap (top-to-bottom)
df_fc_export <- df_fc[rownames(df_fc[out$tree_row[["order"]],]),colnames(df_fc[,out$tree_col[["order"]]])]
#write.csv(df_fc_export, file = "P:/All_20200428_COVID_plasma_multiomics/Proteomics/EAT_unsupervised_analysis/COVID_proteomics_Heatmap_COVIDrelativeNONCOVIDmedian.csv")
# Histogram of Fold Changes
par(mar = c(6.1, 4.1, 4.1, 4.1))
hist(df_fc,
main = "Histogram of Fold Changes COVID relative to NONCOVID",
xlab = "Fold Change in Log2 space",
xlim = c(-15,15),
breaks = 50)
# Protein groups with fold change greater than 1.2
table(rowSums(abs(df_fc)>1.2)>0)
# Protein groups with fold change greater than 2
table(rowSums(abs(df_fc)>2)>0)
# Protein groups with fold change greater than 4
table(rowSums(abs(df_fc)>4)>0)
# Protein groups with fold change greater than 8
table(rowSums(abs(df_fc)>8)>0)
# subset features with fold change greater than 4.2 to cluster
important_features_foldchange <- df_fc[rowSums(abs(df_fc)>6.2)>0,]
out <- pheatmap(important_features_foldchange,
color = scaleRYG,
annotation_colors = my_colour,
annotation_col = patient_annotation[-grep("NONCOVID",rownames(patient_annotation)),-c(1)],
#annotation_row = metabolite_annotation,
cluster_cols = T,
#scale = "column",
show_colnames = F,
show_rownames = F,
#scale = "row",
main = "Heatmap Fold Change COVID/NONCOVID(median)\n Subset Features that have abs(FC) > 6.2 in at least one patient")
##### Fold Change Median of NONCOVID not in ICU ####
# Calculate the median for each column = metabolite/feature
NONCOVID_NOICU_Median.ProteinIntensity <- rowMedians(as.matrix(proteomics_NONCOVID_noICU))
# Foldchange function -> take every row and divide it by the median vector
fold_change.FUNC <- function(x) x-NONCOVID_NOICU_Median.ProteinIntensity
df_fc <- apply(proteomics_COVID,2, fold_change.FUNC)
# Heatmap of the Fold Change calculated from the median in NONCOVID cohort
scaleRYG <- colorRampPalette(c("#3C99B2","#ffffff","#EF2D00"), space = "rgb")(20)
pheatmap(df_fc,
color = scaleRYG,
annotation_colors = my_colour,
annotation_col = patient_annotation[-grep("NONCOVID",rownames(patient_annotation)),-c(1)],
#annotation_row = metabolite_annotation,
cluster_cols = T,
#scale = "column",
show_colnames = F,
show_rownames = F,
#scale = "row",
main = "Heatmap Fold Change COVID/NONCOVID not in the ICU(median)")
##### Fold Change Median of NONCOVID not in ICU ####
# Calculate the median for each column = metabolite/feature
NONCOVID_NOICU_Median <- rowMedians(proteomics_NONCOVID_noICU)
# Foldchange function -> take every row and divide it by the median vector
fold_change.FUNC <- function(x) x-NONCOVID_NOICU_Median
df_fc <- apply(proteomics_COVID,2, fold_change.FUNC)
# Heatmap of the Fold Change calculated from the median in NONCOVID cohort
scaleRYG <- colorRampPalette(c("#3C99B2","#ffffff","#EF2D00"), space = "rgb")(20)
pheatmap(df_fc,
color = scaleRYG,
annotation_colors = my_colour,
annotation_col = patient_annotation[-grep("NONCOVID",rownames(patient_annotation)),-c(1)],
#annotation_row = metabolite_annotation,
cluster_cols = T,
#scale = "column",
show_colnames = F,
show_rownames = F,
#scale = "row",
main = "Heatmap Fold Change COVID/NONCOVID not in the ICU(median)")
##### Fold Change in Coagulation Subset ######
# Calculate the median for each column = metabolite/feature
NONCOVID_Coagulation_Median.ProteinIntensity <- colMedians(as.matrix(coagulation_proteomics_NONCOVID))
# Foldchange function -> take every row and divide it by the median vector
fold_change.FUNC <- function(x) x-NONCOVID_Coagulation_Median.ProteinIntensity
df_fc <- apply(coagulation_proteomics_COVID,1, fold_change.FUNC)
# Heatmap of the Fold Change calculated from the median in NONCOVID cohort
scaleRYG <- colorRampPalette(c("#3C99B2","#ffffff","#EF2D00"), space = "rgb")(20)
out <- pheatmap(df_fc,
color = scaleRYG,
annotation_colors = my_colour,
annotation_col = patient_annotation[-grep("NONCOVID",rownames(patient_annotation)),-c(1)],
#annotation_row = metabolite_annotation,
cluster_cols = T,
#scale = "column",
show_colnames = F,
show_rownames = F,
#scale = "row",
main = "Heatmap Fold Change COVID/NONCOVID(median) in Coagulation Subset Proteins")
plot(out$tree_row,
main = "Dendogram of Proteins\n Fold Change COVID/NONCOVID(median) in Coagulation Subset")
plot(out$tree_col,
main = "Dendogram of COVID Patients\n Fold Change COVID/NONCOVID(median) in Coagulation Subset")
Datasets
Models
