1 xCell attempt, Biopsy: 20170820

Caveat: xCell asks for gene Symbols.

ls()

##   [1] "adjusts"                       "ah"                           
##   [3] "annot_lm"                      "annot_lp"                     
##   [5] "chosen_colors"                 "fig_s1"                       
##   [7] "fig_s1_cds"                    "hs"                           
##   [9] "hs_adjusted"                   "hs_adjusted_fcqsl"            
##  [11] "hs_adjusted_fcqsl_metr"        "hs_annotations"               
##  [13] "hs_cds_expt"                   "hs_cds_macr"                  
##  [15] "hs_expt"                       "hs_fasta"                     
##  [17] "hs_final_annotations"          "hs_gff"                       
##  [19] "hs_gff_annot"                  "hs_go_biomart"                
##  [21] "hs_lengths"                    "hs_macr"                      
##  [23] "hs_macr_batchnorm_metrics"     "hs_macr_cpmlow"               
##  [25] "hs_macr_cpmlow_pca"            "hs_macr_cpmlowcbcb"           
##  [27] "hs_macr_cpmlowcbcb_pca"        "hs_macr_cpmrlelow"            
##  [29] "hs_macr_cpmrlelow_pca"         "hs_macr_filt"                 
##  [31] "hs_macr_l2cpmlow"              "hs_macr_l2cpmlow_pca"         
##  [33] "hs_macr_l2cpmlowcom"           "hs_macr_l2cpmlowcom_pca"      
##  [35] "hs_macr_l2cpmpofa"             "hs_macr_l2cpmpofa_pca"        
##  [37] "hs_macr_l2cpmquantlowcmod"     "hs_macr_l2cpmquantlowcmod_pca"
##  [39] "hs_macr_l2cpmquantlowcom"      "hs_macr_l2cpmquantlowfsva"    
##  [41] "hs_macr_l2cpmquantlowfsva_pca" "hs_macr_l2cpmquantlowruv"     
##  [43] "hs_macr_l2cpmquantlowsva"      "hs_macr_l2cpmrlelow"          
##  [45] "hs_macr_l2cpmrlelow_pca"       "hs_macr_l2cpmrlelowcom"       
##  [47] "hs_macr_l2cpmrlelowcom_pca"    "hs_macr_metrics"              
##  [49] "hs_macr_nonil"                 "hs_macr_pca"                  
##  [51] "hs_macr_pca_info"              "hs_macr_se"                   
##  [53] "hs_macr_vp_cb"                 "hs_macr_vp_cbs3"              
##  [55] "hs_sv_adjusted"                "hs_sv_adjusted_fcqsl"         
##  [57] "hs_sv_adjusted_fcqsl_metr"     "labels"                       
##  [59] "lb_annotations"                "lb_fasta"                     
##  [61] "lb_gff"                        "lb_goids"                     
##  [63] "lb_lengths"                    "lb_tooltips"                  
##  [65] "lp_annotations"                "lp_fasta"                     
##  [67] "lp_gff"                        "lp_goids"                     
##  [69] "lp_lengths"                    "lp_macr"                      
##  [71] "lp_macr_fcqcl"                 "lp_macr_fcqcl_metr"           
##  [73] "lp_macr_fcqfl"                 "lp_macr_fcqfl_metr"           
##  [75] "lp_macr_fcqrr"                 "lp_macr_fcqrr_metr"           
##  [77] "lp_macr_fcrfl"                 "lp_macr_fcrfl_metr"           
##  [79] "lp_macr_limma_fcqlr"           "lp_macr_lpqcf"                
##  [81] "lp_macr_lpqcf_metr"            "lp_macr_rprrf"                
##  [83] "lp_macr_rrrrr_metr"            "lp_macr_sva"                  
##  [85] "lp_macr_sva_norm"              "lp_macr_sva_norm_metr"        
##  [87] "lp_macr_vp_cb"                 "lp_macr_vp_hsfcqsl"           
##  [89] "lp_macr_vp_lpqcf"              "lp_testing"                   
##  [91] "lp_tooltips"                   "lpan"                         
##  [93] "lpan_ah"                       "lpan_org"                     
##  [95] "new_adjusts"                   "new_colors"                   
##  [97] "note"                          "old_options"                  
##  [99] "orgdbs"                        "parasite_expt"                
## [101] "pretty_pca"                    "previous_file"                
## [103] "rmd_file"                      "sv_adjusted"                  
## [105] "test_adjusted"                 "this_save"                    
## [107] "tmp"                           "tt"                           
## [109] "txtdbs"                        "ver"                          
## [111] "written"

xcell_input <- normalize_expt(hs_expt, norm="quant", convert="cpm", transform="log2")

## This function will replace the expt$expressionset slot with:

## log2(cpm(quant(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

## Filter is false, this should likely be set to something, good
##  choices include cbcb, kofa, pofa (anything but FALSE).  If you want this to
##  stay FALSE, keep in mind that if other normalizations are performed, then the
##  resulting libsizes are likely to be strange (potentially negative!)

## Not correcting the count-data for batch effects.  If batch is
##  included in EdgerR/limma's model, then this is probably wise; but in extreme
##  batch effects this is a good parameter to play with.

## Step 1: not doing count filtering.

## Step 2: normalizing the data with quant.

## Using normalize.quantiles.robust due to a thread error in preprocessCore.

## Step 3: converting the data with cpm.

## Step 4: transforming the data with log2.

## transform_counts: Found 1447108 values equal to 0, adding 1 to the matrix.

## Step 5: not doing batch correction.

xcell_input <- exprs(xcell_input)

testing <- load_biomart_annotations()

## The biomart annotations file already exists, loading from it.

xref <- testing[["annotation"]][, c("ensembl_gene_id", "hgnc_symbol")]
## Lets get the full set of downloadable attributes
xcell_test <- merge(xcell_input, xref, by.x="row.names", by.y="ensembl_gene_id")
rownames(xcell_test) <- make.names(xcell_test[["hgnc_symbol"]], unique=TRUE)
xcell_test[["Row.names"]] <- NULL
xcell_test[["hgnc_symbol"]] <- NULL

library(xCell)
data("xCell.data", package="xCell")
summary(xCell.data)

##             Length Class             Mode     
## spill           3  -none-            list     
## spill.array     3  -none-            list     
## signatures    489  GeneSetCollection list     
## genes       10808  -none-            character

library(GSEABase)

## Loading required package: annotate

## Loading required package: AnnotationDbi

## Loading required package: stats4

## Loading required package: IRanges

## Loading required package: S4Vectors

## 
## Attaching package: 'S4Vectors'

## The following object is masked from 'package:base':
## 
##     expand.grid

## Loading required package: XML

## Loading required package: graph

## 
## Attaching package: 'graph'

## The following object is masked from 'package:XML':
## 
##     addNode

details(xCell.data$signatures[[1]])

## setName: aDC%HPCA%1.txt 
## geneIds: C1QA, C1QB, ..., CCL22 (total: 8)
## geneIdType: Null
## collectionType: Null 
## setIdentifier: PEDS-092FVH8-LT:623:Tue Jun  6 14:36:33 2017:2
## description: 
## organism: 
## pubMedIds: 
## urls: 
## contributor: 
## setVersion: 0.0.1
## creationDate:

result <- xCellAnalysis(xcell_test)

## [1] "Num. of genes: 10658"
## Estimating ssGSEA scores for 489 gene sets.
## 
  |                                                                       
  |                                                                 |   0%Using parallel with 4 cores
## 
  |                                                                       
  |=                                                                |   2%
  |                                                                       
  |===                                                              |   4%
  |                                                                       
  |====                                                             |   6%
  |                                                                       
  |=====                                                            |   8%
  |                                                                       
  |======                                                           |  10%
  |                                                                       
  |========                                                         |  12%
  |                                                                       
  |=========                                                        |  14%
  |                                                                       
  |==========                                                       |  16%
  |                                                                       
  |============                                                     |  18%
  |                                                                       
  |=============                                                    |  20%
  |                                                                       
  |==============                                                   |  22%
  |                                                                       
  |================                                                 |  24%
  |                                                                       
  |=================                                                |  26%
  |                                                                       
  |==================                                               |  28%
  |                                                                       
  |====================                                             |  30%
  |                                                                       
  |=====================                                            |  32%
  |                                                                       
  |======================                                           |  34%
  |                                                                       
  |=======================                                          |  36%
  |                                                                       
  |=========================                                        |  38%
  |                                                                       
  |==========================                                       |  40%
  |                                                                       
  |===========================                                      |  42%
  |                                                                       
  |=============================                                    |  44%
  |                                                                       
  |==============================                                   |  46%
  |                                                                       
  |===============================                                  |  48%
  |                                                                       
  |================================                                 |  50%
  |                                                                       
  |==================================                               |  52%
  |                                                                       
  |===================================                              |  54%
  |                                                                       
  |====================================                             |  56%
  |                                                                       
  |======================================                           |  58%
  |                                                                       
  |=======================================                          |  60%
  |                                                                       
  |========================================                         |  62%
  |                                                                       
  |==========================================                       |  64%
  |                                                                       
  |===========================================                      |  66%
  |                                                                       
  |============================================                     |  68%
  |                                                                       
  |==============================================                   |  70%
  |                                                                       
  |===============================================                  |  72%
  |                                                                       
  |================================================                 |  74%
  |                                                                       
  |=================================================                |  76%
  |                                                                       
  |===================================================              |  78%
  |                                                                       
  |====================================================             |  80%
  |                                                                       
  |=====================================================            |  82%
  |                                                                       
  |=======================================================          |  84%
  |                                                                       
  |========================================================         |  86%
  |                                                                       
  |=========================================================        |  88%
  |                                                                       
  |==========================================================       |  90%
  |                                                                       
  |============================================================     |  92%
  |                                                                       
  |=============================================================    |  94%
  |                                                                       
  |==============================================================   |  96%
  |                                                                       
  |================================================================ |  98%
  |                                                                       
  |=================================================================| 100%

pp(file="images/xcell_signatures.png")

## Going to write the image to: images/xcell_signatures.png when dev.off() is called.

heatmap.3(result, trace="none")
dev.off()

## png 
##   2

pp(file="images/xcell_boxplot.png", image=hpgltools::plot_boxplot(t(result)))

## Writing the image to: images/xcell_boxplot.png and calling dev.off().

test_table <- hs_macr_sva$limma$all_tables$sh_vs_chr

## Error in eval(expr, envir, enclos): object 'hs_macr_sva' not found

library(GSEABase)
library(GSVA)
library(GSVAdata)

## Loading required package: hgu95a.db

## Loading required package: org.Hs.eg.db

## 

## 

data(c2BroadSets)

##leukemia_es <- gsva(leukemia_filtered_eset, c2BroadSets,
##                    min.sz=10, max.sz=500, verbose=TRUE)
eset <- hs_expt$expressionset
ensembl_to_entrez <- sm(select(x=hs, keys=rownames(exprs(eset)), keytype="ENSEMBL", columns=c("ENTREZID")))

## Error in get(str): object 'GO.db' not found

## The broad IDs are all Entrez, but mine are Ensembl.
old_ids <- exprs(eset)
new_ids <- merge(old_ids, ensembl_to_entrez, by.x="row.names", by.y="ENSEMBL")

## Error in merge(old_ids, ensembl_to_entrez, by.x = "row.names", by.y = "ENSEMBL"): object 'ensembl_to_entrez' not found

remaining_ids <- new_ids[["Row.names"]]

## Error in eval(expr, envir, enclos): object 'new_ids' not found

new_eset <- eset[remaining_ids, ]

## Error in eset[remaining_ids, ]: object 'remaining_ids' not found

rownames(new_eset) <- make.names(new_ids[["ENTREZID"]], unique=TRUE)

## Error in make.names(new_ids[["ENTREZID"]], unique = TRUE): object 'new_ids' not found

rownames(new_eset) <- gsub(pattern="^X", replacement="", x=rownames(new_eset))

## Error in rownames(new_eset): object 'new_eset' not found

eset <- new_eset

## Error in eval(expr, envir, enclos): object 'new_eset' not found

fData(eset)$ENTREZID <- rownames(fData(eset))
annotation(eset) <- "org.Hs.eg.db"


geneset_test <- GeneSet(eset, setName="annoying")
filtered_eset <- genefilter::nsFilter(eset, require.entrez=TRUE, remove.dupEntrez=FALSE,
                                      var.func=IQR, var.filter=TRUE, var.cutoff=0.5, filterByQuantile=TRUE,
                                      feature.exclude="^AFFX")
gsva_test <- gsva(eset, c2BroadSets, verbose=TRUE)

## Warning in .local(expr, gset.idx.list, ...): 14510 genes with constant
## expression values throuhgout the samples.

## Warning in .local(expr, gset.idx.list, ...): Since argument method!
## ="ssgsea", genes with constant expression values are discarded.

## Mapping identifiers between gene sets and feature names

## Error in .gsva(exprs(expr), mapped.gset.idx.list, method, kcdf, rnaseq, : The gene set list is empty!  Filter may be too stringent.

tsne_test <- plot_tsne(gsva_test)

## Error in plot_pca(..., pc_method = "tsne"): object 'gsva_test' not found

## Testing obnoxious gsva examples which are unbelievably unhelpful
leuk_test <- genefilter::nsFilter(leukemia_eset, require.entrez=TRUE)

## Error in genefilter::nsFilter(leukemia_eset, require.entrez = TRUE): object 'leukemia_eset' not found

data(sample.ExpressionSet) # from Biobase
egs <- GeneSet(sample.ExpressionSet[201:250,], setName="Sample")
gsva_test <- gsva(egs, c2BroadSets, verbose=TRUE)

## Error in (function (classes, fdef, mtable) : unable to find an inherited method for function 'gsva' for signature '"GeneSet", "GeneSetCollection"'

all_gsva_categories <- names(c2BroadSets)
all_prefixes <- gsub(pattern="^(.*?)_.*", replacement="\\1", x=all_gsva_categories)
unique_prefixes <- unique(sort(all_prefixes))

reactome_subset <- grepl(x=rownames(gsva_test), pattern="^REACTOME")

## Error in rownames(gsva_test): object 'gsva_test' not found

reactome_gsva <- gsva_test[reactome_subset, ]

## Error in eval(expr, envir, enclos): object 'gsva_test' not found

pp("/tmp/testing.svg")

## Going to write the image to: /tmp/testing.svg when dev.off() is called.

##tt <- heatmap.3(exprs(reactome_gsva), trace="none", cexRow=0.1, cexCol=0.5)
tt <- heatmap.3(exprs(reactome_gsva), cexRow=0.1, cexCol=0.5)

## Error in exprs(reactome_gsva): object 'reactome_gsva' not found

dev.off()

## png 
##   2

original <- exprs(reactome_gsva)

## Error in exprs(reactome_gsva): object 'reactome_gsva' not found

new_order <- tt$rowInd
new <- original[new_order, ]

## Error in eval(expr, envir, enclos): object 'original' not found

gsva_result <- simple_gsva(hs_expt)

## Error in if (length(annotation(eset)) == 0 | grep(pattern = "Fill me in", : argument is of length zero

new_expt <- gsva_result[["expt"]]

## Error in eval(expr, envir, enclos): object 'gsva_result' not found

gsva_dis <- plot_disheat(new_expt)

## Error in plot_heatmap(expt_data, expt_colors = expt_colors, expt_design = expt_design, : object 'new_expt' not found

pp(file="images/gsva_all.png")

## Going to write the image to: images/gsva_all.png when dev.off() is called.

tt <- heatmap.3(exprs(new_expt), cexRow=0.1, trace="none", cexCol=0.5)

## Error in exprs(new_expt): object 'new_expt' not found

dev.off()

## png 
##   2

all_gsva_categories <- names(c2BroadSets)
all_prefixes <- gsub(pattern="^(.*?)_.*", replacement="\\1", x=all_gsva_categories)
unique_prefixes <- unique(sort(all_prefixes))

reactome_subset <- grepl(x=rownames(new_expt$expressionset), pattern="^REACTOME")

## Error in rownames(new_expt$expressionset): object 'new_expt' not found

reactome_gsva <- new_expt$expressionset[reactome_subset, ]

## Error in eval(expr, envir, enclos): object 'new_expt' not found

pp("images/reactome_gsva.png")

## Going to write the image to: images/reactome_gsva.png when dev.off() is called.

tt <- heatmap.3(exprs(reactome_gsva), cexRow=0.1, cexCol=0.5, trace="none")

## Error in exprs(reactome_gsva): object 'reactome_gsva' not found

dev.off()

## png 
##   2

if (!isTRUE(get0("skip_load"))) {
  pander::pander(sessionInfo())
  message(paste0("This is hpgltools commit: ", get_git_commit()))
  this_save <- paste0(gsub(pattern="\\.Rmd", replace="", x=rmd_file), "-v", ver, ".rda.xz")
  message(paste0("Saving to ", this_save))
  tmp <- sm(saveme(filename=this_save))
}

## If you wish to reproduce this exact build of hpgltools, invoke the following:

## > git clone http://github.com/abelew/hpgltools.git

## > git reset 3bc454f6c7093e2716d7a9bc1a01f9c15eb49101

## This is hpgltools commit: Wed Jan 2 16:23:19 2019 -0500: 3bc454f6c7093e2716d7a9bc1a01f9c15eb49101

## Saving to 02_expression_macrophage-v20170820.rda.xz