data_preparation.Rmd

---
title: "cancer_blood_biomarkers"
author: "Hansen Han"
date: "2024-06-18"
output: html_document
---

# Cancer Blood Biomarker Discovery (Data Preparation)
Analysis to identify diagnostic biomarkers of cancer via meta-analysis of gene expression studies.

## Load Libraries
```{r messages=FALSE errors=FALSE}
library(MetaIntegrator)
library(tidyverse)
library(magrittr)
library(ggpubr)
library(biomaRt)
```

## Prepare Data for Analysis

### Download From GEO 

Download Discovery & Validation Study Data
```{r}
study_ids <- c(
  "GSE198048",
  "GSE138118",
  "GSE136651",
  "GSE74629",
  "GSE125158",
  "GSE20189",
  "GSE16443",
  "GSE39345",
  "GSE120691",
  "GSE164191",
  "GSE49515",
  "GSE27562",
  "GSE10715"
)
```

```{r}
# ran into an issue here downloading data with getGEOData(), fixed by installing GEOquery.  
blood_biomarker_data_raw <- getGEOData(gseVector = study_ids, destdir="/Users/hansenhan/cancer_blood_biomarkers")
blood_biomarker_data_raw %>% save(file = "blood_biomarker_data_raw.RData")
```

Identify studies missing expression data
```{r}
studies_missing_expr_data <- c()
for (study in blood_biomarker_data_raw$originalData) {
  if (study$exp_comment == "Expression data is missing") {
    studies_missing_expr_data <- c(studies_missing_expr_data, study$name)
  }
}

studies_missing_expr_data
```

### Load Manual Data

#### GSE198048
(This is log2 format)
```{r}
GSE198048_data <- read_tsv("GSE198048_Normalized_mRNA_matrix.txt")
GSE198048_data
```
The data is loaded incorrectly, with the columns needing to be shifted to the right. 

```{r}
sample_ids <- colnames(GSE198048_data)
sample_ids
```
```{r}
# set probe ids to index
GSE198048_data %<>% column_to_rownames(var = "100")
GSE198048_data
```


Last column is loaded incorrectly, ie: "0.795771355682911\t3.45972844016865", need to split it out and add on a new last column

```{r}
corrected_last_column_data <- c()
corrected_new_last_column_data <- c()

for (entry in GSE198048_data$`99`) {
  split_data <- strsplit(entry, "\t")
  corrected_last_column_data <- c(corrected_last_column_data, split_data[[1]][1])
  corrected_new_last_column_data <- c(corrected_new_last_column_data, split_data[[1]][2])
}

GSE198048_data$`99` <- corrected_last_column_data
GSE198048_data$newcolumn <- corrected_new_last_column_data
GSE198048_data
```

Rename the data with the correct sample names
```{r}
colnames(GSE198048_data) <- sample_ids
GSE198048_data
```

Assign GEO accession numbers to expression samples
```{r}

new_sample_ids <- c()

# rotate so we can pair rename samples easily
expr_data <- as.data.frame(t(GSE198048_data)) %>% rownames_to_column(var = "patient_id")


# find the matching GSM sample id for each patient
for (patient_id_expr in expr_data$patient_id) {
  for (title_str in blood_biomarker_data_raw$originalData$GSE198048$pheno %>% pull(title)) {
    patient_id_meta <- strsplit(title_str, "_")[[1]][2]
    if (patient_id_expr == patient_id_meta) {
      gsm_id <- blood_biomarker_data_raw$originalData$GSE198048$pheno %>% filter(title == title_str) %>% pull(geo_accession)
      new_sample_ids <- c(new_sample_ids, gsm_id)
    }
  }
    
}

# save them
expr_data$patient_id <- new_sample_ids

# rotate it back
expr_data %<>% column_to_rownames(var = "patient_id") %>% t() %>% as.data.frame()
expr_data
```
Save to MetaIntegrator object
```{r}
blood_biomarker_data_raw$originalData$GSE198048$expr <- expr_data
blood_biomarker_data_raw$originalData$GSE198048$exp_comment <- "Expression data manually loaded"
```

```{r}
blood_biomarker_data_raw$originalData$GSE198048$pheno
```


#### GSE136651
We have GPL18573 - Illumina NextSeq 5000
We also have GPL24106 - MinION (ignore)

Remove platform-specific study data, make it so that it is only Illumina platform....there are only 2 nanopore samples
```{r}
blood_biomarker_data_raw$originalData$GSE136651 <- blood_biomarker_data_raw$originalData$GSE136651_GPL18573
blood_biomarker_data_raw$originalData$GSE136651_GPL18573 <- NULL
blood_biomarker_data_raw$originalData$GSE136651_GPL24106 <- NULL
```

Load TPMs
Originally tried to load "GSE136651_norm_counts_TPM_GRCh38.p13_NCBI.tsv.gz" from NCBI, but it only has 16 samples when there should be around 30+ (must have gotten snagged on the nanostrining samples embedded in there...)

Instead, given the BioProject ID (PRJNA562932), I used my ```ena_rnaseq_pipeline pipeline``` (https://github.com/hansenrhan/ena_rnaseq_quantification) to download the FASTQs, quantify them with Salmon and combine them into a single csv.


```{r}
GSE136651_expr <- read_csv("PRJNA562932_TPMs.csv")
GSE136651_expr$...1 <- NULL
GSE136651_expr
```

```{r}
blood_biomarker_data_raw$originalData$GSE136651$expr <- GSE136651_expr
blood_biomarker_data_raw$originalData$GSE136651$exp_comment <- "Expression data manually loaded"
```



#### GSE120691

Decided to also use my rnaseq quant pipeline instead of the file from GEO (GSE120691_norm_counts_TPM_GRCh38.p13_NCBI.tsv.gz)

```{r}
GSE120691_expr <- read_csv("PRJNA494155_TPMs.csv")
GSE120691_expr
```

```{r}
blood_biomarker_data_raw$originalData$GSE120691$expr <- GSE120691_expr
blood_biomarker_data_raw$originalData$GSE120691$exp_comment <- "Expression data manually loaded"
```


```{r}
blood_biomarker_data_raw$originalData$GSE120691$pheno
```


### Label Study Metadata
Here we'll use 1 to denote cancer and 0 to denote healthy control

#### GSE136651
```{r}
pheno_data <- blood_biomarker_data_raw$originalData$GSE136651$pheno
for (sample_id in pheno_data$geo_accession) {
  variable_to_examine <- pheno_data %>% filter(geo_accession == sample_id) %>% pull(`subject status:ch1`)
  if (variable_to_examine == "normal healthy donor") {
    label_value <- 0
  } else if (variable_to_examine == "pancreatic cancer patient") {
    label_value <- 1
  } else {
    label_value <- 2
  }
  
  # assign the new label value
  blood_biomarker_data_raw$originalData$GSE136651$class[sample_id] <- label_value
}

# remove filtered smaples from labels
blood_biomarker_data_raw$originalData$GSE136651$class <- blood_biomarker_data_raw$originalData$GSE136651$class[blood_biomarker_data_raw$originalData$GSE136651$class != 2]

# remove case excluded from metadata
blood_biomarker_data_raw$originalData$GSE136651$pheno %<>% filter(geo_accession %in% names(blood_biomarker_data_raw$originalData$GSE136651$class))

#Get a count of how many healthy and diseased samples there are
table(blood_biomarker_data_raw$originalData$GSE136651$class)
```



#### GSE120691
```{r}
pheno_data <- blood_biomarker_data_raw$originalData$GSE120691$pheno
for (sample_id in pheno_data$geo_accession) {
  variable_to_examine <- pheno_data %>% filter(geo_accession == sample_id) %>% pull(`disease state:ch1`)
  if (variable_to_examine == "healthy control") {
    label_value <- 0
  } else if (variable_to_examine == "cervical cancer") {
    label_value <- 1
  } else {
    label_value <- 2
  }
  
  # assign the new label value
  blood_biomarker_data_raw$originalData$GSE120691$class[sample_id] <- label_value
}

# remove filtered smaples from labels
blood_biomarker_data_raw$originalData$GSE120691$class <- blood_biomarker_data_raw$originalData$GSE120691$class[blood_biomarker_data_raw$originalData$GSE120691$class != 2]

# remove case excluded from metadata
blood_biomarker_data_raw$originalData$GSE120691$pheno %<>% filter(geo_accession %in% names(blood_biomarker_data_raw$originalData$GSE120691$class))

#Get a count of how many healthy and diseased samples there are
table(blood_biomarker_data_raw$originalData$GSE120691$class)
```



#### GSE198048
```{r}
pheno_data <- blood_biomarker_data_raw$originalData$GSE198048$pheno
for (sample_id in pheno_data$geo_accession) {
  if (pheno_data %>% filter(geo_accession == sample_id) %>% pull(`histology:ch1`) == "NA") {
    # if it is NA, it is false positive and a healthy control
    label_value <- 0
  } else {
    label_value <- 1
  }
  
  # assign the new label value
  blood_biomarker_data_raw$originalData$GSE198048$class[sample_id] <- label_value
}

#Get a count of how many healthy and diseased samples there are
table(blood_biomarker_data_raw$originalData$GSE198048$class)
```


#### GSE138118
```{r}
pheno_data <- blood_biomarker_data_raw$originalData$GSE138118$pheno
for (sample_id in pheno_data$geo_accession) {
  if (pheno_data %>% filter(geo_accession == sample_id) %>% pull(`description`) == "Health volunteer Blood") {
    # if it is NA, it is false positive and a healthy control
    label_value <- 0
  } else {
    label_value <- 1
  }
  
  # assign the new label value
  blood_biomarker_data_raw$originalData$GSE138118$class[sample_id] <- label_value
}

#Get a count of how many healthy and diseased samples there are
table(blood_biomarker_data_raw$originalData$GSE138118$class)
```


#### GSE74629
```{r}
pheno_data <- blood_biomarker_data_raw$originalData$GSE74629$pheno
for (sample_id in pheno_data$geo_accession) {
  if (pheno_data %>% filter(geo_accession == sample_id) %>% pull(`diagnosis:ch1`) == "healthy") {
    # if it is NA, it is false positive and a healthy control
    label_value <- 0
  } else {
    label_value <- 1
  }
  
  # assign the new label value
  blood_biomarker_data_raw$originalData$GSE74629$class[sample_id] <- label_value
}

#Get a count of how many healthy and diseased samples there are
table(blood_biomarker_data_raw$originalData$GSE74629$class)
```

#### GSE125158
```{r}
pheno_data <- blood_biomarker_data_raw$originalData$GSE125158$pheno
for (sample_id in pheno_data$geo_accession) {
  if (pheno_data %>% filter(geo_accession == sample_id) %>% pull(`diagnosis:ch1`) == "healthy") {
    # if it is NA, it is false positive and a healthy control
    label_value <- 0
  } else {
    label_value <- 1
  }
  
  # assign the new label value
  blood_biomarker_data_raw$originalData$GSE125158$class[sample_id] <- label_value
}

#Get a count of how many healthy and diseased samples there are
table(blood_biomarker_data_raw$originalData$GSE125158$class)
```


#### GSE20189
Case excluded should be dropped.
```{r}
pheno_data <- blood_biomarker_data_raw$originalData$GSE20189$pheno
for (sample_id in pheno_data$geo_accession) {
  variable_to_examine <- pheno_data %>% filter(geo_accession == sample_id) %>% pull(`description`)
  if (variable_to_examine == "Control") {
    label_value <- 0
  } else if (variable_to_examine == "Case") {
    label_value <- 1
  } else {
    label_value <- 2
  }
  
  # assign the new label value
  blood_biomarker_data_raw$originalData$GSE20189$class[sample_id] <- label_value
}

# remove from labels
blood_biomarker_data_raw$originalData$GSE20189$class <- blood_biomarker_data_raw$originalData$GSE20189$class[blood_biomarker_data_raw$originalData$GSE20189$class != 2]

# remove samples from metadata 
blood_biomarker_data_raw$originalData$GSE20189$pheno %<>% filter(geo_accession %in% names(blood_biomarker_data_raw$originalData$GSE20189$class))


#Get a count of how many healthy and diseased samples there are
table(blood_biomarker_data_raw$originalData$GSE20189$class)
```

#### GSE16443
```{r}
pheno_data <- blood_biomarker_data_raw$originalData$GSE16443$pheno
for (sample_id in pheno_data$geo_accession) {
  if (pheno_data %>% filter(geo_accession == sample_id) %>% pull(`status:ch1`) == "Healthy") {
    # if it is NA, it is false positive and a healthy control
    label_value <- 0
  } else {
    label_value <- 1
  }
  
  # assign the new label value
  blood_biomarker_data_raw$originalData$GSE16443$class[sample_id] <- label_value
}

#Get a count of how many healthy and diseased samples there are
table(blood_biomarker_data_raw$originalData$GSE16443$class)
```

#### GSE39345
```{r}
pheno_data <- blood_biomarker_data_raw$originalData$GSE39345$pheno
for (sample_id in pheno_data$geo_accession) {
  variable_to_examine <- pheno_data %>% filter(geo_accession == sample_id) %>% pull(`sample type:ch1`)
  if (variable_to_examine == "healthy") {
    label_value <- 0
  } else if (variable_to_examine == "before chemotherapy") {
    label_value <- 1
  } else {
    label_value <- 2
  }
  
  # assign the new label value
  blood_biomarker_data_raw$originalData$GSE39345$class[sample_id] <- label_value
}

# remove case excluded from metadata
blood_biomarker_data_raw$originalData$GSE39345$pheno %<>% filter(`sample type:ch1` != "after chemotherapy")

# remove from labels
blood_biomarker_data_raw$originalData$GSE39345$class <- blood_biomarker_data_raw$originalData$GSE39345$class[blood_biomarker_data_raw$originalData$GSE39345$class != 2]

#Get a count of how many healthy and diseased samples there are
table(blood_biomarker_data_raw$originalData$GSE39345$class)
```

#### GSE164191
```{r}
pheno_data <- blood_biomarker_data_raw$originalData$GSE164191$pheno
for (sample_id in pheno_data$geo_accession) {
  if (pheno_data %>% filter(geo_accession == sample_id) %>% pull(`disease status:ch1`) == "normal") {
    # if it is NA, it is false positive and a healthy control
    label_value <- 0
  } else {
    label_value <- 1
  }
  
  # assign the new label value
  blood_biomarker_data_raw$originalData$GSE164191$class[sample_id] <- label_value
}

#Get a count of how many healthy and diseased samples there are
table(blood_biomarker_data_raw$originalData$GSE164191$class)
```

#### GSE49515
```{r}
pheno_data <- blood_biomarker_data_raw$originalData$GSE49515$pheno
for (sample_id in pheno_data$geo_accession) {
  if (pheno_data %>% filter(geo_accession == sample_id) %>% pull(`disease status:ch1`) == "healthy") {
    # if it is NA, it is false positive and a healthy control
    label_value <- 0
  } else {
    label_value <- 1
  }
  
  # assign the new label value
  blood_biomarker_data_raw$originalData$GSE49515$class[sample_id] <- label_value
}

#Get a count of how many healthy and diseased samples there are
table(blood_biomarker_data_raw$originalData$GSE49515$class)
```




#### GSE27562
```{r}
pheno_data <- blood_biomarker_data_raw$originalData$GSE27562$pheno
for (sample_id in pheno_data$geo_accession) {
  variable_to_examine <- pheno_data %>% filter(geo_accession == sample_id) %>% pull(`phenotype:ch1`)
  if (variable_to_examine == "Normal") {
    label_value <- 0
  } else if (variable_to_examine %in% c("Malignant", "Ectopic", "Pre-Surgery (aka Malignant)")) {
    label_value <- 1
  } else {
    # remove Post-Surgery" and "Bennign" samples
    label_value <- 2
  }
  
  # remove non PBMCs
  variable_to_examine <- pheno_data %>% filter(geo_accession == sample_id) %>% pull(`tissue:ch1`)
  if (variable_to_examine == "peripheral blood leukocytes") {
    label_value <- 2
  }
  
  
  

  # assign the new label value
  blood_biomarker_data_raw$originalData$GSE27562$class[sample_id] <- label_value
}

# remove filtered smaples from labels
blood_biomarker_data_raw$originalData$GSE27562$class <- blood_biomarker_data_raw$originalData$GSE27562$class[blood_biomarker_data_raw$originalData$GSE27562$class != 2]

# remove case excluded from metadata
blood_biomarker_data_raw$originalData$GSE27562$pheno %<>% filter(geo_accession %in% names(blood_biomarker_data_raw$originalData$GSE27562$class))

#Get a count of how many healthy and diseased samples there are
table(blood_biomarker_data_raw$originalData$GSE27562$class)
```


#### GSE10715
```{r}
pheno_data <- blood_biomarker_data_raw$originalData$GSE10715$pheno
for (sample_id in pheno_data$geo_accession) {
  if (pheno_data %>% filter(geo_accession == sample_id) %>% pull(`Blood:ch1`) == "Normal") {
    # if it is NA, it is false positive and a healthy control
    label_value <- 0
  } else {
    label_value <- 1
  }
  
  # assign the new label value
  blood_biomarker_data_raw$originalData$GSE10715$class[sample_id] <- label_value
}

#Get a count of how many healthy and diseased samples there are
table(blood_biomarker_data_raw$originalData$GSE10715$class)
```





### Convert Expression to Gene Names & Filter Expression Data
- Use conversion tables downloaded from GEO in biomarkers/platform_conversion_tables/ and bioMart to convert gene names
- Also remove any samples that are not found in the metadata

```{r}
library(biomaRt)
ensembl <- useEnsembl(biomart = "ensembl", 
                   dataset = "hsapiens_gene_ensembl", 
                   mirror = "useast")
attributes = listAttributes(ensembl)
attributes[1:5,]
```

```{r}
attributes %>% filter(grepl("abi", name))
```

```{r}
agilent_gene_map <- getBM(attributes = c('agilent_wholegenome_4x44k_v2', 'entrezgene_id', 'ensembl_transcript_id' , "ensembl_transcript_id_version", 'ensembl_gene_id', 'external_gene_name'), mart = ensembl)
agilent_gene_map
```


```{r}
illumina_gene_map <- getBM(attributes = c('illumina_humanref_8_v3', 'illumina_humanwg_6_v3', 'entrezgene_id', 'ensembl_transcript_id' , "ensembl_transcript_id_version", 'ensembl_gene_id', 'external_gene_name'), mart = ensembl)
illumina_gene_map
```


```{r}
ensembl_entrezgene_ensembl_gene_map <- getBM(attributes = c('affy_hg_u133_plus_2', 'entrezgene_id', 'ensembl_transcript_id' , "ensembl_transcript_id_version", 'ensembl_gene_id', 'external_gene_name'), mart = ensembl)
ensembl_entrezgene_ensembl_gene_map
```

#### Identify Metadata/Expression Mismatches
```{r}
samples_in_metadata_not_in_expr <- function(dataset_object) {
  missing_samples <- c()
  studies_with_issues <- c()

  sample_names <- dataset_object$pheno %>% pull(geo_accession) %>% unique()
  for (sample_name in sample_names) {
    if (!sample_name %in% colnames(data.frame(dataset_object$expr))) {
      missing_samples <- c(missing_samples, sample_name)
      studies_with_issues <- c(studies_with_issues, dataset_object$formattedName)
    }
  }
  print("Samples w/ metadata and missing expression data")
  print(missing_samples %>% unique())
}

samples_in_expr_not_in_metadata <- function(dataset_object) {
  missing_samples <- c()
  studies_with_issues <- c()
  sample_names <- colnames(data.frame(dataset_object$expr))
  for (sample_name in sample_names) {
    if (!sample_name %in% (dataset_object$pheno %>% pull(geo_accession) %>% unique())) {
      missing_samples <- c(missing_samples, sample_name)
      studies_with_issues <- c(studies_with_issues, dataset_object$formattedName)
    }
  }
  print("Samples w/ expression data and missing metadata")
  print(missing_samples %>% unique()) 
  
}



for (dataset_obj in blood_biomarker_data_raw$originalData) {
  print(dataset_obj$formattedName)
  samples_in_metadata_not_in_expr(dataset_obj)
  samples_in_expr_not_in_metadata(dataset_obj)
}

```

Note to self: 
The samples you see above in GSE120691 need to be removed (they are neoplasmia not cancer...)



#### GSE136651
```{r}
expr_data <- blood_biomarker_data_raw$originalData$GSE136651$expr
expr_data$ensembl_transcript_id_version <- expr_data$Name
expr_data$Name <- NULL
expr_data %<>% left_join(ensembl_entrezgene_ensembl_gene_map %>% dplyr::select(ensembl_transcript_id_version, external_gene_name), by=c("ensembl_transcript_id_version"="ensembl_transcript_id_version")) %>% filter(!is.na(external_gene_name)) %>% dplyr::select(-ensembl_transcript_id_version) %>%
  mutate(across(-external_gene_name, as.numeric)) %>% group_by(external_gene_name) %>% summarise(across(where(is.numeric), mean, na.rm = TRUE)) %>% ungroup() %>% filter(external_gene_name != "") %>% column_to_rownames(var="external_gene_name") 

# remove any samples that have been filtered out of the metadata
samples_in_metadata <- blood_biomarker_data_raw$originalData$GSE136651$pheno %>% pull(geo_accession)
expr_data %<>% select(!!!syms(samples_in_metadata))


# set keys
blood_biomarker_data_raw$originalData$GSE136651$keys <- rownames(expr_data)
blood_biomarker_data_raw$originalData$GSE136651$key_comment <- "Keys manually loaded"

# log transform the data since it is in TPM format 
blood_biomarker_data_raw$originalData$GSE136651$expr <- log10(as.matrix(expr_data) + 1)
```




#### GSE120691
```{r}
expr_data <- blood_biomarker_data_raw$originalData$GSE120691$expr
expr_data$ensembl_transcript_id_version <- expr_data$Name
expr_data$Name <- NULL
expr_data$`...1` <- NULL
expr_data %<>% left_join(ensembl_entrezgene_ensembl_gene_map %>% dplyr::select(ensembl_transcript_id_version, external_gene_name), by=c("ensembl_transcript_id_version"="ensembl_transcript_id_version")) %>% filter(!is.na(external_gene_name)) %>% dplyr::select(-ensembl_transcript_id_version) %>%
  mutate(across(-external_gene_name, as.numeric)) %>% group_by(external_gene_name) %>% summarise(across(where(is.numeric), mean, na.rm = TRUE)) %>% ungroup() %>% filter(external_gene_name != "") %>% column_to_rownames(var="external_gene_name") 

# remove any samples that have been filtered out of the metadata
samples_in_metadata <- blood_biomarker_data_raw$originalData$GSE120691$pheno %>% pull(geo_accession)
expr_data %<>% select(!!!syms(samples_in_metadata))

# set keys
blood_biomarker_data_raw$originalData$GSE120691$keys <- rownames(expr_data)
blood_biomarker_data_raw$originalData$GSE120691$key_comment <- "Keys manually loaded"

# log transform the data since it is in TPM format
blood_biomarker_data_raw$originalData$GSE120691$expr <- log2(as.matrix(expr_data) + 1)
```



#### GSE198048
```{r}
expr_data <- blood_biomarker_data_raw$originalData$GSE198048$expr
expr_data %<>% rownames_to_column(var="ensembl_gene_id") %>% left_join(ensembl_entrezgene_ensembl_gene_map %>% dplyr::select(ensembl_gene_id, external_gene_name), by=c("ensembl_gene_id"="ensembl_gene_id")) %>% filter(!is.na(external_gene_name)) %>% dplyr::select(-ensembl_gene_id) %>%
  mutate(across(-external_gene_name, as.numeric)) %>% group_by(external_gene_name) %>% summarise(across(where(is.numeric), mean, na.rm = TRUE)) %>% ungroup() %>% filter(external_gene_name != "") %>% column_to_rownames(var="external_gene_name") 



# remove any samples that have been filtered out of the metadata
samples_in_metadata <- blood_biomarker_data_raw$originalData$GSE198048$pheno %>% pull(geo_accession)
expr_data %<>% select(!!!syms(samples_in_metadata))

# set keys
blood_biomarker_data_raw$originalData$GSE198048$keys <- rownames(expr_data)
blood_biomarker_data_raw$originalData$GSE198048$key_comment <- "Keys manually loaded"

# data is already log2 transformed according to GEO
blood_biomarker_data_raw$originalData$GSE198048$expr <- as.matrix(expr_data)
```


#### GSE138118

Load the conversion table for the platform (GPL17692	[HuGene-2_1-st] Affymetrix Human Gene 2.1 ST Array [transcript (gene) version])
```{r}
GPL17692_conversion_table <- read_tsv("/Users/hansenhan/biomarkers/platform_conversion_tables/GPL17692.txt", skip = 20)

gene_list <- c()
for (gene_str in GPL17692_conversion_table$gene_assignment) {
  if (grepl("//", gene_str)) {
    extracted_gene <- strsplit(gene_str, "//")[[1]][2]
    extracted_gene <- gsub(" ", "", extracted_gene)
    gene_list <- c(gene_list, extracted_gene)
     
  } else {
    gene_list <- c(gene_list, "")
  }
}

GPL17692_conversion_table$external_gene_name <- gene_list

GPL17692_conversion_table %<>% select(probeset_id, external_gene_name)
GPL17692_conversion_table$probeset_id <- as.character(GPL17692_conversion_table$probeset_id)
GPL17692_conversion_table
```

```{r}
expr_data <- as.data.frame(blood_biomarker_data_raw$originalData$GSE138118$expr)
expr_data %<>% rownames_to_column(var = "probeset_id") %>% left_join(GPL17692_conversion_table, by = c("probeset_id"="probeset_id")) %>% select(-probeset_id) %>%
  mutate(across(-external_gene_name, as.numeric)) %>% group_by(external_gene_name) %>% summarise(across(where(is.numeric), mean, na.rm = TRUE)) %>% ungroup() %>% filter(external_gene_name != "") %>% column_to_rownames(var="external_gene_name") 



# remove any samples that have been filtered out of the metadata
samples_in_metadata <- blood_biomarker_data_raw$originalData$GSE138118$pheno %>% pull(geo_accession)
expr_data %<>% select(!!!syms(samples_in_metadata))

# set keys
blood_biomarker_data_raw$originalData$GSE138118$keys <- rownames(expr_data)
blood_biomarker_data_raw$originalData$GSE138118$key_comment <- "Keys manually loaded"

# it says the samples were processed with RMA which defaults to log2, so these should be in log2 scale already
blood_biomarker_data_raw$originalData$GSE138118$expr <- as.matrix(expr_data)
```

#### GSE74629 
Tried loading the conversion table for the platform (	GPL10558	Illumina HumanHT-12 V4.0 expression beadchip) but got really weird data points (genes called "1-Mar" and "4-Mar" for example, which are clearly dates)


Using illumina wg_6_v3 to map probes to gene names
```{r}
ilm_gene_map <- illumina_gene_map
ilm_gene_map$probeset_id <- ilm_gene_map$illumina_humanwg_6_v3
ilm_gene_map %<>% select(probeset_id, external_gene_name) %>% drop_na()
```

```{r}
expr_data <- as.data.frame(blood_biomarker_data_raw$originalData$GSE74629$expr)
expr_data %<>% rownames_to_column(var = "probeset_id") %>% left_join(ilm_gene_map, by = c("probeset_id"="probeset_id")) %>% select(-probeset_id) %>%
  mutate(across(-external_gene_name, as.numeric)) %>% group_by(external_gene_name) %>% summarise(across(where(is.numeric), mean, na.rm = TRUE)) %>% ungroup() %>% filter(external_gene_name != "") %>% column_to_rownames(var="external_gene_name") 


# remove any samples that have been filtered out of the metadata
samples_in_metadata <- blood_biomarker_data_raw$originalData$GSE74629$pheno %>% pull(geo_accession)
expr_data %<>% select(!!!syms(samples_in_metadata))

# set keys
blood_biomarker_data_raw$originalData$GSE74629$keys <- rownames(expr_data)
blood_biomarker_data_raw$originalData$GSE74629$key_comment <- "Keys manually loaded"

# data was normalized using the neqc package which outputs in log2, doesn't need log transform.
blood_biomarker_data_raw$originalData$GSE74629$expr <- as.matrix(expr_data)

```


#### GSE125158

Set up probe <-> gene map
```{r}
agilent_gene_map$probeset_id <- agilent_gene_map$agilent_wholegenome_4x44k_v2
agilent_gene_map %<>% select(probeset_id, external_gene_name) %>% drop_na()
agilent_gene_map
```


```{r}
expr_data <- as.data.frame(blood_biomarker_data_raw$originalData$GSE125158$expr)
expr_data %<>% rownames_to_column(var = "probeset_id") %>% left_join(agilent_gene_map, by = c("probeset_id"="probeset_id")) %>% select(-probeset_id) %>%
  mutate(across(-external_gene_name, as.numeric)) %>% group_by(external_gene_name) %>% summarise(across(where(is.numeric), mean, na.rm = TRUE)) %>% ungroup() %>% filter(external_gene_name != "") %>% column_to_rownames(var="external_gene_name") 
expr_data

```

Remove Study (Values are normalized to GADPH by gene, not valid for this)
```{r}
blood_biomarker_data_raw$originalData$GSE125158 <- NULL
```


#### GSE20189
```{r}

study <- blood_biomarker_data_raw$originalData$GSE20189
samples_in_metadata <- study$pheno %>% pull(geo_accession)
expr_data <- as.data.frame(study$expr)

expr_data %<>% rownames_to_column(var="affy_hg_u133_plus_2") %>% left_join(ensembl_entrezgene_ensembl_gene_map %>% dplyr::select(affy_hg_u133_plus_2, external_gene_name), by=c("affy_hg_u133_plus_2"="affy_hg_u133_plus_2")) %>% filter(!is.na(external_gene_name)) %>% dplyr::select(-affy_hg_u133_plus_2) %>%
  mutate(across(-external_gene_name, as.numeric)) %>% group_by(external_gene_name) %>% summarise(across(where(is.numeric), mean, na.rm = TRUE)) %>% ungroup() %>% filter(external_gene_name != "")  %>% column_to_rownames(var="external_gene_name") 

# remove any samples that have been filtered out of the metadata
samples_in_metadata <- blood_biomarker_data_raw$originalData$GSE20189$pheno %>% pull(geo_accession)
expr_data %<>% select(!!!syms(samples_in_metadata))

blood_biomarker_data_raw$originalData$GSE20189$keys <- rownames(expr_data)
blood_biomarker_data_raw$originalData$GSE20189$key_comment <- "Keys manually loaded"

# save to the dataset
# doesn't need log transform "values are log scale RMA intensities"
blood_biomarker_data_raw$originalData$GSE20189$expr <- as.matrix(expr_data)
```


#### GSE16443
- Platform: GPL2986	ABI Human Genome Survey Microarray Version 2
- Data is log2 transformed

Load the conversion table 
```{r}
GPL2986_conversion_table <- read_tsv("/Users/hansenhan/biomarkers/platform_conversion_tables/GPL2986.txt", skip = 17)
GPL2986_conversion_table

GPL2986_conversion_table$external_gene_name <- GPL2986_conversion_table$`Gene Symbol`
GPL2986_conversion_table$probeset_id <-as.character(GPL2986_conversion_table$`ID`)
GPL2986_conversion_table %<>% select(external_gene_name, probeset_id) %>% mutate(external_gene_name = ifelse(external_gene_name == "null", NA, external_gene_name))

GPL2986_conversion_table
```

```{r}
expr_data <- as.data.frame(blood_biomarker_data_raw$originalData$GSE16443$expr)
expr_data %<>% rownames_to_column(var = "probeset_id") %>% left_join(GPL2986_conversion_table, by = c("probeset_id"="probeset_id")) %>% select(-probeset_id) %>%
  mutate(across(-external_gene_name, as.numeric)) %>% group_by(external_gene_name) %>% summarise(across(where(is.numeric), mean, na.rm = TRUE)) %>% ungroup() %>% filter(external_gene_name != "") %>% column_to_rownames(var="external_gene_name") 
expr_data


# remove any samples that have been filtered out of the metadata
samples_in_metadata <- blood_biomarker_data_raw$originalData$GSE16443$pheno %>% pull(geo_accession)
expr_data %<>% select(!!!syms(samples_in_metadata))

blood_biomarker_data_raw$originalData$GSE16443$keys <- rownames(expr_data)
blood_biomarker_data_raw$originalData$GSE16443$key_comment <- "Keys manually loaded"

# save to the dataset
blood_biomarker_data_raw$originalData$GSE16443$expr <- as.matrix(expr_data)

```


#### GSE164191
```{r}
expr_data <- as.data.frame(blood_biomarker_data_raw$originalData$GSE164191$expr)

expr_data %<>% rownames_to_column(var="affy_hg_u133_plus_2") %>% left_join(ensembl_entrezgene_ensembl_gene_map %>% dplyr::select(affy_hg_u133_plus_2, external_gene_name), by=c("affy_hg_u133_plus_2"="affy_hg_u133_plus_2")) %>% filter(!is.na(external_gene_name)) %>% dplyr::select(-affy_hg_u133_plus_2) %>%
  mutate(across(-external_gene_name, as.numeric)) %>% group_by(external_gene_name) %>% summarise(across(where(is.numeric), mean, na.rm = TRUE)) %>% ungroup() %>% filter(external_gene_name != "") %>% column_to_rownames(var="external_gene_name") 


# remove samples not in the metadata
samples_in_metadata <- blood_biomarker_data_raw$originalData$GSE164191$pheno %>% pull(geo_accession)
expr_data %<>% select(!!!syms(samples_in_metadata))


# save to the dataset
blood_biomarker_data_raw$originalData$GSE164191$expr <- as.matrix(expr_data)
blood_biomarker_data_raw$originalData$GSE164191$keys <- gene_names
blood_biomarker_data_raw$originalData$GSE164191$key_comment <- "Keys manually loaded"

```

#### GSE39345
- Data was quantile normalized
```{r}
expr_data <- as.data.frame(blood_biomarker_data_raw$originalData$GSE39345$expr)
expr_data %<>% rownames_to_column(var = "probeset_id") %>% left_join(ilm_gene_map, by = c("probeset_id"="probeset_id")) %>% select(-probeset_id) %>%
  mutate(across(-external_gene_name, as.numeric)) %>% group_by(external_gene_name) %>% summarise(across(where(is.numeric), mean, na.rm = TRUE)) %>% ungroup() %>% filter(external_gene_name != "") %>% column_to_rownames(var="external_gene_name") 


# remove any samples that have been filtered out of the metadata
samples_in_metadata <- blood_biomarker_data_raw$originalData$GSE39345$pheno %>% pull(geo_accession)
expr_data %<>% select(!!!syms(samples_in_metadata))

# set keys
blood_biomarker_data_raw$originalData$GSE39345$keys <- rownames(expr_data)
blood_biomarker_data_raw$originalData$GSE39345$key_comment <- "Keys manually loaded"

# log2 transform the data 
blood_biomarker_data_raw$originalData$GSE39345$expr <- log2(as.matrix(expr_data) + 1)

```

#### GSE49515
```{r}
expr_data <- as.data.frame(blood_biomarker_data_raw$originalData$GSE49515$expr)

expr_data %<>% rownames_to_column(var="affy_hg_u133_plus_2") %>% left_join(ensembl_entrezgene_ensembl_gene_map %>% dplyr::select(affy_hg_u133_plus_2, external_gene_name), by=c("affy_hg_u133_plus_2"="affy_hg_u133_plus_2")) %>% filter(!is.na(external_gene_name)) %>% dplyr::select(-affy_hg_u133_plus_2) %>%
  mutate(across(-external_gene_name, as.numeric)) %>% group_by(external_gene_name) %>% summarise(across(where(is.numeric), mean, na.rm = TRUE)) %>% ungroup() %>% filter(external_gene_name != "") %>% column_to_rownames(var="external_gene_name") 


# remove samples not in the metadata
samples_in_metadata <- blood_biomarker_data_raw$originalData$GSE49515$pheno %>% pull(geo_accession)
expr_data %<>% select(!!!syms(samples_in_metadata))


# save to the dataset
blood_biomarker_data_raw$originalData$GSE49515$expr <- as.matrix(expr_data)
blood_biomarker_data_raw$originalData$GSE49515$keys <- gene_names
blood_biomarker_data_raw$originalData$GSE49515$key_comment <- "Keys manually loaded"
```

#### GSE27562
```{r}
expr_data <- as.data.frame(blood_biomarker_data_raw$originalData$GSE27562$expr)

expr_data %<>% rownames_to_column(var="affy_hg_u133_plus_2") %>% left_join(ensembl_entrezgene_ensembl_gene_map %>% dplyr::select(affy_hg_u133_plus_2, external_gene_name), by=c("affy_hg_u133_plus_2"="affy_hg_u133_plus_2")) %>% filter(!is.na(external_gene_name)) %>% dplyr::select(-affy_hg_u133_plus_2) %>%
  mutate(across(-external_gene_name, as.numeric)) %>% group_by(external_gene_name) %>% summarise(across(where(is.numeric), mean, na.rm = TRUE)) %>% ungroup() %>% filter(external_gene_name != "") %>% column_to_rownames(var="external_gene_name") 


# remove samples not in the metadata
samples_in_metadata <- blood_biomarker_data_raw$originalData$GSE27562$pheno %>% pull(geo_accession)
expr_data %<>% select(!!!syms(samples_in_metadata))


# save to the dataset
blood_biomarker_data_raw$originalData$GSE27562$expr <- as.matrix(expr_data)
blood_biomarker_data_raw$originalData$GSE27562$keys <- gene_names
blood_biomarker_data_raw$originalData$GSE27562$key_comment <- "Keys manually loaded"
```

#### GSE10715
```{r}
expr_data <- as.data.frame(blood_biomarker_data_raw$originalData$GSE10715$expr)

expr_data %<>% rownames_to_column(var="affy_hg_u133_plus_2") %>% left_join(ensembl_entrezgene_ensembl_gene_map %>% dplyr::select(affy_hg_u133_plus_2, external_gene_name), by=c("affy_hg_u133_plus_2"="affy_hg_u133_plus_2")) %>% filter(!is.na(external_gene_name)) %>% dplyr::select(-affy_hg_u133_plus_2) %>%
  mutate(across(-external_gene_name, as.numeric)) %>% group_by(external_gene_name) %>% summarise(across(where(is.numeric), mean, na.rm = TRUE)) %>% ungroup() %>% filter(external_gene_name != "") %>% column_to_rownames(var="external_gene_name") 


# remove samples not in the metadata
samples_in_metadata <- blood_biomarker_data_raw$originalData$GSE10715$pheno %>% pull(geo_accession)
expr_data %<>% select(!!!syms(samples_in_metadata))


# save to the dataset
blood_biomarker_data_raw$originalData$GSE10715$expr <- as.matrix(expr_data)
blood_biomarker_data_raw$originalData$GSE10715$keys <- gene_names
blood_biomarker_data_raw$originalData$GSE10715$key_comment <- "Keys manually loaded"
```


## Check Dataset

```{r}
checkDataObject(blood_biomarker_data_raw, "Meta", "Pre-Analysis")
```

## Export Data
```{r}
save(blood_biomarker_data_raw, file="processed_cancer_blood_datasets.Rdata")
```