<jupyter_start><jupyter_text>---output: html_notebooktitle: "Exploratory Analysis"---<jupyter_code>source(rprojroot::is_git_root$find_file('lib/R/utils.R'))
config <- setup(packages=c('dplyr', 'GEOquery', 'ggplot2', 'jetset'))<jupyter_output><empty_output><jupyter_text>Choose dataset<jupyter_code>config$geo_datasets %>% paste(seq_along(.), ., sep=': ') %>% cat(sep=', ')
gse_id <- config$geo_datasets[6]

cat('Selected dataset:', gse_id)<jupyter_output><empty_output><jupyter_text>Load data<jupyter_code>eset <- load_data(gse_id)

print(eset)
pdata <- read.delim(sprintf('data/meta/%s.tsv', gse_id), check.names=FALSE, colClasses='character')
annot <- gse_id %>%
    parse_annotation(pdata) %>%
        as.data.frame() %>%
        cbind(pdata) %>%

annot %>%
    filter(her2 == 'HER2+') %>%
    select(treatment, outcome) %>%
    table()
platform <- c(GPL96='hgu133a', GPL570='hgu133plus2')
probes <- jmap(platform[attr(eset, 'annotation')], symbol=config$genes)

samples <- annot %>%
    filter(her2 == 'HER2+' & treatment != 'none' & !is.na(outcome)) %>%
    getElement('geo_accession')

dat <- gse_id %>%
    load_data() %>%
    exprs() %>%
    extract(probes, samples) %>%
    t() %>%
    as.data.frame() %>%
    setNames(config$genes) %>%
    cbind(select(annot, geo_accession, treatment, outcome))

print(head(dat))<jupyter_output><empty_output><jupyter_text>Expression levels<jupyter_code>plot_expr <- function(gene, dat) {
    p_value <- with(dat, {
        x <- dat[outcome == 'pCR', gene]
        y <- dat[outcome == 'RD', gene]
        t.test(x, y)$p.value
    })
    ggplot(dat) +
        geom_boxplot() +
        geom_jitter(position=position_jitter(width=0.1)) +
        aes_string(x='outcome', y=gene, fill='outcome') +
        ggtitle(sprintf('P-value for t-test: %.3f', p_value))
}

for (gene in config$genes) {
    print(plot_expr(gene, dat))
}<jupyter_output><empty_output><jupyter_text>Correlation between gene expression levels<jupyter_code>plot_corr <- function(pair, dat) {
    ggplot(dat) + geom_point() +
        aes_string(pair[1], pair[2], color='outcome') +
        ggtitle(sprintf('Correlation: %.2f', cor(dat[[pair[1]]], dat[[pair[2]]])))
}

gene_pairs <- combn(config$genes, 2, simplify=FALSE)

for (pair in gene_pairs) {
    print(plot_corr(pair, dat))
}<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1919758412424504, 'fraction_numerical': 0.011647972389991372, 'mean_word_length': 4.81203007518797, 'pattern_counts': {'":': 0, '<': 36, '<?xml version=': 0, '>': 39, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '653752', 'license': 'Unlicense', 'n_tokens_mistral': 936, 'n_tokens_neox': 868, 'n_words': 125, 'path': 'report/exploratory.ipynb', 'repo_name': 'leogama/phlda', 'size': '2318'}