\r\n",
"\r\n",
"
\r\n",
"
\r\n",
"
\r\n",
"\r\n",
"
\r\n",
"\r\n",
"
\r\n",
"\r\n",
"\r\n",
"
\r\n",
"
\r\n",
"
\r\n",
"\r\n",
"
\r\n",
"\r\n",
"
\r\n",
"\r\n",
"\r\n",
"\r\n",
"
\r\n",
"\r\n",
"
\r\n",
"\r\n",
"
\r\n",
"\r\n",
"
\r\n",
"\r\n",
"
\r\n",
"\r\n",
"
\r\n",
"\r\n",
"
\r\n",
"\r\n",
"\r\n",
"\r\n",
"\r\n",
"\r\n",
"\r\n",
"\r\n",
"\r\n",
"\r\n"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"#Run the following command to take the quiz\n",
"IRdisplay::display_html('
')"
]
},
{
"cell_type": "markdown",
"id": "74400910",
"metadata": {},
"source": [
"### Samples and groups selection\n",
"In this learning course, we focus on analyzing two groups: *disease* and *control* from the *entorhinal cortex* region of the ```GSE48350``` dataset . Therefore, we need to select the samples\n",
"that belong to *entorhinal cortex* region and form a group that has two classes. To select the samples that belong to the *entorhinal cortex*, we can use the following command:"
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "b805a146",
"metadata": {
"vscode": {
"languageId": "r"
}
},
"outputs": [],
"source": [
"# Select sample from Entorhinal Cortex region\n",
"idx <- grep(\"entorhinal\", samples$`brain region:ch1`)\n",
"samples = samples[idx,]\n",
"# Get expression and normalized expression data for samples collected from entorhinal cortex region\n",
"expression_data = expression_data[,idx]\n",
"norm_expression_data = norm_expression_data[,idx]"
]
},
{
"cell_type": "markdown",
"id": "33963c0d",
"metadata": {},
"source": [
"To check how many samples belong to the *entorhinal cortex* region, we can use the following command:"
]
},
{
"cell_type": "code",
"execution_count": 15,
"id": "5640f260",
"metadata": {
"vscode": {
"languageId": "r"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[1] \"#Genes: 54675 - #Samples: 54\"\n"
]
}
],
"source": [
"# Print out new number genes and samples\n",
"print(paste0(\"#Genes: \",dim(norm_expression_data)[1],\" - #Samples: \",dim(norm_expression_data)[2]))"
]
},
{
"cell_type": "markdown",
"id": "66b15c60",
"metadata": {},
"source": [
"In order to perform Differential Expression analysis in the later submodule, we need to group patients into two groups: “disease” and “control”. From the samples source name, patients diagnosed with Alzheimer’s Disease are annotated with ‘AD’.\n",
"Therefore, any samples that contain string ‘AD’ are labelled by ‘d’ (condition/disease) and the remaining samples are assigned to ‘c’ (control/normal).\n",
"The code to group the samples is presented below:"
]
},
{
"cell_type": "code",
"execution_count": 16,
"id": "9ccaca6f",
"metadata": {
"vscode": {
"languageId": "r"
}
},
"outputs": [],
"source": [
"# Select disease samples\n",
"disease_idx <- grep(\"AD\", samples$source_name_ch1)\n",
"# Create a vector to store label\n",
"groups <- rep(\"X\", nrow(samples))\n",
"# Annotate disease samples as \"d\"\n",
"groups[disease_idx] <- \"d\"\n",
"# Control samples are labeled as \"c\"\n",
"groups[which(groups!=\"d\")] <- \"c\"\n",
"groups <- factor(groups)"
]
},
{
"cell_type": "markdown",
"id": "a9bfe875",
"metadata": {},
"source": [
"To see the number of *control* and *disease*, we can use the following command:"
]
},
{
"cell_type": "code",
"execution_count": 17,
"id": "bff2353e",
"metadata": {
"vscode": {
"languageId": "r"
}
},
"outputs": [
{
"data": {
"text/plain": [
"groups\n",
" c d \n",
"39 15 "
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"table(groups)"
]
},
{
"cell_type": "markdown",
"id": "b0209857-88b4-422e-b64b-fbea32ebf1d1",
"metadata": {},
"source": [
"## Exporting the data\n",
"When we have successfully retrieved expression data from GEO, we can export the expression data to a `.csv` file format for inspection in other software such as Excel using the `write_csv` function from `readr` package.\n",
"In the code below, we will save the raw expression matrix, normalized expression matrix, and grouping information to `.csv` files."
]
},
{
"cell_type": "code",
"execution_count": 18,
"id": "0f39641e",
"metadata": {
"vscode": {
"languageId": "r"
}
},
"outputs": [],
"source": [
"# Convert raw and normalized expression matrix to datafames and save them to the csv file\n",
"expression_data <- as.data.frame(expression_data)\n",
"norm_expression_data <- as.data.frame(norm_expression_data)\n",
"# Convert group to a dataframe\n",
"groups <- as.data.frame(groups)\n",
"# Create a sub-directory data folder to save the expression matrix if it is not available\n",
"dir <- getwd()\n",
"subDir <- \"/data\"\n",
"path <- paste0(dir,subDir)\n",
"if (!file.exists(path)){\n",
" dir.create(file.path(path))\n",
"}\n",
"# Save expresion values and group to the csv files format in the local folder\n",
"write.csv(expression_data, file=\"./data/raw_GSE48350.csv\")\n",
"write.csv(norm_expression_data, file=\"./data/normalized_GSE48350.csv\")\n",
"write.csv(groups, file=\"./data/groups_GSE48350.csv\")"
]
},
{
"cell_type": "markdown",
"id": "23a731f6",
"metadata": {},
"source": [
"The `.csv` format is a very simple format that might not suitable to store big datasets.\n",
"We can export the expression data to `.rds` format, which is more memory efficient for loading and saving the data.\n",
"We can save all the relevant data in a `list` and write to the disk using the built in `saveRDS` function."
]
},
{
"cell_type": "code",
"execution_count": 19,
"id": "665869cc",
"metadata": {
"vscode": {
"languageId": "r"
}
},
"outputs": [],
"source": [
"# Putting raw data, normalized data, and groups into a list\n",
"dat <- list(expression_data = expression_data, norm_expression_data = norm_expression_data, groups = groups)\n",
"# Save the data to the local disk using rds format\n",
"saveRDS(dat, file=\"./data/GSE48350.rds\")\n"
]
},
{
"cell_type": "markdown",
"id": "9ade9695",
"metadata": {},
"source": [
"The user can copy the data from the local disk to the Bucket using the `gsutil cp` command on the Vertex AI Workbench terminal.\n",
"To do this, the user has to click on the `+` symbol beside the name of the current tab in the vertex AI workbench interface to open a new Launcher tab.\n",
"The next step is to click on `$_Terminal` open a Terminal and then type out the `gsutil cp` command below. The command to copy data to the Google Cloud Bucket is as follows:\n",
"\n",
"`gsutil cp PATH_TO_LOCAL_FOLDER gs://PATH_IN_GOOGLE_BUCKET`\n",
"\n",
"The steps to create a Google Bucket storage is shown in project GIT repository `https://github.com/tinnlab/NOSI-Google-Cloud-Training.git`. For example, if we want to copy the file `./data/GSE48350.rds` to our Google Bucket named `cpa-output`, we can use the following command."
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": false
},
"source": [
"```{admonition}
Tip: Saving data to the Google Cloud Bucket\n",
"gsutil cp ./data/GSE48350.rds gs://cpa-output\n",
"```"
]
},
{
"cell_type": "markdown",
"id": "1ef6c72b",
"metadata": {},
"source": [
"## Gene ID mapping\n",
"\n",
"\n",
"## Learning Objectives:\n",
"1. Understand different probe set ID.\n",
"2. Map probe IDs into gene identifiers and symbols."
]
},
{
"cell_type": "markdown",
"id": "1e8f26b4",
"metadata": {},
"source": [
"## Understanding different probe set ID\n",
"\n",
"Gene set or pathway analysis requires that gene sets and expression data use the same type of gene ID\n",
"(Entrez Gene ID, Gene symbol or probe set ID etc). However, this is frequently not true for the data we have.\n",
"For example, our gene sets mostly use Entrez Gene ID, but microarray datasets are frequently labeled by probe\n",
"set ID (or RefSeq transcript ID, etc.). Therefore, we need to map the probe set IDs to Entrez gene ID.\n",
"Here, we will use `GSE48350` dataset that we have used in the previous section for demonstration of gene ID mapping.\n",
"In order to know what kind of probe set ID, we need to navigate to GEO record page of `GSE48350`. Under the `Platform`\n",
"tab, we can find the probe ID information.\n",
""
]
},
{
"cell_type": "markdown",
"id": "3f2e5579",
"metadata": {},
"source": [
"From the record page, we know that the dataset was generated from 1 platform using the Affymetrix Human Genome U133 Plus 2.0 Array.\n",
"To convert or map the probe set IDs to Entrez gene ID, we need to find the corresponding annotation package from\n",
"Bioconductor. For analyzed data, we need to use the `hgu133plus2.db` and `AnnotationDbi` databases.\n",
"We can install the package using following R command:"
]
},
{
"cell_type": "code",
"execution_count": 22,
"id": "5fb6d5b0",
"metadata": {
"vscode": {
"languageId": "r"
}
},
"outputs": [],
"source": [
"suppressMessages({if (!require(\"BiocManager\", quietly = TRUE))\n",
" suppressWarnings(install.packages(\"BiocManager\"))\n",
" suppressWarnings(BiocManager::install(\"hgu133plus2.db\", update = F))\n",
" suppressWarnings(BiocManager::install(\"AnnotationDbi\", update = F))\n",
"})"
]
},
{
"cell_type": "markdown",
"id": "8bee6afb",
"metadata": {},
"source": [
"In the data processing section, we successfully downloaded the dataset `GSE48350` and saved it to the `data` sub-directory.\n",
"Now we can load the dataset and check for the Probe Set ID names by using the following command."
]
},
{
"cell_type": "code",
"execution_count": 23,
"id": "bae8d0ed",
"metadata": {
"vscode": {
"languageId": "r"
}
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"- '1007_s_at'
- '1053_at'
- '117_at'
- '121_at'
- '1255_g_at'
- '1294_at'
\n"
],
"text/latex": [
"\\begin{enumerate*}\n",
"\\item '1007\\_s\\_at'\n",
"\\item '1053\\_at'\n",
"\\item '117\\_at'\n",
"\\item '121\\_at'\n",
"\\item '1255\\_g\\_at'\n",
"\\item '1294\\_at'\n",
"\\end{enumerate*}\n"
],
"text/markdown": [
"1. '1007_s_at'\n",
"2. '1053_at'\n",
"3. '117_at'\n",
"4. '121_at'\n",
"5. '1255_g_at'\n",
"6. '1294_at'\n",
"\n",
"\n"
],
"text/plain": [
"[1] \"1007_s_at\" \"1053_at\" \"117_at\" \"121_at\" \"1255_g_at\" \"1294_at\" "
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"data = readRDS(\"./data/GSE48350.rds\")\n",
"expression_data = data$expression_data\n",
"probeIDs = rownames(expression_data)\n",
"head(probeIDs)"
]
},
{
"cell_type": "markdown",
"id": "7d6c37a9",
"metadata": {},
"source": [
"## Mapping probe IDs into gene identifiers and symbols.\n",
"For pathway analysis, databases like Gene Ontology and KEGG use gene symbol annotation. Therefore, Probe set IDs should be mapped to gene symbol to be used in the later submodules. We can use the pre-annotated databases which are available from Bioconductor to map the probe IDs to the gene symbol. `AnnotationDbi` is an R package that provides an interface for connecting and querying various annotation databases using SQLite data storage while `hgu133plus2.db` is database to perform Probe IDs to gene symbol mapping for human data. We can use the following command to map the Probe set IDs:"
]
},
{
"cell_type": "code",
"execution_count": 24,
"id": "e75ae09a",
"metadata": {
"vscode": {
"languageId": "r"
}
},
"outputs": [],
"source": [
"suppressMessages({\n",
" library(hgu133plus2.db)\n",
" library(AnnotationDbi)\n",
"})"
]
},
{
"cell_type": "code",
"execution_count": 25,
"id": "eb5a9160",
"metadata": {
"vscode": {
"languageId": "r"
}
},
"outputs": [],
"source": [
"suppressMessages({\n",
"annotLookup <- AnnotationDbi::select(hgu133plus2.db, keys = probeIDs, columns = c('PROBEID','GENENAME','SYMBOL'))\n",
"})"
]
},
{
"cell_type": "code",
"execution_count": 26,
"id": "ecea1cd0",
"metadata": {
"vscode": {
"languageId": "r"
}
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"A data.frame: 6 × 3\n",
"\n",
"\t | PROBEID | GENENAME | SYMBOL |
|---|
\n",
"\t | <chr> | <chr> | <chr> |
|---|
\n",
"\n",
"\n",
"\t| 1 | 1007_s_at | discoidin domain receptor tyrosine kinase 1 | DDR1 |
|---|
\n",
"\t| 2 | 1053_at | replication factor C subunit 2 | RFC2 |
|---|
\n",
"\t| 3 | 117_at | heat shock protein family A (Hsp70) member 6 | HSPA6 |
|---|
\n",
"\t| 4 | 121_at | paired box 8 | PAX8 |
|---|
\n",
"\t| 5 | 1255_g_at | guanylate cyclase activator 1A | GUCA1A |
|---|
\n",
"\t| 6 | 1294_at | ubiquitin like modifier activating enzyme 7 | UBA7 |
|---|
\n",
"\n",
"
\n"
],
"text/latex": [
"A data.frame: 6 × 3\n",
"\\begin{tabular}{r|lll}\n",
" & PROBEID & GENENAME & SYMBOL\\\\\n",
" & & & \\\\\n",
"\\hline\n",
"\t1 & 1007\\_s\\_at & discoidin domain receptor tyrosine kinase 1 & DDR1 \\\\\n",
"\t2 & 1053\\_at & replication factor C subunit 2 & RFC2 \\\\\n",
"\t3 & 117\\_at & heat shock protein family A (Hsp70) member 6 & HSPA6 \\\\\n",
"\t4 & 121\\_at & paired box 8 & PAX8 \\\\\n",
"\t5 & 1255\\_g\\_at & guanylate cyclase activator 1A & GUCA1A\\\\\n",
"\t6 & 1294\\_at & ubiquitin like modifier activating enzyme 7 & UBA7 \\\\\n",
"\\end{tabular}\n"
],
"text/markdown": [
"\n",
"A data.frame: 6 × 3\n",
"\n",
"| | PROBEID <chr> | GENENAME <chr> | SYMBOL <chr> |\n",
"|---|---|---|---|\n",
"| 1 | 1007_s_at | discoidin domain receptor tyrosine kinase 1 | DDR1 |\n",
"| 2 | 1053_at | replication factor C subunit 2 | RFC2 |\n",
"| 3 | 117_at | heat shock protein family A (Hsp70) member 6 | HSPA6 |\n",
"| 4 | 121_at | paired box 8 | PAX8 |\n",
"| 5 | 1255_g_at | guanylate cyclase activator 1A | GUCA1A |\n",
"| 6 | 1294_at | ubiquitin like modifier activating enzyme 7 | UBA7 |\n",
"\n"
],
"text/plain": [
" PROBEID GENENAME SYMBOL\n",
"1 1007_s_at discoidin domain receptor tyrosine kinase 1 DDR1 \n",
"2 1053_at replication factor C subunit 2 RFC2 \n",
"3 117_at heat shock protein family A (Hsp70) member 6 HSPA6 \n",
"4 121_at paired box 8 PAX8 \n",
"5 1255_g_at guanylate cyclase activator 1A GUCA1A\n",
"6 1294_at ubiquitin like modifier activating enzyme 7 UBA7 "
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"head(annotLookup)"
]
},
{
"cell_type": "code",
"execution_count": 27,
"id": "6aae161e",
"metadata": {
"tags": [
"remove-input"
],
"vscode": {
"languageId": "r"
}
},
"outputs": [
{
"data": {
"text/html": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"#Run the following command to take the quiz\n",
"IRdisplay::display_html('')"
]
},
{
"cell_type": "markdown",
"id": "c3f0cef5",
"metadata": {},
"source": [
"We can check the number of probe IDs that mapped to the unique gene symbols using the following command:"
]
},
{
"cell_type": "code",
"execution_count": 28,
"id": "0379ec35",
"metadata": {
"vscode": {
"languageId": "r"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[1] \"There are 54675 probe IDs that mapped to 22243 gene symbols\"\n"
]
}
],
"source": [
"# Check number of the probe IDs\n",
"print(paste0(\"There are \",length(unique(annotLookup$PROBEID)),\" probe IDs that mapped to \",length(unique(annotLookup$SYMBOL)),\" gene symbols\"))"
]
},
{
"cell_type": "markdown",
"id": "6247efff",
"metadata": {},
"source": [
"From the lookup table we can spot that a single gene symbol can be mapped to multiple probe IDs. In the next submodule, we will discuss how these mapped gene symbols can be used to identify which genes are differentially expressed. \n"
]
}
],
"metadata": {
"celltoolbar": "Tags",
"environment": {
"kernel": "ir",
"name": "r-cpu.4-1.m93",
"type": "gcloud",
"uri": "gcr.io/deeplearning-platform-release/r-cpu.4-1:m93"
},
"kernelspec": {
"display_name": "R",
"language": "R",
"name": "ir"
},
"language_info": {
"codemirror_mode": "r",
"file_extension": ".r",
"mimetype": "text/x-r-source",
"name": "R",
"pygments_lexer": "r",
"version": "4.2.1"
},
"toc-showcode": true
},
"nbformat": 4,
"nbformat_minor": 5
}