## ----setup, echo=FALSE, include=FALSE-----------------------------------------
pkg <- read.dcf("../DESCRIPTION", fields = "Package")[1]
library(pkg, character.only = TRUE)

## ----install, eval = FALSE----------------------------------------------------
# if (!requireNamespace("BiocManager", quietly = TRUE)) {
#     install.packages("BiocManager")
# }
# 
# BiocManager::install("AnnotatedBCGEData")

## ----install_dev, eval = FALSE------------------------------------------------
# BiocManager::install("AnnotatedBCGEData", version="devel")

## ----load_this_package, eval = FALSE------------------------------------------
# library(AnnotatedBCGEData)

## ----load_other_packages, include = FALSE-------------------------------------
library(ggplot2)
library(dplyr)
library(tibble)
library(tidyr)
library(zen4R)
library(AnnotatedBCGEData)
library(SummarizedExperiment)

## ----search_PRS---------------------------------------------------------------
searchOntologyTerms("Progesterone Receptor Status")

## ----search_RS----------------------------------------------------------------
searchOntologyTerms("receptor status")

## ----search_PgR---------------------------------------------------------------
searchOntologyTerms("PgR Status")

## ----search_ihc---------------------------------------------------------------
searchOntologyTerms("immunohistochemical", term_type = "Definition")

## ----search_by_field----------------------------------------------------------
searchForDatasetsByField(c("C16149", "C16150"))

## ----search_by_value----------------------------------------------------------
searchForDatasetsByValue("C15496")

## ----get_metadata-------------------------------------------------------------
getDatasetMetadata(c("GSE41197", "GSE59772", "GSE10797"))

## ----get_dataset, eval = requireNamespace('AnnotatedBCGEData')----------------
GSE41197_SE = getDataset("GSE41197")

expression_data = assay(GSE41197_SE)
sample_metadata = colData(GSE41197_SE)
feature_data = rowData(GSE41197_SE)

## ----get_dataset_version------------------------------------------------------
GSE41197_v2 = getDataset("GSE41197", v=2)

v2_expression_data = assay(GSE41197_v2)
v2_metadata = colData(GSE41197_v2)
v2_feature_data = rowData(GSE41197_v2)

## ----boxplot------------------------------------------------------------------
# Set the random seed so the graph is reproducible.
set.seed(0)

as_tibble(expression_data, rownames='Ensembl_Gene_ID') %>%
  head(n = 10) %>%
  pivot_longer(-Ensembl_Gene_ID, names_to = "Sample_ID", values_to = "Expression_Level") %>%
  ggplot(aes(x=Ensembl_Gene_ID, y=Expression_Level)) +
  geom_boxplot() +
  geom_jitter(color = "red", alpha = 0.2) +
  labs(x='Ensembl Gene ID', y='Expression level') +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

## ----heatmap_1, message=FALSE, warning=FALSE----------------------------------
# Retrieve the data from Zenodo
GSE10797_SE <- getDataset("GSE10797")
GSE59772_SE <- getDataset("GSE59772")

# Retrieve stroma sample IDs as vectors
GSE10797_stroma_sample_ids = colData(GSE10797_SE) %>%
  as_tibble(rownames = 'Sample_ID') %>%
  filter(tissue_source == 'stromal cells from breast cancer patient') %>%
  pull(Sample_ID)

GSE59772_stroma_sample_ids = colData(GSE59772_SE) %>%
  as_tibble(rownames = 'Sample_ID') %>%
  filter(tissue == 'Stroma') %>%
  pull(Sample_ID)

GSE10797_SE = GSE10797_SE[,colnames(GSE10797_SE) %in% GSE10797_stroma_sample_ids]
GSE59772_SE = GSE59772_SE[,colnames(GSE59772_SE) %in% GSE59772_stroma_sample_ids]

## ----heatmap_2----------------------------------------------------------------
common_genes <- intersect(rownames(GSE10797_SE), rownames(GSE59772_SE))

GSE10797_SE <- GSE10797_SE[common_genes,]
GSE59772_SE <- GSE59772_SE[common_genes,]
# This also has the benefit that the genes will be in the same order in both.

## ----heatmap_3----------------------------------------------------------------
combined_data = cbind(assay(GSE10797_SE), assay(GSE59772_SE)) %>%
  as_tibble(rownames = 'Gene_ID') %>%
  pivot_longer(-Gene_ID, names_to = 'Sample_ID', values_to = 'Expression_Level')

top_variance_genes = group_by(combined_data, Gene_ID) %>%
  summarize(Variance = var(Expression_Level)) %>%
  arrange(desc(Variance)) %>%
  head(n = 10) %>%
  pull(Gene_ID)

combined_data = filter(combined_data, Gene_ID %in% top_variance_genes)

ggplot(combined_data, aes(x = Sample_ID, y = Gene_ID, fill = Expression_Level)) +
  geom_tile() +
  labs(x = 'Sample', y = 'Gene', fill = 'Expression level', 
       title = 'Genes in stromal samples with highest variance') + 
  theme_bw() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1), 
        plot.title = element_text(hjust = 0.5, face = "bold"))

## ----boxplot2_1, warning=FALSE------------------------------------------------
gene_ids = rowData(GSE10797_SE) %>%
  as_tibble(rownames = 'Ensembl_Gene_ID') %>%
  filter(Chromosome %in% 1:2 & Gene_Biotype == "lncRNA") %>%
  pull(Ensembl_Gene_ID)

GSE10797_SE <- GSE10797_SE[gene_ids,]
GSE59772_SE <- GSE59772_SE[gene_ids,]

## ----boxplot2_2---------------------------------------------------------------
gse10797_data = assay(GSE10797_SE) %>%
  as_tibble(rownames = 'Gene_ID') %>%
  pivot_longer(-Gene_ID, names_to = 'Sample_ID', values_to = 'Expression_Level') %>%
  mutate(Dataset_ID = "GSE10797")

gse59772_data = assay(GSE59772_SE) %>%
  as_tibble(rownames = 'Gene_ID') %>%
  pivot_longer(-Gene_ID, names_to = 'Sample_ID', values_to = 'Expression_Level') %>%
  mutate(Dataset_ID = "GSE59772")

combined_data = bind_rows(gse10797_data, gse59772_data)

## ----boxplot2_3---------------------------------------------------------------
ggplot(combined_data, aes(x = Gene_ID, y = Expression_Level, fill = Dataset_ID)) +
  geom_boxplot() +
  theme_bw() +
  labs(x = 'Gene', y='Expression Level', 
       title = 'lncRNA expression levels on chromosomes 1 & 2') +
  coord_flip()

## ----Session Info-------------------------------------------------------------
utils::sessionInfo()