I.scapularis 2017: Tick Sample Estimation.

1 Sample Estimation, Mus musculus: 20180119

This document is concerned with analyzing RNAseq data of iscapularis ixodes (tick).

2 Initial tick estimation

isc_plots <- graph_metrics(ixo_expt)

## Graphing number of non-zero genes with respect to CPM by library.

## Graphing library sizes.

## Graphing a boxplot.

## This data will benefit from being displayed on the log scale.

## If this is not desired, set scale='raw'

## Some entries are 0.  We are on log scale, adding 1 to the data.

## Changed 27079 zero count features.

## Graphing a correlation heatmap.

## Graphing a standard median correlation.

## Performing correlation.

## Graphing a distance heatmap.

## Graphing a standard median distance.

## Performing distance.

## Graphing a PCA plot.

## Graphing a T-SNE plot.

## Plotting a density plot.

## This data will benefit from being displayed on the log scale.

## If this is not desired, set scale='raw'

## Some entries are 0.  We are on log scale, setting them to 0.5.

## Changed 27079 zero count features.

## Plotting the representation of the top-n genes.

## Printing a color to condition legend.

isc_plots$tsneplot

## Pretty much all the raw-data variance is between batches.
## This might be fine, lets normalize and look again.
ixo_norm <- normalize_expt(ixo_expt, transform="log2", norm="quant", convert="cpm", filter=TRUE)

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

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

## 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: performing count filter with option: hpgl

## Removing 8908 low-count genes (11578 remaining).

## Step 2: normalizing the data with quant.

## Step 3: converting the data with cpm.

## Step 4: transforming the data with log2.

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

## Step 5: not doing batch correction.

ixo_nplots <- graph_metrics(ixo_norm)

## Graphing number of non-zero genes with respect to CPM by library.

## Graphing library sizes.

## Graphing a boxplot.

## Graphing a correlation heatmap.

## Graphing a standard median correlation.

## Performing correlation.

## Graphing a distance heatmap.

## Graphing a standard median distance.

## Performing distance.

## Graphing a PCA plot.

## Graphing a T-SNE plot.

## Plotting a density plot.

## Plotting the representation of the top-n genes.

## Printing a color to condition legend.

ixo_nplots$pcaplot

## hmm ok
ixo_nbatch <- normalize_expt(ixo_expt, transform="log2", norm="quant", convert="cpm", filter=TRUE, batch="svaseq")

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

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

## Warning in normalize_expt(ixo_expt, transform = "log2", norm = "quant", :
## Quantile normalization and sva do not always play well together.

## Step 1: performing count filter with option: hpgl

## Removing 8908 low-count genes (11578 remaining).

## Step 2: normalizing the data with quant.

## Step 3: converting the data with cpm.

## Step 4: transforming the data with log2.

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

## Step 5: doing batch correction with svaseq.

## In norm_batch, after testing logic of surrogate method/number, the
## number of surrogates is:  and the method is: be.

## 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, 10521 entries 0<x<1.

## batch_counts: Before batch correction, 29 entries are >= 0.

## After checking/setting the number of surrogates, it is: 2.

## batch_counts: Using sva::svaseq for batch correction.

## Note to self:  If you feed svaseq a data frame you will get an error like:

## data %*% (Id - mod %*% blah blah requires numeric/complex arguments.

## The number of elements which are < 0 after batch correction is: 61

## The variable low_to_zero sets whether to change <0 values to 0 and is: FALSE

ixo_nbplots <- graph_metrics(ixo_nbatch)

## Graphing number of non-zero genes with respect to CPM by library.

## Graphing library sizes.

## Graphing a boxplot.

## Graphing a correlation heatmap.

## Graphing a standard median correlation.

## Performing correlation.

## Graphing a distance heatmap.

## Graphing a standard median distance.

## Performing distance.

## Graphing a PCA plot.

## Graphing a T-SNE plot.

## Plotting a density plot.

## Plotting the representation of the top-n genes.

## Printing a color to condition legend.

ixo_nbplots$pcaplot

## Sample 613 is so ridiculously alone that it needs to go.
ixo_removed <- subset_expt(ixo_expt, subset="sampleid!='HPGL0613'")

## There were 10, now there are 9 samples.

ixo_nbatch <- normalize_expt(ixo_removed, transform="log2", norm="quant", convert="cpm", filter=TRUE, batch="svaseq")

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

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

## Warning in normalize_expt(ixo_removed, transform = "log2", norm =
## "quant", : Quantile normalization and sva do not always play well together.

## Step 1: performing count filter with option: hpgl

## Removing 9184 low-count genes (11302 remaining).

## Step 2: normalizing the data with quant.

## Step 3: converting the data with cpm.

## Step 4: transforming the data with log2.

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

## Step 5: doing batch correction with svaseq.

## In norm_batch, after testing logic of surrogate method/number, the
## number of surrogates is:  and the method is: be.

## 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, 7254 entries 0<x<1.

## batch_counts: Before batch correction, 24 entries are >= 0.

## After checking/setting the number of surrogates, it is: 1.

## batch_counts: Using sva::svaseq for batch correction.

## Note to self:  If you feed svaseq a data frame you will get an error like:

## data %*% (Id - mod %*% blah blah requires numeric/complex arguments.

## The number of elements which are < 0 after batch correction is: 69

## The variable low_to_zero sets whether to change <0 values to 0 and is: FALSE

ixo_nbplots <- graph_metrics(ixo_nbatch)

## Graphing number of non-zero genes with respect to CPM by library.

## Graphing library sizes.

## Graphing a boxplot.

## Graphing a correlation heatmap.

## Graphing a standard median correlation.

## Performing correlation.

## Graphing a distance heatmap.

## Graphing a standard median distance.

## Performing distance.

## Graphing a PCA plot.

## Graphing a T-SNE plot.

## Plotting a density plot.

## Plotting the representation of the top-n genes.

## Printing a color to condition legend.

ixo_nbplots$pcaplot

ixo_nbplots$tsneplot

ixo_filt <- normalize_expt(ixo_removed, filter=TRUE)

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

## hpgl(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

## Leaving the data in its current base format, keep in mind that
##  some metrics are easier to see when the data is log2 transformed, but
##  EdgeR/DESeq do not accept transformed data.

## Leaving the data unconverted.  It is often advisable to cpm/rpkm
##  the data to normalize for sampling differences, keep in mind though that rpkm
##  has some annoying biases, and voom() by default does a cpm (though hpgl_voom()
##  will try to detect this).

## Leaving the data unnormalized.  This is necessary for DESeq, but
##  EdgeR/limma might benefit from normalization.  Good choices include quantile,
##  size-factor, tmm, etc.

## 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: performing count filter with option: hpgl

## Removing 9184 low-count genes (11302 remaining).

## Step 2: not normalizing the data.

## Step 3: not converting the data.

## Step 4: not transforming the data.

## Step 5: not doing batch correction.

ixo_testing <- all_pairwise(ixo_filt, model_batch="sva")

## The be method chose 1 surrogate variable(s).

## Estimate type 'sva' is shorthand for 'sva_unsupervised'.

## Other sva options include: sva_supervised and svaseq.

## Attempting sva unsupervised surrogate estimation with 1 surrogates.

## Using sva to visualize before/after batch inclusion.

## Performing a test normalization with: raw

## Finished running DE analyses, collecting outputs.

## Comparing analyses 1/1: uninfected_gut_vs_infected_gut

ixo_tables <- combine_de_tables(
  ixo_testing,
  sig_excel=paste0("excel/sig_ixo_two_batches-v", ver, ".xlsx"),
  excel=paste0("excel/ixo_two_batches-v", ver, ".xlsx"))

## Deleting the file excel/ixo_two_batches-v20180119.xlsx before writing the tables.

## Writing a legend of columns.

## Printing a pca plot before/after surrogates/batch estimation.

## Working on table 1/1: uninfected_gut_vs_infected_gut

## Adding venn plots for uninfected_gut_vs_infected_gut.

## Limma expression coefficients for uninfected_gut_vs_infected_gut; R^2: 0.991; equation: y = 0.982x + 0.106

## Edger expression coefficients for uninfected_gut_vs_infected_gut; R^2: 0.992; equation: y = 1x - 0.0281

## DESeq2 expression coefficients for uninfected_gut_vs_infected_gut; R^2: 0.992; equation: y = 1x + 0.0123

## Writing summary information.

## Attempting to add the comparison plot to pairwise_summary at row: 19 and column: 1

## 
  |                                                                       
  |                                                                 |   0%
  |                                                                       
  |=================================================================| 100%

## Performing save of the workbook.

## Invoking extract_significant_genes().

## Writing a legend of columns.

## Writing excel data for uninfected_gut_vs_infected_gut: 1/4.

## After (adj)p filter, the up genes table has 48 genes.

## After (adj)p filter, the down genes table has 8 genes.

## After fold change filter, the up genes table has 42 genes.

## After fold change filter, the down genes table has 4 genes.

## Printing significant genes to the file: excel/sig_ixo_two_batches-v20180119.xlsx

## 1/1: Creating significant table up_1limma_uninfected_gut_vs_infected_gut

## The sheet name was too long for Excel, truncating it by removing vowels.
## The sheet name was too long for Excel, truncating it by removing vowels.

## Writing excel data for uninfected_gut_vs_infected_gut: 2/4.

## After (adj)p filter, the up genes table has 201 genes.

## After (adj)p filter, the down genes table has 64 genes.

## After fold change filter, the up genes table has 149 genes.

## After fold change filter, the down genes table has 12 genes.

## Printing significant genes to the file: excel/sig_ixo_two_batches-v20180119.xlsx

## 1/1: Creating significant table up_1edger_uninfected_gut_vs_infected_gut

## The sheet name was too long for Excel, truncating it by removing vowels.
## The sheet name was too long for Excel, truncating it by removing vowels.

## Writing excel data for uninfected_gut_vs_infected_gut: 3/4.

## After (adj)p filter, the up genes table has 81 genes.

## After (adj)p filter, the down genes table has 87 genes.

## After fold change filter, the up genes table has 64 genes.

## After fold change filter, the down genes table has 6 genes.

## Printing significant genes to the file: excel/sig_ixo_two_batches-v20180119.xlsx

## 1/1: Creating significant table up_1deseq_uninfected_gut_vs_infected_gut

## The sheet name was too long for Excel, truncating it by removing vowels.
## The sheet name was too long for Excel, truncating it by removing vowels.

## Writing excel data for uninfected_gut_vs_infected_gut: 4/4.

## After (adj)p filter, the up genes table has 0 genes.

## After (adj)p filter, the down genes table has 0 genes.

## After fold change filter, the up genes table has 0 genes.

## After fold change filter, the down genes table has 0 genes.

## Printing significant genes to the file: excel/sig_ixo_two_batches-v20180119.xlsx

## 1/1: Creating significant table up_1basic_uninfected_gut_vs_infected_gut

## The sheet name was too long for Excel, truncating it by removing vowels.
## The sheet name was too long for Excel, truncating it by removing vowels.

## Adding significance bar plots.

ixo_sig_shared <- write_intersect_significant(ixo_tables, extra_annot=fData(ixo_testing))

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

For the moment I will continue to assume a limited/no batch effect – but that might be wrong.

index.html

pander::pander(sessionInfo())

R version 3.4.4 (2018-03-15)

**Platform:** x86_64-pc-linux-gnu (64-bit)

locale: LC_CTYPE=en_US.utf8, LC_NUMERIC=C, LC_TIME=en_US.utf8, LC_COLLATE=en_US.utf8, LC_MONETARY=en_US.utf8, LC_MESSAGES=en_US.utf8, LC_PAPER=en_US.utf8, LC_NAME=C, LC_ADDRESS=C, LC_TELEPHONE=C, LC_MEASUREMENT=en_US.utf8 and LC_IDENTIFICATION=C

attached base packages: stats, graphics, grDevices, utils, datasets, methods and base

other attached packages: foreach(v.1.4.4), Vennerable(v.3.1.0.9000), ruv(v.0.9.7) and hpgltools(v.2018.03)

loaded via a namespace (and not attached): nlme(v.3.1-137), bitops(v.1.0-6), matrixStats(v.0.53.1), devtools(v.1.13.5), bit64(v.0.9-7), doParallel(v.1.0.11), RColorBrewer(v.1.1-2), rprojroot(v.1.3-2), GenomeInfoDb(v.1.14.0), tools(v.3.4.4), backports(v.1.1.2), R6(v.2.2.2), rpart(v.4.1-13), Hmisc(v.4.1-1), DBI(v.0.8), lazyeval(v.0.2.1), BiocGenerics(v.0.24.0), mgcv(v.1.8-23), colorspace(v.1.3-2), nnet(v.7.3-12), withr(v.2.1.2), gridExtra(v.2.3), DESeq2(v.1.18.1), bit(v.1.1-12), compiler(v.3.4.4), preprocessCore(v.1.40.0), graph(v.1.56.0), Biobase(v.2.38.0), htmlTable(v.1.11.2), xml2(v.1.2.0), DelayedArray(v.0.4.1), labeling(v.0.3), checkmate(v.1.8.5), scales(v.0.5.0.9000), DEoptimR(v.1.0-8), robustbase(v.0.92-8), genefilter(v.1.60.0), quadprog(v.1.5-5), RBGL(v.1.54.0), commonmark(v.1.4), stringr(v.1.3.0), digest(v.0.6.15), foreign(v.0.8-70), rmarkdown(v.1.9), XVector(v.0.18.0), base64enc(v.0.1-3), htmltools(v.0.3.6), limma(v.3.34.9), htmlwidgets(v.1.2), rlang(v.0.2.0.9001), rstudioapi(v.0.7), RSQLite(v.2.1.0), BiocParallel(v.1.12.0), gtools(v.3.5.0), acepack(v.1.4.1), RCurl(v.1.95-4.10), magrittr(v.1.5), GenomeInfoDbData(v.1.0.0), Formula(v.1.2-2), Matrix(v.1.2-14), Rcpp(v.0.12.16), munsell(v.0.4.3), S4Vectors(v.0.16.0), stringi(v.1.1.7), yaml(v.2.1.18), edgeR(v.3.20.9), MASS(v.7.3-49), SummarizedExperiment(v.1.8.1), zlibbioc(v.1.24.0), Rtsne(v.0.13), plyr(v.1.8.4), grid(v.3.4.4), blob(v.1.1.1), parallel(v.3.4.4), ggrepel(v.0.7.0), doSNOW(v.1.0.16), lattice(v.0.20-35), splines(v.3.4.4), pander(v.0.6.1), annotate(v.1.56.2), locfit(v.1.5-9.1), knitr(v.1.20), pillar(v.1.2.1), GenomicRanges(v.1.30.3), corpcor(v.1.6.9), geneplotter(v.1.56.0), reshape2(v.1.4.3), codetools(v.0.2-15), stats4(v.3.4.4), XML(v.3.98-1.11), evaluate(v.0.10.1), latticeExtra(v.0.6-28), data.table(v.1.10.4-3), gtable(v.0.2.0), ggplot2(v.2.2.1), openxlsx(v.4.0.17), xtable(v.1.8-2), roxygen2(v.6.0.1), survival(v.2.42-3), snow(v.0.4-2), tibble(v.1.4.2), iterators(v.1.0.9), AnnotationDbi(v.1.40.0), memoise(v.1.1.0), IRanges(v.2.12.0), cluster(v.2.0.7-1), sva(v.3.26.0) and directlabels(v.2017.03.31)

message(paste0("This is hpgltools commit: ", get_git_commit()))

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

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

## > git reset 1b009834267dea125ee94934203413fbd606e783

## R> packrat::restore()

## This is hpgltools commit: Mon Apr 23 14:59:56 2018 -0400: 1b009834267dea125ee94934203413fbd606e783

this_save <- paste0(gsub(pattern="\\.Rmd", replace="", x=rmd_file), "-v", ver, ".rda.xz")
message(paste0("Saving to ", this_save))

## Saving to 01_annotation_iscapularis-v20180119.rda.xz

tmp <- sm(saveme(filename=this_save))