index.html

cpsy <- set_expt_colors(cpsy)
cpsy_metrics <- sm(graph_metrics(cpsy))
cpsy_norm <- sm(normalize_expt(cpsy, transform="log2", convert="cpm", norm="quant", filter=TRUE))
cpsy_norm_metrics <- sm(graph_metrics(cpsy_norm))
tt = set_expt_batch(cpsy, fact="replicate")
cpsy_normv2 <- normalize_expt(tt, transform="log2", convert="cpm", norm="quant", filter=TRUE, batch="combat_scale")
## This function will replace the expt$expressionset slot with:
## log2(combat_scale(cpm(quant(filter(data)))))
## It backs up the current data into a slot named:
##  expt$backup_expressionset. It will also save copies of each step along the way
##  in expt$normalized with the corresponding libsizes. Keep the libsizes in mind
##  when invoking limma.  The appropriate libsize is the non-log(cpm(normalized)).
##  This is most likely kept at:
##  'new_expt$normalized$intermediate_counts$normalization$libsizes'
##  A copy of this may also be found at:
##  new_expt$best_libsize
## Step 1: performing count filter with option: cbcb
## Removing 43 low-count genes (1771 remaining).
## Step 2: normalizing the data with quant.
## Step 3: converting the data with cpm.
## Step 4: transforming the data with log2.
## Performing batch correction at step 5.
## Step 5: doing batch correction with combat_scale.
## Note to self:  If you get an error like 'x contains missing values'; I think this means that the data has too many 0's and needs to have a better low-count filter applied.
## batch_counts: Before batch correction, 431 entries 0<x<1.
## batch_counts: Before batch correction, 712 entries are >= 0.
## batch_counts: Using sva::combat with a prior for batch correction and with scaling.
## Warning in do_batch(count_table, ...): The batch_counts call failed. Returning non-batch
## reduced data.
cpsy_normv2_metrics <- sm(graph_metrics(cpsy_normv2))
cpsy_normv3 <- normalize_expt(tt, transform="log2", convert="cpm", norm="quant", filter=TRUE, batch="sva")
## This function will replace the expt$expressionset slot with:
## log2(sva(cpm(quant(filter(data)))))
## It backs up the current data into a slot named:
##  expt$backup_expressionset. It will also save copies of each step along the way
##  in expt$normalized with the corresponding libsizes. Keep the libsizes in mind
##  when invoking limma.  The appropriate libsize is the non-log(cpm(normalized)).
##  This is most likely kept at:
##  'new_expt$normalized$intermediate_counts$normalization$libsizes'
##  A copy of this may also be found at:
##  new_expt$best_libsize
## Step 1: performing count filter with option: cbcb
## Removing 43 low-count genes (1771 remaining).
## Step 2: normalizing the data with quant.
## Step 3: converting the data with cpm.
## Step 4: transforming the data with log2.
## Performing batch correction at step 5.
## Step 5: doing batch correction with sva.
## Note to self:  If you get an error like 'x contains missing values'; I think this means that the data has too many 0's and needs to have a better low-count filter applied.
## batch_counts: Before batch correction, 431 entries 0<x<1.
## batch_counts: Before batch correction, 712 entries are >= 0.
## batch_counts: Using sva::fsva for batch correction.
## Number of significant surrogate variables is:  2 
## Iteration (out of 5 ):1  2  3  4  5
## The number of elements which are < 0 after batch correction is: 630
## The variable low_to_zero sets whether to change <0 values to 0 and is: FALSE
cpsy_normv3_metrics <- sm(graph_metrics(cpsy_normv3))
cpsy_normv3_metrics$pcaplot
cpsy_metrics$libsize

cpsy_metrics$density

cpsy_norm_metrics$corheat

cpsy_norm_metrics$pcaplot

cpsy_norm_metrics$smc

cpsy_normv2_metrics$corheat

cpsy_normv2_metrics$pcaplot

cpsy_normv2_metrics$smc

Holy ass crackers, the data is beautiful.

## It turns out that doing these analyses in parallel does weird and bad things to knitr.
## When running interactively, parallel may be true, otherwise no.
parallel <- FALSE

## First try a normal condition + batch model
ver <- "02"
keepers <- list(
    "cpsyvswt_thy" = c("cpsy", "wt"),
    "cpsyvswt_plasma" = c("cpsy_plasma","wt_plasma"),  ## I think this one will not work well
    "cpsy_plasmavsthy" = c("cpsy_plasma", "cpsy"),
    "wt_plasmavsthy" = c("wt_plasma", "wt"))


## Try without any batch awareness, this is terribad
cpsy_pairwise_nobatch <- sm(all_pairwise(cpsy, model_batch=FALSE, parallel=parallel))
cpsy_tables <- sm(combine_de_tables(cpsy_pairwise, keepers=keepers, extra_annot=Biobase::exprs(cpsy_norm$expressionset),
                                 excel=paste0("excel/cpsy_nobatch_tables-", ver, ".xlsx")))
cpsy_sig <- sm(extract_significant_genes(cpsy_tables, according_to="all",
                                         excel=paste0("excel/cpsy_nobatch_sig-", ver, ".xlsx")))
## Then try with batch in the model
cpsy_pairwise <- sm(all_pairwise(cpsy, model_batch=TRUE, parallel=parallel))
cpsy_tables <- sm(combine_de_tables(cpsy_pairwise, keepers=keepers, extra_annot=Biobase::exprs(cpsy_norm$expressionset),
                                    excel=paste0("excel/cpsy_batchmodel_tables-", ver, ".xlsx"))
cpsy_sig <- sm(extract_significant_genes(cpsy_tables, according_to="all",
                                         excel=paste0("excel/cpsy_batchmodel_sig-", ver, ".xlsx"))))
## Be a little more fussy with the data via sva
cpsy_pairwise <- sm(all_pairwise(cpsy, model_batch="sva", parallel=parallel))
cpsy_tables <- sm(combine_de_tables(cpsy_pairwise, keepers=keepers, extra_annot=Biobase::exprs(cpsy_norm$expressionset),
                                    excel=paste0("excel/cpsy_sva_tables-", ver, ".xlsx")))
cpsy_sig <- sm(extract_significant_genes(cpsy_tables, according_to="all",
                                         excel=paste0("excel/cpsy_sva_sig-", ver, ".xlsx")))
cpsy <- set_expt_batch(expt=cpsy, fact="replicate")
## Now be an asshat and fuck with the data
cpsy_combat <- sm(normalize_expt(cpsy, convert="cpm", norm="quant", filter=TRUE, batch="combat_scale"))
cpsy_pairwise <- sm(all_pairwise(cpsy_combat, model_batch=FALSE, force=TRUE, parallel=parallel)
                    cpsy_tables <- combine_de_tables(cpsy_pairwise, keepers=keepers, extra_annot=Biobase::exprs(cpsy_norm$expressionset),
                                                     excel=paste0("excel/cpsy_combat_tables-", ver, ".xlsx")))
cpsy_sig <- sm(extract_significant_genes(cpsy_tables, according_to="all",
                                         excel=paste0("excel/cpsy_combat_sig-", ver, ".xlsx")))
## Error: <text>:24:1: unexpected symbol
## 23:                                     excel=paste0("excel/cpsy_batchmodel_tables-", ver, ".xlsx"))
## 24: cpsy_sig
##     ^