CL-Brener specific material with Fernanda

Pulling together the various strands of material used in Fernanda’s manuscript. I think this should prove to just be a repetition of work already performed by Edson.

1 Collect annotation data

Actually, since this was performed with a previous release of the T.cruzi genome, it is probably smarter for me to load the annotations from my earlier savefile. Because the tritrypdb does not save the older annotation data with new releases, there is no guarantee that I will match the annotations completely if I make my own (they do provide a mapping of changed IDs which seems to work, but why waste the time messing with it when I still have the saved annotations from when I did it originally).

ver <- "20170810"
previous_file <- "01_annotation.Rmd"
tmp <- try(sm(loadme(filename=paste0(gsub(pattern="\\.Rmd", replace="", x=previous_file), "-v", ver, ".rda.xz"))))
ver <- "202005"

2 Create a CL-Brener exclusive data set excluding the problematic sample

all_expt <- sm(create_expt(metadata="sample_sheets/cl14clbr_samples_combined_all.xlsx",
                           file_column="clbrenerfile", gene_info=all_genes))
all_expt <- exclude_genes_expt(all_expt)

## Before removal, there were 25100 entries.

## Now there are 23305 entries.

## Percent kept: 99.781, 99.815, 99.814, 99.818, 97.848, 98.584, 99.776, 99.449, 99.331, 99.457, 99.321, 99.395, 99.157, 98.435, 98.543, 99.306, 99.396, 99.653, 99.663, 99.249, 99.436, 98.720, 99.886, 99.761

## Percent removed: 0.219, 0.185, 0.186, 0.182, 2.152, 1.416, 0.224, 0.551, 0.669, 0.543, 0.679, 0.605, 0.843, 1.565, 1.457, 0.694, 0.604, 0.347, 0.337, 0.751, 0.564, 1.280, 0.114, 0.239

outlier_subset <- "sampleid!='HPGL0490'"
all_minus_one <- subset_expt(all_expt, subset=outlier_subset)

## Using a subset expression.

## There were 24, now there are 23 samples.

all_minus_one <- subset_expt(all_minus_one, subset="type!='CL14'")

## Using a subset expression.

## There were 23, now there are 11 samples.

all_minus_cl14 <- subset_expt(all_expt, subset="type!='CL14'")

## Using a subset expression.

## There were 24, now there are 12 samples.

tmp_expt <- sm(normalize_expt(all_minus_cl14, filter=TRUE, convert="cpm", norm="quant", transform="log2"))
pp(file="images/all_clbr_pca.pdf", image=plot_pca(tmp_expt)$plot)

## Writing the image to: images/all_clbr_pca.pdf and calling dev.off().

tmp_expt <- sm(normalize_expt(all_minus_one, filter=TRUE, convert="cpm", norm="quant", transform="log2"))
pp(file="images/all_clbr_minus_one_pca.pdf", image=plot_pca(tmp_expt)$plot)

## Writing the image to: images/all_clbr_minus_one_pca.pdf and calling dev.off().

3 Write count tables

written <- write_expt(all_minus_one,
                      excel=glue::glue("excel/fernanda_written-v{ver}.xlsx"))

## Writing the first sheet, containing a legend and some summary data.

## Writing the raw reads.

## Graphing the raw reads.

## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'

## Warning in doTryCatch(return(expr), name, parentenv, handler): display list
## redraw incomplete

## Warning in doTryCatch(return(expr), name, parentenv, handler): display list
## redraw incomplete

## varpart sees only 1 batch, adjusting the model accordingly.

## Attempting mixed linear model with: ~  (1|condition)

## Fitting the expressionset to the model, this is slow.

## Dividing work into 100 chunks...

## 
## Total:66 s

## Placing factor: condition at the beginning of the model.

## Writing the normalized reads.

## Graphing the normalized reads.

## `geom_smooth()` using formula 'y ~ x'

## `geom_smooth()` using formula 'y ~ x'

## Warning in doTryCatch(return(expr), name, parentenv, handler): display list
## redraw incomplete

## Warning in doTryCatch(return(expr), name, parentenv, handler): display list
## redraw incomplete

## varpart sees only 1 batch, adjusting the model accordingly.

## Attempting mixed linear model with: ~  (1|condition)

## Fitting the expressionset to the model, this is slow.

## Dividing work into 100 chunks...

## 
## Total:96 s

## Placing factor: condition at the beginning of the model.

## Writing the median reads by factor.

4 Perform a DE analysis

all_minus_norm <- normalize_expt(all_minus_one, filter=TRUE)

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

## cbcb(data)

## It will save copies of each step along the way
##  in expt$normalized with the corresponding libsizes. Keep libsizes in mind
##  when invoking limma.  The appropriate libsize is 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: cbcb

## Removing 1245 low-count genes (22060 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.

all_minus_de <- all_pairwise(all_minus_norm, model_batch="svaseq")

## batch_counts: Before batch/surrogate estimation, 241564 entries are x>1: 100%.

## batch_counts: Before batch/surrogate estimation, 396 entries are x==0: 0%.

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

## Attempting svaseq estimation with 1 surrogates.

## Plotting a PCA before surrogates/batch inclusion.

## Not putting labels on the plot.

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

## Performing a test normalization with: raw

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

## log2(svaseq(cpm(cbcb(data))))

## It will save copies of each step along the way
##  in expt$normalized with the corresponding libsizes. Keep libsizes in mind
##  when invoking limma.  The appropriate libsize is 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 unnormalized.  This is necessary for DESeq, but
##  EdgeR/limma might benefit from normalization.  Good choices include quantile,
##  size-factor, tmm, etc.

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

## Removing 0 low-count genes (22060 remaining).

## Step 2: not normalizing the data.

## Step 3: converting the data with cpm.

## Step 4: transforming the data with log2.

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

## Step 5: doing batch correction with svaseq.

## Note to self:  If you get an error like 'x contains missing values' The data has too many 0's and needs a stronger low-count filter applied.

## Passing off to all_adjusters.

## batch_counts: Before batch/surrogate estimation, 237917 entries are x>1: 98%.

## batch_counts: Before batch/surrogate estimation, 396 entries are x==0: 0%.

## batch_counts: Before batch/surrogate estimation, 4347 entries are 0<x<1: 2%.

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

## Attempting svaseq estimation with 2 surrogates.

## There are 44 (0%) elements which are < 0 after batch correction.

## Setting low elements to zero.

## Not putting labels on the plot.

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

keepers <- list(
    "epi_over_t96h" = c("CLBrEpi", "CLBrA96"),
    "epi_over_t60h" = c("CLBrEpi", "CLBrA60"),
    "epi_over_tryp" = c("CLBrEpi", "CLBrTryp"))
all_minus_tables <- combine_de_tables(
    all_minus_de, keepers=keepers,
    excel=glue::glue("excel/clbr_analyses-v{ver}.xlsx"))

## Deleting the file excel/clbr_analyses-v202005.xlsx before writing the tables.

## Writing a legend of columns.

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

## Working on 1/3: epi_over_t96h which is: CLBrEpi/CLBrA96.

## Found table with CLBrEpi_vs_CLBrA96

## Working on 2/3: epi_over_t60h which is: CLBrEpi/CLBrA60.

## Found table with CLBrEpi_vs_CLBrA60

## Working on 3/3: epi_over_tryp which is: CLBrEpi/CLBrTryp.

## Found inverse table with CLBrTryp_vs_CLBrEpi

## Adding venn plots for epi_over_t96h.

## Limma expression coefficients for epi_over_t96h; R^2: 0.895; equation: y = 0.971x + 0.0236

## Deseq expression coefficients for epi_over_t96h; R^2: 0.894; equation: y = 0.838x + 1.33

## Edger expression coefficients for epi_over_t96h; R^2: 0.894; equation: y = 0.839x + 1.31

## Adding venn plots for epi_over_t60h.

## Limma expression coefficients for epi_over_t60h; R^2: 0.895; equation: y = 0.971x + 0.0236

## Deseq expression coefficients for epi_over_t60h; R^2: 0.894; equation: y = 0.838x + 1.33

## Edger expression coefficients for epi_over_t60h; R^2: 0.894; equation: y = 0.839x + 1.31

## Adding venn plots for epi_over_tryp.

## Limma expression coefficients for epi_over_tryp; R^2: 0.895; equation: y = 0.971x + 0.0236

## Deseq expression coefficients for epi_over_tryp; R^2: 0.894; equation: y = 0.838x + 1.33

## Edger expression coefficients for epi_over_tryp; R^2: 0.894; equation: y = 0.839x + 1.31

## Writing summary information, compare_plot is: TRUE.

## Performing save of excel/clbr_analyses-v202005.xlsx.

all_minus_sig <- extract_significant_genes(
    all_minus_tables,
    excel=glue::glue("excel/clbr_analyses_sig-v{ver}.xlsx"))

## Writing a legend of columns.

## Printing significant genes to the file: excel/clbr_analyses_sig-v202005.xlsx

## 1/3: Creating significant table up_limma_epi_over_t96h

## 2/3: Creating significant table up_limma_epi_over_t60h

## 3/3: Creating significant table up_limma_epi_over_tryp

## Printing significant genes to the file: excel/clbr_analyses_sig-v202005.xlsx

## 1/3: Creating significant table up_edger_epi_over_t96h

## 2/3: Creating significant table up_edger_epi_over_t60h

## 3/3: Creating significant table up_edger_epi_over_tryp

## Printing significant genes to the file: excel/clbr_analyses_sig-v202005.xlsx

## 1/3: Creating significant table up_deseq_epi_over_t96h

## 2/3: Creating significant table up_deseq_epi_over_t60h

## 3/3: Creating significant table up_deseq_epi_over_tryp

## Printing significant genes to the file: excel/clbr_analyses_sig-v202005.xlsx

## 1/3: Creating significant table up_ebseq_epi_over_t96h

## 2/3: Creating significant table up_ebseq_epi_over_t60h

## 3/3: Creating significant table up_ebseq_epi_over_tryp

## Printing significant genes to the file: excel/clbr_analyses_sig-v202005.xlsx

## 1/3: Creating significant table up_basic_epi_over_t96h

## 2/3: Creating significant table up_basic_epi_over_t60h

## 3/3: Creating significant table up_basic_epi_over_tryp

## Adding significance bar plots.

5 Epimastigote vs. not Epimastigote

epi_vs <- all_minus_one
epi_vs_meta <- pData(epi_vs)
epi_vs_meta[["epimastigote"]] <- "no"
epi_samples <- epi_vs_meta[["stage"]] == "Epi"
epi_vs_meta[epi_samples, "epimastigote"] <- "yes"
pData(epi_vs[["expressionset"]]) <- epi_vs_meta
epi_vs <- set_expt_conditions(epi_vs, fact="epimastigote")

epi_vs_de <- all_pairwise(epi_vs, model_batch="svaseq", filter=TRUE)

## batch_counts: Before batch/surrogate estimation, 241564 entries are x>1: 100%.

## batch_counts: Before batch/surrogate estimation, 396 entries are x==0: 0%.

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

## Attempting svaseq estimation with 1 surrogates.

## Plotting a PCA before surrogates/batch inclusion.

## Not putting labels on the plot.

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

## Performing a test normalization with: raw

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

## log2(svaseq(cpm(cbcb(data))))

## It will save copies of each step along the way
##  in expt$normalized with the corresponding libsizes. Keep libsizes in mind
##  when invoking limma.  The appropriate libsize is 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 unnormalized.  This is necessary for DESeq, but
##  EdgeR/limma might benefit from normalization.  Good choices include quantile,
##  size-factor, tmm, etc.

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

## Removing 0 low-count genes (22060 remaining).

## Step 2: not normalizing the data.

## Step 3: converting the data with cpm.

## Step 4: transforming the data with log2.

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

## Step 5: doing batch correction with svaseq.

## Note to self:  If you get an error like 'x contains missing values' The data has too many 0's and needs a stronger low-count filter applied.

## Passing off to all_adjusters.

## batch_counts: Before batch/surrogate estimation, 237917 entries are x>1: 98%.

## batch_counts: Before batch/surrogate estimation, 396 entries are x==0: 0%.

## batch_counts: Before batch/surrogate estimation, 4347 entries are 0<x<1: 2%.

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

## Attempting svaseq estimation with 2 surrogates.

## There are 7 (0%) elements which are < 0 after batch correction.

## Setting low elements to zero.

## Not putting labels on the plot.

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

epi_vs_table <- combine_de_tables(epi_vs_de)

## Writing a legend of columns.

## Working on table 1/1: yes_vs_no

epi_vs_sig <- extract_significant_genes(epi_vs_table)

## Writing a legend of columns.

## Printing significant genes to the file: excel/significant_genes.xlsx

## 1/1: Creating significant table up_limma_yes_vs_no

## Printing significant genes to the file: excel/significant_genes.xlsx

## 1/1: Creating significant table up_edger_yes_vs_no

## Printing significant genes to the file: excel/significant_genes.xlsx

## 1/1: Creating significant table up_deseq_yes_vs_no

## Printing significant genes to the file: excel/significant_genes.xlsx

## 1/1: Creating significant table up_ebseq_yes_vs_no

## Printing significant genes to the file: excel/significant_genes.xlsx

## 1/1: Creating significant table up_basic_yes_vs_no

## Adding significance bar plots.

6 Ontology searching

For these contrasts, we are doing ontology searches of a couple contrasts:

1. epimastigote vs. trypsomastigote
2. 60 hours vs. epimastigote

But first we need to get the correct ontology data frame and gene lengths. There is one aspect of this worth considering, in response to a query on the bioconductor email lists, I made a change to how I collect ontology information for the species of the eupathdb. Previously it was only the GOSLIM table, but now I include both that and the GO.db table. I do not know if this will affect these results.

Also, I am going to assume the default metric of ‘significant’ from the perspective of deseq2.

I might try something new, too, I also now have pretty good tables from reactome for the tryps.

all_go <- all_go[, c("GID", "GO")]
colnames(all_go) <- c("ID", "GO")
all_go <- as.data.frame(all_go)

colnames(all_lengths) <- c("ID", "length")

6.1 Epimastigote vs. Trypomastigote

sig_up <- all_minus_sig[["deseq"]][["ups"]][["epi_over_tryp"]]
sig_down <- all_minus_sig[["deseq"]][["downs"]][["epi_over_tryp"]]
epi_vs_tryp_up_goseq <- simple_goseq(
    sig_genes=sig_up,
    go_db=all_go,
    length_db=all_lengths,
    excel=glue::glue("excel/clbr_up_epi_vs_tryp_goseq-v{ver}.xlsx"))

## Using the row names of your table.

## Found 1132 genes out of 2404 from the sig_genes in the go_db.

## Found 2404 genes out of 2404 from the sig_genes in the length_db.

## Using manually entered categories.

## Calculating the p-values...

## 'select()' returned 1:1 mapping between keys and columns

## simple_goseq(): Calculating q-values

## simple_goseq(): Filling godata with terms, this is slow.

## Testing that go categories are defined.

## Removing undefined categories.

## Gathering synonyms.

## Gathering category definitions.

## simple_goseq(): Making pvalue plots for the ontologies.

## Writing data to: excel/clbr_up_epi_vs_tryp_goseq-v202005.xlsx.

epi_vs_tryp_down_goseq <- simple_goseq(
    sig_genes=sig_down,
    go_db=all_go,
    length_db=all_lengths,
    excel=glue::glue("excel/clbr_down_epi_vs_tryp_goseq-v{ver}.xlsx"))

## Using the row names of your table.

## Found 385 genes out of 3525 from the sig_genes in the go_db.

## Found 3525 genes out of 3525 from the sig_genes in the length_db.

## Using manually entered categories.

## Calculating the p-values...

## 'select()' returned 1:1 mapping between keys and columns

## simple_goseq(): Calculating q-values

## simple_goseq(): Filling godata with terms, this is slow.

## Testing that go categories are defined.

## Removing undefined categories.

## Gathering synonyms.

## Gathering category definitions.

## simple_goseq(): Making pvalue plots for the ontologies.

## Writing data to: excel/clbr_down_epi_vs_tryp_goseq-v202005.xlsx.

6.1.1 Plots

## up
epi_vs_tryp_up_goseq[["pvalue_plots"]][["bpp_plot_over"]]

epi_vs_tryp_up_goseq[["pvalue_plots"]][["mfp_plot_over"]]

## down
epi_vs_tryp_down_goseq[["pvalue_plots"]][["bpp_plot_over"]]

epi_vs_tryp_down_goseq[["pvalue_plots"]][["mfp_plot_over"]]

6.2 Epimastigote vs. 60 hours

sig_up <- all_minus_sig[["deseq"]][["ups"]][["epi_over_t60h"]]
sig_down <- all_minus_sig[["deseq"]][["downs"]][["epi_over_t60h"]]
epi_vs_t60_up_goseq <- simple_goseq(
    sig_genes=sig_up,
    go_db=all_go,
    length_db=all_lengths,
    excel=glue::glue("excel/clbr_up_epi_vs_t60_goseq-v{ver}.xlsx"))

## Using the row names of your table.

## Found 854 genes out of 2489 from the sig_genes in the go_db.

## Found 2489 genes out of 2489 from the sig_genes in the length_db.

## Using manually entered categories.

## Calculating the p-values...

## 'select()' returned 1:1 mapping between keys and columns

## simple_goseq(): Calculating q-values

## simple_goseq(): Filling godata with terms, this is slow.

## Testing that go categories are defined.

## Removing undefined categories.

## Gathering synonyms.

## Gathering category definitions.

## simple_goseq(): Making pvalue plots for the ontologies.

## Writing data to: excel/clbr_up_epi_vs_t60_goseq-v202005.xlsx.

epi_vs_t60_down_goseq <- simple_goseq(
    sig_genes=sig_down,
    go_db=all_go,
    length_db=all_lengths,
    excel=glue::glue("excel/clbr_down_epi_vs_t60_goseq-v{ver}.xlsx"))

## Using the row names of your table.

## Found 210 genes out of 624 from the sig_genes in the go_db.

## Found 624 genes out of 624 from the sig_genes in the length_db.

## Using manually entered categories.

## Calculating the p-values...

## 'select()' returned 1:1 mapping between keys and columns

## simple_goseq(): Calculating q-values

## simple_goseq(): Filling godata with terms, this is slow.

## Testing that go categories are defined.

## Removing undefined categories.

## Gathering synonyms.

## Gathering category definitions.

## simple_goseq(): Making pvalue plots for the ontologies.

## Writing data to: excel/clbr_down_epi_vs_t60_goseq-v202005.xlsx.

6.2.1 Plots

## up
epi_vs_t60_up_goseq[["pvalue_plots"]][["bpp_plot_over"]]

epi_vs_t60_up_goseq[["pvalue_plots"]][["mfp_plot_over"]]

## down
epi_vs_t60_down_goseq[["pvalue_plots"]][["bpp_plot_over"]]

epi_vs_t60_down_goseq[["pvalue_plots"]][["mfp_plot_over"]]

7 Epi vs everything else

sig_up <- epi_vs_sig[["deseq"]][["ups"]][["yes_vs_no"]]
sig_down <- epi_vs_sig[["deseq"]][["downs"]][["yes_vs_no"]]
epi_vs_all_up_goseq <- simple_goseq(
    sig_genes=sig_up,
    go_db=all_go,
    length_db=all_lengths,
    excel=glue::glue("excel/clbr_up_epi_vs_all_goseq-v{ver}.xlsx"))

## Using the row names of your table.

## Found 744 genes out of 1701 from the sig_genes in the go_db.

## Found 1701 genes out of 1701 from the sig_genes in the length_db.

## Using manually entered categories.

## Calculating the p-values...

## 'select()' returned 1:1 mapping between keys and columns

## simple_goseq(): Calculating q-values

## simple_goseq(): Filling godata with terms, this is slow.

## Testing that go categories are defined.

## Removing undefined categories.

## Gathering synonyms.

## Gathering category definitions.

## simple_goseq(): Making pvalue plots for the ontologies.

## Writing data to: excel/clbr_up_epi_vs_all_goseq-v202005.xlsx.

epi_vs_all_down_goseq <- simple_goseq(
    sig_genes=sig_down,
    go_db=all_go,
    length_db=all_lengths,
    excel=glue::glue("excel/clbr_down_epi_vs_all_goseq-v{ver}.xlsx"))

## Using the row names of your table.

## Found 222 genes out of 1852 from the sig_genes in the go_db.

## Found 1852 genes out of 1852 from the sig_genes in the length_db.

## Using manually entered categories.

## Calculating the p-values...

## 'select()' returned 1:1 mapping between keys and columns

## simple_goseq(): Calculating q-values

## simple_goseq(): Filling godata with terms, this is slow.

## Testing that go categories are defined.

## Removing undefined categories.

## Gathering synonyms.

## Gathering category definitions.

## simple_goseq(): Making pvalue plots for the ontologies.

## Writing data to: excel/clbr_down_epi_vs_all_goseq-v202005.xlsx.

7.0.1 Plots

## up
epi_vs_all_up_goseq[["pvalue_plots"]][["bpp_plot_over"]]

epi_vs_all_up_goseq[["pvalue_plots"]][["mfp_plot_over"]]

## down
epi_vs_all_down_goseq[["pvalue_plots"]][["bpp_plot_over"]]

epi_vs_all_down_goseq[["pvalue_plots"]][["mfp_plot_over"]]

8 A new request

The following, slightly mysterious request has appeared in my TODO list:

“NEW Stuff with Santuza” 1. DE analyses of CL Brener 1) Epi vs Trypo; 2) Epi vs ama60 2. GO of above

This is peculiar because of how unsatisfactory everyone seemed to find the epimastigote samples previously. I will therefore resurrect the expressionset containing all of the data and seek first to show why this is either a good or bad idea. Maybe it is fine, I do not recall.

The data structure which still includes the epimastigote samples was called ‘all_expt_all’ because it contained all samples for all haplotypes.

all_norm <- sm(normalize_expt(all_expt_all, filter=TRUE, transform="log2", convert="cpm",
                              norm="quant"))
all_plots <- sm(graph_metrics(all_norm))
all_plots$pcaplot

Well, the above PCA plot shows pretty clearly why we might be concerned about the CLBr Tryp samples, one of them looks very much like an A96. The CLBr Epis look ok, I think, but the plot is pretty crowded, I am going to pull only the 3 types of interest and see if it is clearer.

tmp <- subset_expt(all_expt_all, subset="type=='CLBr'&stage!='A96'")
tmp_norm <- sm(normalize_expt(tmp, filter=TRUE, transform="log2", convert="cpm", norm="quant"))
tt <- plot_pca(tmp_norm)
tt$plot

Ok, so I am pretty sure that the only screwball is hpgl0490. I am going to leave it in place for the moment even though it looks funky. I am also a bit torn about what is more appropriate: perform the differential expression analysis with all samples or just the 9 of interest to Santuza for this question. I think it is more appropriate to do all of them.

## This will take a while...
all_norm_all <- sm(normalize_expt(all_expt_all, filter=TRUE))
all_de <- all_pairwise(all_norm_all, model_batch="svaseq")

epi_tryp_60_keepers <- list(
  "clbrepi_vs_tryp" = c("CLBr.Epi", "CLBr.Tryp"),
  "clbrepi_vs_a60" = c("CLBr.Epi", "CLBr.A60"),
  "cl14epi_vs_tryp" = c("CL14.Epi", "CL14.Tryp"),
  "cl14epi_vs_a60" = c("CL14.Epi", "CL14.A60"))
peculiar_data <- combine_de_tables(
  all_de,
  keepers=epi_tryp_60_keepers,
  excel=paste0("excel/epi_tryp_60-v", ver, ".xlsx"))
peculiar_sig <- extract_significant_genes(
  peculiar_data,
  excel=paste0("excel/epi_tryp_60_sig-v", ver, ".xlsx"))

9 Perform goseq with this

all_go <- all_go[, c("GID", "GO")]
colnames(all_go) <- c("ID", "GO")
all_go <- as.data.frame(all_go)
sig <- peculiar_sig$limma$ups$CLBr.Epi_vs_CLBr.Tryp
colnames(all_lengths) <- c("ID", "length")

clbr_up_epi_vs_tryp_goseq <- simple_goseq(
    sig_genes=sig,
    go_db=all_go,
    length_db=all_lengths,
    excel=paste0("excel/clbr_up_epi_vs_tryp_goseq-v", ver, ".xlsx"))
clbr_up_epi_vs_tryp_goseq$pvalue_plots$bpp_plot_over

clbr_down_epi_vs_tryp_goseq <- simple_goseq(
  sig_genes=sig,
  go_db=all_go,
  length_db=all_lengths,
  excel=paste0("excel/clbr_down_epi_vs_tryp_goseq-v", ver, ".xlsx"))
clbr_down_epi_vs_tryp_goseq$pvalue_plots$bpp_plot_over

clbr_up_epi_vs_a60_goseq <- simple_goseq(
  sig_genes=sig,
  go_db=all_go,
  length_db=all_lengths,
  excel=paste0("excel/clbr_up_epi_vs_a60_goseq-v", ver, ".xlsx"))
clbr_up_epi_vs_a60_goseq$pvalue_plots$bpp_plot_over

clbr_down_epi_vs_a60_goseq <- simple_goseq(
  sig_genes=sig,
  go_db=all_go,
  length_db=all_lengths,
  excel=paste0("excel/clbr_down_epi_vs_a60_goseq-v", ver, ".xlsx"))
clbr_down_epi_vs_a60_goseq$pvalue_plots$bpp_plot_over

pander::pander(sessionInfo())

R version 3.6.1 (2019-07-05)

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

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

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

other attached packages: Rgraphviz(v.2.30.0), graph(v.1.64.0), SparseM(v.1.78), topGO(v.2.38.1), ruv(v.0.9.7.1), BiocParallel(v.1.20.1), variancePartition(v.1.16.1), hpgltools(v.1.0), Biobase(v.2.46.0) and BiocGenerics(v.0.32.0)

loaded via a namespace (and not attached): rappdirs(v.0.3.1), rtracklayer(v.1.46.0), R.methodsS3(v.1.8.0), tidyr(v.1.0.2), ggplot2(v.3.3.0), acepack(v.1.4.1), bit64(v.0.9-7), knitr(v.1.28), DelayedArray(v.0.12.3), R.utils(v.2.9.2), data.table(v.1.12.8), rpart(v.4.1-15), RCurl(v.1.98-1.2), doParallel(v.1.0.15), GenomicFeatures(v.1.38.2), preprocessCore(v.1.48.0), callr(v.3.4.3), cowplot(v.1.0.0), usethis(v.1.6.1), RSQLite(v.2.2.0), europepmc(v.0.3), bit(v.1.1-15.2), enrichplot(v.1.6.1), xml2(v.1.3.2), SummarizedExperiment(v.1.16.1), assertthat(v.0.2.1), viridis(v.0.5.1), xfun(v.0.13), hms(v.0.5.3), evaluate(v.0.14), DEoptimR(v.1.0-8), fansi(v.0.4.1), progress(v.1.2.2), caTools(v.1.18.0), dbplyr(v.1.4.3), igraph(v.1.2.5), DBI(v.1.1.0), geneplotter(v.1.64.0), htmlwidgets(v.1.5.1), stats4(v.3.6.1), purrr(v.0.3.4), ellipsis(v.0.3.0), dplyr(v.0.8.5), backports(v.1.1.6), annotate(v.1.64.0), biomaRt(v.2.42.1), vctrs(v.0.2.4), remotes(v.2.1.1), withr(v.2.2.0), ggforce(v.0.3.1), triebeard(v.0.3.0), robustbase(v.0.93-6), checkmate(v.2.0.0), GenomicAlignments(v.1.22.1), prettyunits(v.1.1.1), cluster(v.2.1.0), DOSE(v.3.12.0), crayon(v.1.3.4), genefilter(v.1.68.0), edgeR(v.3.28.1), pkgconfig(v.2.0.3), labeling(v.0.3), tweenr(v.1.0.1), GenomeInfoDb(v.1.22.1), nlme(v.3.1-147), pkgload(v.1.0.2), nnet(v.7.3-14), devtools(v.2.3.0), rlang(v.0.4.6), lifecycle(v.0.2.0), BiocFileCache(v.1.10.2), directlabels(v.2020.1.31), rprojroot(v.1.3-2), polyclip(v.1.10-0), matrixStats(v.0.56.0), Matrix(v.1.2-18), urltools(v.1.7.3), boot(v.1.3-25), base64enc(v.0.1-3), geneLenDataBase(v.1.22.0), ggridges(v.0.5.2), processx(v.3.4.2), png(v.0.1-7), viridisLite(v.0.3.0), bitops(v.1.0-6), R.oo(v.1.23.0), KernSmooth(v.2.23-17), pander(v.0.6.3), Biostrings(v.2.54.0), blob(v.1.2.1), stringr(v.1.4.0), qvalue(v.2.18.0), jpeg(v.0.1-8.1), gridGraphics(v.0.5-0), S4Vectors(v.0.24.4), scales(v.1.1.0), memoise(v.1.1.0), magrittr(v.1.5), plyr(v.1.8.6), gplots(v.3.0.3), gdata(v.2.18.0), zlibbioc(v.1.32.0), compiler(v.3.6.1), RColorBrewer(v.1.1-2), lme4(v.1.1-23), DESeq2(v.1.26.0), Rsamtools(v.2.2.3), cli(v.2.0.2), XVector(v.0.26.0), ps(v.1.3.2), htmlTable(v.1.13.3), Formula(v.1.2-3), MASS(v.7.3-51.6), mgcv(v.1.8-31), tidyselect(v.1.0.0), stringi(v.1.4.6), yaml(v.2.2.1), GOSemSim(v.2.12.1), askpass(v.1.1), locfit(v.1.5-9.4), latticeExtra(v.0.6-29), ggrepel(v.0.8.2), fastmatch(v.1.1-0), tools(v.3.6.1), rstudioapi(v.0.11), foreach(v.1.5.0), foreign(v.0.8-76), gridExtra(v.2.3), farver(v.2.0.3), Rtsne(v.0.15), ggraph(v.2.0.2), digest(v.0.6.25), rvcheck(v.0.1.8), BiocManager(v.1.30.10), quadprog(v.1.5-8), Rcpp(v.1.0.4.6), GenomicRanges(v.1.38.0), httr(v.1.4.1), AnnotationDbi(v.1.48.0), colorspace(v.1.4-1), XML(v.3.99-0.3), fs(v.1.4.1), IRanges(v.2.20.2), splines(v.3.6.1), RBGL(v.1.62.1), statmod(v.1.4.34), graphlayouts(v.0.7.0), ggplotify(v.0.0.5), sessioninfo(v.1.1.1), xtable(v.1.8-4), jsonlite(v.1.6.1), nloptr(v.1.2.2.1), tidygraph(v.1.1.2), corpcor(v.1.6.9), testthat(v.2.3.2), R6(v.2.4.1), Vennerable(v.3.1.0.9000), Hmisc(v.4.4-0), pillar(v.1.4.3), htmltools(v.0.4.0), glue(v.1.4.0), minqa(v.1.2.4), clusterProfiler(v.3.14.3), codetools(v.0.2-16), fgsea(v.1.12.0), pkgbuild(v.1.0.7), lattice(v.0.20-41), tibble(v.3.0.1), sva(v.3.34.0), pbkrtest(v.0.4-8.6), curl(v.4.3), BiasedUrn(v.1.07), colorRamps(v.2.3), gtools(v.3.8.2), zip(v.2.0.4), GO.db(v.3.10.0), openxlsx(v.4.1.4), openssl(v.1.4.1), survival(v.3.1-12), limma(v.3.42.2), rmarkdown(v.2.1), desc(v.1.2.0), munsell(v.0.5.0), DO.db(v.2.9), fastcluster(v.1.1.25), GenomeInfoDbData(v.1.2.2), iterators(v.1.0.12), goseq(v.1.38.0), reshape2(v.1.4.4) and gtable(v.0.3.0)

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 5763684c1dc9aaca47da1768d30388741268c02c

## This is hpgltools commit: Sat Apr 25 14:42:57 2020 -0400: 5763684c1dc9aaca47da1768d30388741268c02c

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

## Saving to 01_annotation-v202005.rda.xz

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