M. musculus 3 cell types, 1 timepoint, 3 genotypes, and 1 replicate.

1 M. musculus

This will be a very minimal analysis until we get some replicates.

1.1 Annotations

I am using mm38_95.

## My load_biomart_annotations() function defaults to human, so that will be quick.
mm_annot <- load_biomart_annotations(species="mmusculus")

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

mm_annot <- mm_annot[["annotation"]]
mm_annot[["txid"]] <- paste0(mm_annot[["ensembl_transcript_id"]], ".", mm_annot[["version"]])
rownames(mm_annot) <- make.names(mm_annot[["ensembl_gene_id"]], unique=TRUE)

tx_gene_map <- mm_annot[, c("txid", "ensembl_gene_id")]

So, I now have 2 data frames of parasite annotations and 1 human.

1.2 Metadata

I am going to write a quick sample sheet in the current working directory called ‘all_samples.xlsx’ and put the names of the count tables in it.

1.3 Create expressionsets

Here I combine the metadata, count data, and annotations.

It is worth noting that the gene IDs from htseq-count probably do not match the annotations retrieved because they are likely exon-based rather than gene based. This is not really a problem, but don’t forget it!

mm38_salmon <- sm(create_expt("sample_sheets/all_samples.xlsx", tx_gene_map=tx_gene_map,
                              gene_info=mm_annot, file_column="salmonfile"))

mmtx_annot <- mm_annot
rownames(mmtx_annot) <- mm_annot[["txid"]]
mm38_saltx <- sm(create_expt("sample_sheets/all_samples.xlsx",
                             gene_info=mmtx_annot, file_column="salmonfile"))

hisat_annot <- mm_annot
##mm38_hisat <- create_expt("sample_sheets/all_samples.xlsx",
##                          gene_info=hisat_annot)

1.4 Query expressionsets

In this block I will calculate all the diagnostic plots, but not show them. I will show them next with a little annotation.

I will leave the output for the first of each invocation and silence it for the second.

mm38_plots <- sm(graph_metrics(mm38_salmon))

mm38_norm <- normalize_expt(mm38_salmon, norm="quant", convert="cpm",
                            transform="log2", filter=TRUE)

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

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

## 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 2794 low-count genes (3970 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 1105 values equal to 0, adding 1 to the matrix.

## Step 5: not doing batch correction.

mm38n_plots <- sm(graph_metrics(mm38_norm))

1.5 Show some plots

mm38_plots$legend

mm38_plots$libsize

mm38_plots$nonzero

mm38n_plots$density

mm38n_plots$pc_plot

1.6 Do a simple DE

The only interesting DE I see in this is to compare the retinas to the dlgns. I can treat them as replicates and compare.

mm <- set_expt_conditions(mm38_salmon, fact="celltype")
mm_norm <- sm(normalize_expt(mm, norm="quant",
                             convert="cpm", transform="log2", filter=TRUE))
plot_pca(mm_norm)$plot

mm_de <- all_pairwise(mm, model_batch=FALSE)

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

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

## 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 2794 low-count genes (3970 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 1105 values equal to 0, adding 1 to the matrix.

## Step 5: not doing batch correction.

## Plotting a PCA before surrogates/batch inclusion.

## Not putting labels on the plot.

## Assuming no batch in model for testing pca.

## Not putting labels on the plot.

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

1.7 Set up for initial analysis

Since we don’t have replicates, I am going to do simple subtractions. The most important caveat is that I must subtract each baseline from its ht/ko samples before considering the contrasts of interest.

dlgnbaseline <- "iprgc_01"
retbaseline <- "iprgc_02"
scnbaseline <- "iprgc_03"

dlgnnorm <- "iprgc_04"
retnorm <- "iprgc_05"
scnnorm <- NULL  ## Does not yet exist.

dlgnko <- "iprgc_06"
retko <- "iprgc_07"
scnko <- "iprgc_08"

normret <- exprs(mm38_norm)[, retnorm] - exprs(mm38_norm)[, retbaseline]
koret <- exprs(mm38_norm)[, retko] - exprs(mm38_norm)[, retbaseline]
koscn <- exprs(mm38_norm)[, scnko] - exprs(mm38_norm)[, scnbaseline]
normdlgn <- exprs(mm38_norm)[, dlgnnorm] - exprs(mm38_norm)[, dlgnbaseline]
kodlgn <- exprs(mm38_norm)[, dlgnko] - exprs(mm38_norm)[, dlgnbaseline]

1.8 Contrasts of interest

Now that we have the 5 available ‘conditions’, perform the subtractions to query the data.

## het vs het
normret_vs_normdlgn <- normret - normdlgn

## ko vs ko
koret_vs_kodlgn <- koret - kodlgn
koret_vs_koscn <- koret - koscn
kodlgn_vs_koscn <- kodlgn - koscn

## ko vs het
koret_vs_normret <- koret - normret
kodlgn_vs_normdlgn <- kodlgn - normdlgn
koscn_vs_normdlgn <- koscn - normdlgn

## ratio of ratios
normko_retdlgn <- normret_vs_normdlgn - koret_vs_kodlgn

pair_mtrx <- cbind(
  ## Baseline subtractions
  normdlgn, normret, kodlgn, koret, koscn,
  ## het_vs_het, of which there is only 1 because we do not have hetscn
  normret_vs_normdlgn,
  ## ko_vs_ko, of which we have 3
  koret_vs_kodlgn, koret_vs_koscn, kodlgn_vs_koscn,
  ## ko_vs_het, 3 including one getting around missing hetscn
  koret_vs_normret, kodlgn_vs_normdlgn, koscn_vs_normdlgn,
  ## ratio of ratios
  normko_retdlgn)

mm_tables <- sm(combine_de_tables(
  mm_de, extra_annot=pair_mtrx,
  excel=glue::glue("excel/{rundate}mm_tables-v{ver}.xlsx")))

2 Some pictures

As I understand it, there is some interest in an ontology search using the ratio of ratios.

ror <- normko_retdlgn
up_idx <- ror >= 1
down_idx <- ror <= -1
ror_up <- ror[up_idx]
length(ror_up)

## [1] 599

ror_down <- ror[down_idx]
length(ror_down)

## [1] 550

ror_gprofiler_up <- simple_gprofiler(sig_genes=ror_up, species="mmusculus",
                                     excel=glue::glue("excel/{rundate}mm_ror_gpfoiler_up-v{ver}.xlsx"))

## Performing gProfiler GO search of 599 genes against mmusculus.

## GO search found 95 hits.

## Performing gProfiler KEGG search of 599 genes against mmusculus.

## KEGG search found 2 hits.

## Performing gProfiler REAC search of 599 genes against mmusculus.

## REAC search found 0 hits.

## Performing gProfiler MI search of 599 genes against mmusculus.

## MI search found 0 hits.

## Performing gProfiler TF search of 599 genes against mmusculus.

## TF search found 38 hits.

## Performing gProfiler CORUM search of 599 genes against mmusculus.

## CORUM search found 0 hits.

## Performing gProfiler HP search of 599 genes against mmusculus.

## HP search found 4 hits.

## Writing data to: excel/20200110mm_ror_gpfoiler_up-v20200110.xlsx.

## Finished writing data.

ror_gprofiler_up$pvalue_plots$mfp_plot_over

ror_gprofiler_up$pvalue_plots$bpp_plot_over

ror_gprofiler_up$pvalue_plots$ccp_plot_over

ror_gprofiler_up$pvalue_plots$tf_plot_over

ror_gprofiler_up$pvalue_plots$hp_plot_over

ror_gprofiler_down <- simple_gprofiler(sig_genes=ror_down, species="mmusculus",
                                       excel=glue::glue("excel/{rundate}mm_ror_gpfoiler_down-v{ver}.xlsx"))

## Performing gProfiler GO search of 550 genes against mmusculus.

## GO search found 243 hits.

## Performing gProfiler KEGG search of 550 genes against mmusculus.

## KEGG search found 12 hits.

## Performing gProfiler REAC search of 550 genes against mmusculus.

## REAC search found 5 hits.

## Performing gProfiler MI search of 550 genes against mmusculus.

## MI search found 0 hits.

## Performing gProfiler TF search of 550 genes against mmusculus.

## TF search found 37 hits.

## Performing gProfiler CORUM search of 550 genes against mmusculus.

## CORUM search found 0 hits.

## Performing gProfiler HP search of 550 genes against mmusculus.

## HP search found 0 hits.

## Writing data to: excel/20200110mm_ror_gpfoiler_down-v20200110.xlsx.

## Finished writing data.

ror_gprofiler_down$pvalue_plots$mfp_plot_over

ror_gprofiler_down$pvalue_plots$bpp_plot_over

ror_gprofiler_down$pvalue_plots$reactome_plot_over

ror_gprofiler_down$pvalue_plots$ccp_plot_over

ror_gprofiler_down$pvalue_plots$tf_plot_over

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))

loadme(filename=this_save)