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
rownames(hisat_annot) <- paste0("gene.", rownames(hisat_annot))
mm38_hisat <- create_expt("sample_sheets/all_samples.xlsx",
                          gene_info=hisat_annot)

## Reading the sample metadata.

## The sample definitions comprises: 8 rows(samples) and 8 columns(metadata fields).

## Reading count tables.

## Reading count tables with read.table().

## /mnt/sshfs/cbcbsub/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/mmusculus_iprgc_2019/preprocessing/iprgc_01/outputs/hisat2_mm38_95/r1_trimmed.count.xz contains 25788 rows.

## /mnt/sshfs/cbcbsub/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/mmusculus_iprgc_2019/preprocessing/iprgc_02/outputs/hisat2_mm38_95/r1_trimmed.count.xz contains 25788 rows and merges to 25788 rows.

## /mnt/sshfs/cbcbsub/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/mmusculus_iprgc_2019/preprocessing/iprgc_03/outputs/hisat2_mm38_95/r1_trimmed.count.xz contains 25788 rows and merges to 25788 rows.

## /mnt/sshfs/cbcbsub/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/mmusculus_iprgc_2019/preprocessing/iprgc_04/outputs/hisat2_mm38_95/r1_trimmed.count.xz contains 25788 rows and merges to 25788 rows.

## /mnt/sshfs/cbcbsub/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/mmusculus_iprgc_2019/preprocessing/iprgc_05/outputs/hisat2_mm38_95/r1_trimmed.count.xz contains 25788 rows and merges to 25788 rows.

## /mnt/sshfs/cbcbsub/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/mmusculus_iprgc_2019/preprocessing/iprgc_06/outputs/hisat2_mm38_95/r1_trimmed.count.xz contains 25788 rows and merges to 25788 rows.

## /mnt/sshfs/cbcbsub/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/mmusculus_iprgc_2019/preprocessing/iprgc_07/outputs/hisat2_mm38_95/r1_trimmed.count.xz contains 25788 rows and merges to 25788 rows.

## /mnt/sshfs/cbcbsub/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/mmusculus_iprgc_2019/preprocessing/iprgc_08/outputs/hisat2_mm38_95/r1_trimmed.count.xz contains 25788 rows and merges to 25788 rows.

## Finished reading count tables.

## Matched 25554 annotations and counts.

## Bringing together the count matrix and gene information.

## Some annotations were lost in merging, setting them to 'undefined'.

## The final expressionset has 25783 rows and 8 columns.

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.

1.4.1 Initial salmon plots

mm38_plots_sa <- sm(graph_metrics(mm38_salmon))

mm38_norm_sa <- 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_sa <- sm(graph_metrics(mm38_norm_sa))

mm38_plots_sa$legend

mm38_plots_sa$libsize

mm38_plots_sa$nonzero

mm38n_plots_sa$density

mm38n_plots_sa$pc_plot

1.4.2 Initial hisat2 plots

mm38_plots_hi <- sm(graph_metrics(mm38_hisat))

mm38_norm_hi <- normalize_expt(mm38_hisat, 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 12233 low-count genes (13550 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 19 values equal to 0, adding 1 to the matrix.

## Step 5: not doing batch correction.

mm38n_plots_hi <- sm(graph_metrics(mm38_norm_hi))

mm38_plots_hi$libsize

mm38_plots_hi$nonzero

mm38n_plots_hi$density

mm38n_plots_hi$pc_plot

1.5 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.

These differential expression analyses are EXPLICITLY NOT what you care about. However, they are useful for two purposes:

1. Seeing that the three tissue types are indeed different.
2. Setting up the table of results with appropriate rows/columns of (rows)genes and (columns) annotations. We will therefore add to these tables the results of the expression analyses that you actually do care about.

When we receive full replicate sets, this cheater method of encapsulating the data will not longer be required.

1.5.1 With salmon

mm_sa <- set_expt_conditions(mm38_salmon, fact="celltype")
mm_norm_sa <- sm(normalize_expt(mm_sa, convert="rpkm", transform="log2", column="cds_length"))
plot_pca(mm_norm_sa)$plot

mm_de_sa <- all_pairwise(mm_sa, 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.5.2 With hisat2

mm_hi <- set_expt_conditions(mm38_hisat, fact="celltype")
mm_norm_hi <- sm(normalize_expt(mm_hi, convert="rpkm", transform="log2", column="cds_length"))
plot_pca(mm_norm_hi)$plot

mm_de_hi <- sm(all_pairwise(mm_hi, model_batch=FALSE))

1.6 Set up for initial analysis

Until we get full replicates, I will do simple subtractions.

1.6.1 Term definition

In an attempt to keep some clarity in the terms used, I want to define them now. There are three contexts in which we will consider the data:

1. The individual sample type. When considering individual samples, I will use three terms in this and only this context: wild-type (wt), het, and mut.

2. The individual translatome. These are defines as something / baseline. I will exclusively call the wt samples ‘baseline’ when speaking in this context. I will exclusively state ‘normal’ when referring to het / wt samples, and I will state ‘ko’ when referring to mut / wt samples in the translatome context.

3. Translatome vs. translatome. Whenever comparing translatomes, I will use the names as in #2 and always put the numerator first when writing the name of a comparison.

The most complex example of the above nomenclature is:

“normko_retdlgn is defined as normret_vs_normdlgn - koret_vs_kodlgn”

This states we are examining at the translatome context: (norm(retina translatome) - norm(dlgn translatome)) - (ko(retina translatome) - ko(dlgn translatome))

Which in turn is synonymous to the following at the sample context: ((rethet - retwt) - (dlgnhet - dlgnwt)) - ((retko - retwt) - (dlgnko - dlgnwt))

Now let us associate the various variable names with the appropriate samples:

dlgnwt <- "iprgc_01"
retwt <- "iprgc_02"
scnwt <- "iprgc_03"

dlgnhet <- "iprgc_04"
rethet <- "iprgc_05"
scnhet <- NULL  ## Does not yet exist.

dlgnmut <- "iprgc_06"
retmut <- "iprgc_07"
scnmut <- "iprgc_08"

Give these variable names, now lets associate columns of the expression data with them. These are at the sample context, so the appropriate names are: ‘wt’, ‘het’, and ‘mut’. In each case I will prefix the genotype with the tissue type: ‘ret’, ‘dlgn’, and ‘scn’. Thus ‘retwt’ refers to the sample used to calculate the translatome retina baseline; in contrast ‘dlgnmut’ is the sample which provides the dlgn knockout.

## Sample context
mm38_norm <- mm_norm_sa
dlgnwt <- exprs(mm38_norm)[, dlgnwt]
retwt <- exprs(mm38_norm)[, retwt]
scnwt <- exprs(mm38_norm)[, scnwt]
dlgnhet <- exprs(mm38_norm)[, dlgnhet]
rethet <- exprs(mm38_norm)[, rethet]
dlgnmut <- exprs(mm38_norm)[, dlgnmut]
retmut <- exprs(mm38_norm)[, retmut]
scnmut <- exprs(mm38_norm)[, scnmut]

Each of the above 8 variables provides 1 column of information. We have 3 baseline comparisons available to us. In each of these we compare one wt sample to another.

## Baseline comparisons
wt_dlgnret <- dlgnwt - retwt
wt_scnret <- scnwt - retwt
wt_dlgnscn <- dlgnwt - scnwt

Simultaneously, we have 5 available translatomes. This are provided by comparing each het or mut to the associated wt. These will therefore receive names: ‘norm’ and ‘ko’ instead of ‘het’ and ‘mut’.

## Translatome context
normret <- rethet - retwt
koret <- retmut - retwt
koscn <- scnmut - scnwt
normdlgn <- dlgnhet - dlgnwt
kodlgn <- dlgnmut - dlgnwt

Given these translatomes, there are a few contrasts of likely interest. These are performed by comparing the relevant translatomes.

Will will split these into 4 separate categories: het vs het, ko vs ko, ko vs het, and ratio vs ratio.

Finally, note that we are being explicitly redundant in these definitions. I am making variable names for both the a/b ratio and the b/a ratio. Thus we have some redundantly redundant (haha) flexibility when deciding on what we want to plot.

## norm vs norm
normdlgn_vs_normret <- normdlgn - normret
normret_vs_normdlgn <- normret - normdlgn

## ko vs ko
koret_vs_kodlgn <- koret - kodlgn
kodlgn_vs_koret <- kodlgn - koret

koret_vs_koscn <- koret - koscn
koscn_vs_koret <- koscn - koret

kodlgn_vs_koscn <- kodlgn - koscn
koscn_vs_kodlgn <- koscn - kodlgn

On the other hand, I am assuming we always want the normals as denominators and kos as numerators.

## ko vs norm
koret_vs_normret <- koret - normret

kodlgn_vs_normdlgn <- kodlgn - normdlgn

Finally, here is the ratio of ratios example I printed above:

I named it ‘normko_retdlgn’ in an attempt to make clear that it is actually: (normret/normdlgn)/(koret/kodlgn)

or stated differently: “norm divided by ko for ret divided by dlgn.”

## ratio of ratios
normko_retdlgn <- normret_vs_normdlgn - koret_vs_kodlgn

1.7 Define a matrix of these values.

My matrix of data will now contain 1 column for each of the above 27 samples/comparisons.

pair_mtrx <- cbind(
  ## Individual samples
  dlgnwt, retwt, scnwt, dlgnhet, rethet, dlgnmut, retmut, scnmut,
  ## Baseline comparisons
  wt_dlgnret, wt_scnret, wt_dlgnscn,
  ## Baseline subtractions
  normdlgn, normret, kodlgn, koret, koscn,
  ## het_vs_het, of which there is only 1 because we do not have hetscn
  normdlgn_vs_normret, normret_vs_normdlgn,
  ## ko_vs_ko, of which we have 3
  koret_vs_kodlgn, kodlgn_vs_koret,
  koret_vs_koscn, koscn_vs_koret,
  kodlgn_vs_koscn, koscn_vs_kodlgn,
  ## ko_vs_het, 3 including one getting around missing hetscn
  koret_vs_normret, kodlgn_vs_normdlgn,
  ## ratio of ratios
  normko_retdlgn)

1.8 Cutoffs

I am not sure if we will use these indexes, but I am writing these out as subsets of genes to look at. These indexes are stating that, given a cutoff (0), we want to look at only the genes which have higher x / baseline values than the cutoff.

## Queries about gene subsets.
## These are all in the context of translatomes.
cutoff <- 0
ret_kept_idx <- normret > cutoff & koret > cutoff
scn_kept_idx <- koscn > cutoff
dlgn_kept_idx <- normdlgn > cutoff & kodlgn > cutoff
ret_dlgn_kept_idx <- ret_kept_idx & dlgn_kept_idx
ret_scn_kept_idx <- ret_kept_idx & scn_kept_idx
dlgn_scn_kept_idx <- dlgn_kept_idx & scn_kept_idx

##normdlgn_vs_normret[!ret_dlgn_kept_idx] <- NA
##normret_vs_normdlgn[!ret_dlgn_kept_idx] <- NA
##koret_vs_kodlgn[!ret_dlgn_kept_idx] <- NA
##kodlgn_vs_koret[!ret_dlgn_kept_idx] <- NA
##koret_vs_koscn[!ret_scn_kept_idx] <- NA
##koscn_vs_koret[!ret_scn_kept_idx] <- NA
##kodlgn_vs_koscn[!dlgn_scn_kept_idx] <- NA
##koscn_vs_kodlgn[!dlgn_scn_kept_idx] <- NA
##koret_vs_normret[!ret_kept_idx] <- NA
##kodlgn_vs_normdlgn[!dlgn_kept_idx] <- NA
##normko_retdlgn <- normko_retdlgn[!ret_dlgn_kept_idx] <- NA

1.9 Add the matrix to the differential expression

I will use my function combine_de_tables() to add this information to my existing annotation data along with the results from the statistically valid comparison of the three tissue types.

{r{r add_matrix_de, fig.show=“hide”} mm_tables <- sm(combine_de_tables( mm_de_sa, extra_annot=pair_mtrx, excel=glue::glue(“excel/{rundate}mm_salmon_tables-v{ver}.xlsx”))) ```

2 Plots of interesting comparisons

## Put retina baseline on y axis as black, retina het on x axis as black.
## Then recolor a subset of these as red, the reds are when normret > 0
library(ggplot2)
plotted <- as.data.frame(pair_mtrx[, c("rethet", "retwt")])
red_idx <- normret > 0
plotted[, "color"] <- ifelse(red_idx, "red", "black")
ret_hetwt <- ggplot2::ggplot(plotted, aes_string(x="rethet",
                                                 y="retwt",
                                                 color="color")) +
  geom_point(alpha=0.5) +
  scale_color_manual(values=c("black", "red"))
ret_hetwt

## Warning: Removed 37 rows containing missing values (geom_point).

plotted <- pair_mtrx[, c("retmut", "retwt")]
ret_mutwt <- plot_scatter(plotted)
ret_mutwt

plotted <- pair_mtrx[, c("dlgnhet", "dlgnwt")]
dlgn_hetwt <- plot_scatter(plotted)
dlgn_hetwt

plotted <- pair_mtrx[, c("dlgnmut", "dlgnwt")]
dlgn_mutwt <- plot_scatter(plotted)
dlgn_mutwt

plotted <- pair_mtrx[, c("scnmut", "scnwt")]
scn_mutwt <- plot_scatter(plotted)
scn_mutwt

## Axon translatome specific
##  x-axis: normdlgn_vs_normret or normret_vs_normdlgn,
##              ^^^^
##  y-axis: dlgnwt-retwt (baseline dlgn - baseline retina)
plotted <- as.data.frame(pair_mtrx[, c("normdlgn_vs_normret", "wt_dlgnret")])
red_idx <- normret > 0
## Note that this order is opposite of above.
plotted[, "color"] <- ifelse(red_idx, "black", "red")
axon_trans_ret_target <- ggplot2::ggplot(plotted, aes_string(x="normdlgn_vs_normret",
                                                             y="wt_dlgnret",
                                                             color="color")) +
  geom_point(alpha=0.5) +
  scale_color_manual(values=c("black", "red"))
axon_trans_ret_target

## Warning: Removed 37 rows containing missing values (geom_point).

## DLGN translatome wrt. Retina translatome
plotted <- pair_mtrx[, c("normret", "normdlgn")]
normret_normdlgn <- plot_scatter(plotted)
normret_normdlgn

plotted <- pair_mtrx[, c("koret", "kodlgn")]
koret_kodlgn <- plot_scatter(plotted)
koret_kodlgn

plotted <- pair_mtrx[, c("koret", "koscn")]
koret_koscn_plot <- plot_scatter(plotted)
koret_koscn_plot

plotted <- pair_mtrx[, c("normdlgn", "kodlgn")]
normdlgn_kodlgn_plot <- plot_scatter(plotted)
normdlgn_kodlgn_plot

plotted <- pair_mtrx[, c("normret", "koret")]
normret_koret_plot <- plot_scatter(plotted)
normret_koret_plot

plotted <- pair_mtrx[, c("normret_vs_normdlgn", "koret_vs_kodlgn")]
normal_ko_axon_translatome <- plot_scatter(plotted)
normal_ko_axon_translatome

3 Translatome level comparisons

3.1 Make a generic plotter for this stuff.

translatome_plotter <- function(pair_mtrx, x_axis="koret", y_axis="koscn", lfc=NULL, fc=NULL, linewidth=1.5,
                                up_color="red", down_color="#098534", line_color="#fcba03", alpha=0.5,
                                x_limit=c(-7, 7), y_limit=c(-7, 7)) {
  if (is.null(fc) & is.null(lfc)) {
    message("No fc/lfc was provided, defaulting to 10 fold.")
    lfc <- log2(10)
  } else if (is.null(lfc)) {
    lfc <- log2(fc)
  }
  plotted <- as.data.frame(pair_mtrx[, c(x_axis, y_axis)])
  na_idx <- is.na(plotted)
  plotted[na_idx] <- 0
  up_idx <- plotted[, y_axis] - plotted[, x_axis] >= lfc
  down_idx <- plotted[, y_axis] - plotted[, x_axis] <= (lfc * -1)
  up_genes <- rownames(plotted)[up_idx]
  down_genes <- rownames(plotted)[down_idx]
  ## Note that this order is opposite of above.
  plotted[["color"]] <- "black"
  plotted[up_idx, "color"] <- up_color
  plotted[down_idx, "color"] <- down_color
  plotted[["color"]] <- as.factor(plotted[["color"]])
  levels(plotted[["color"]]) <- c("black", up_color, down_color)
  plt <- ggplot2::ggplot(plotted, aes_string(x=x_axis,
                                             y=y_axis,
                                             color="color")) +
    geom_abline(size=1.1, slope=1, intercept=lfc, color="orange") +
    geom_abline(size=1.1, slope=1, intercept=(-1 * lfc), color="orange") +
    scale_x_continuous(limits=c(x_limit)) +
    scale_y_continuous(limits=c(y_limit)) +
    geom_point(alpha=alpha) +
    scale_color_manual(values=c(down_color, "black", up_color))
  retlist <- list(
    "mtrx" = plotted,
    "ups" = up_genes,
    "downs" = down_genes,
    "plot" = plt)
  return(retlist)
}

1. X-axis: Always retina.
2. Y-axis: Always a target tissue.

First plot: KO scn translatome on y axis vs. KO retina translatome on x axis.

## Scn knockout translatome vs retina knockout translatome.
scnko_wrt_retko_translatome <- translatome_plotter(pair_mtrx)

## No fc/lfc was provided, defaulting to 10 fold.

scnko_wrt_retko_translatome$plot

## Warning: Removed 1 rows containing missing values (geom_point).

up_go <- simple_gprofiler(sig_genes=scnko_wrt_retko_translatome$ups, species="mmusculus")

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

## GO search found 0 hits.

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

## KEGG search found 0 hits.

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

## REAC search found 0 hits.

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

## MI search found 0 hits.

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

## TF search found 0 hits.

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

## CORUM search found 0 hits.

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

## HP search found 0 hits.

down_go <- simple_gprofiler(sig_genes=scnko_wrt_retko_translatome$downs, species="mmusculus")

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

## GO search found 2 hits.

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

## KEGG search found 1 hits.

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

## REAC search found 0 hits.

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

## MI search found 0 hits.

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

## TF search found 10 hits.

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

## CORUM search found 1 hits.

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

## HP search found 0 hits.

## dlgn normal translatome vs retina normal translatome.
dlgnnorm_wrt_retnorm_translatome <- translatome_plotter(pair_mtrx,
                                                        x_axis="normret", y_axis="normdlgn")

## No fc/lfc was provided, defaulting to 10 fold.

dlgnnorm_wrt_retnorm_translatome$plot

## Warning: Removed 1 rows containing missing values (geom_point).

up_go <- simple_gprofiler(sig_genes=dlgnnorm_wrt_retnorm_translatome$ups, species="mmusculus")

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

## GO search found 3 hits.

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

## KEGG search found 0 hits.

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

## REAC search found 2 hits.

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

## MI search found 0 hits.

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

## TF search found 0 hits.

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

## CORUM search found 1 hits.

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

## HP search found 0 hits.

down_go <- simple_gprofiler(sig_genes=dlgnnorm_wrt_retnorm_translatome$downs, species="mmusculus")

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

## GO search found 14 hits.

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

## KEGG search found 1 hits.

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

## REAC search found 3 hits.

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

## MI search found 0 hits.

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

## TF search found 4 hits.

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

## CORUM search found 1 hits.

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

## HP search found 47 hits.

down_go$go

##    query.number significant  p.value term.size query.size overlap.size
## 1             1        TRUE 1.39e-02       115         23            4
## 2             1        TRUE 1.34e-02       114         23            4
## 3             1        TRUE 3.77e-04        47         23            4
## 4             1        TRUE 1.67e-02        36         23            3
## 5             1        TRUE 2.11e-02       366         45            7
## 6             1        TRUE 1.97e-02       362         45            7
## 7             1        TRUE 1.97e-02       362         45            7
## 8             1        TRUE 2.65e-08       142         45            9
## 9             1        TRUE 2.04e-08       138         45            9
## 10            1        TRUE 1.32e-03       148         45            6
## 11            1        TRUE 1.39e-02       126         45            5
## 12            1        TRUE 3.64e-03        96         45            5
## 13            1        TRUE 3.45e-03        95         45            5
## 14            1        TRUE 8.05e-04        71         45            5
##    precision recall    term.id domain subgraph.number
## 1      0.174  0.035 GO:0009581     BP               2
## 2      0.174  0.035 GO:0009582     BP               2
## 3      0.174  0.085 GO:0009583     BP               2
## 4      0.130  0.083 GO:0009584     BP               2
## 5      0.156  0.019 GO:0048880     BP               4
## 6      0.156  0.019 GO:0150063     BP               4
## 7      0.156  0.019 GO:0001654     BP               4
## 8      0.200  0.063 GO:0050953     BP               5
## 9      0.200  0.065 GO:0007601     BP               5
## 10     0.133  0.041 GO:0060041     BP               3
## 11     0.111  0.040 GO:0097730     CC               1
## 12     0.111  0.052 GO:0097731     CC               1
## 13     0.111  0.053 GO:0097733     CC               1
## 14     0.111  0.070 GO:0001750     CC               1
##                                term.name relative.depth
## 1         detection of external stimulus              1
## 2          detection of abiotic stimulus              1
## 3            detection of light stimulus              2
## 4             detection of visible light              3
## 5             sensory system development              1
## 6              visual system development              2
## 7                        eye development              1
## 8   sensory perception of light stimulus              1
## 9                      visual perception              2
## 10 retina development in camera-type eye              1
## 11                     non-motile cilium              1
## 12                 9+0 non-motile cilium              2
## 13             photoreceptor cell cilium              3
## 14           photoreceptor outer segment              4
##                                                                                                                                                                  intersection
## 1                                                                                                 ENSMUSG00000021804,ENSMUSG00000029491,ENSMUSG00000031293,ENSMUSG00000037446
## 2                                                                                                 ENSMUSG00000021804,ENSMUSG00000029491,ENSMUSG00000031293,ENSMUSG00000037446
## 3                                                                                                 ENSMUSG00000021804,ENSMUSG00000029491,ENSMUSG00000031293,ENSMUSG00000037446
## 4                                                                                                                    ENSMUSG00000021804,ENSMUSG00000031293,ENSMUSG00000037446
## 5                                        ENSMUSG00000021099,ENSMUSG00000023978,ENSMUSG00000029491,ENSMUSG00000031293,ENSMUSG00000037446,ENSMUSG00000048015,ENSMUSG00000071648
## 6                                        ENSMUSG00000021099,ENSMUSG00000023978,ENSMUSG00000029491,ENSMUSG00000031293,ENSMUSG00000037446,ENSMUSG00000048015,ENSMUSG00000071648
## 7                                        ENSMUSG00000021099,ENSMUSG00000023978,ENSMUSG00000029491,ENSMUSG00000031293,ENSMUSG00000037446,ENSMUSG00000048015,ENSMUSG00000071648
## 8  ENSMUSG00000021804,ENSMUSG00000023978,ENSMUSG00000025386,ENSMUSG00000029491,ENSMUSG00000031293,ENSMUSG00000037446,ENSMUSG00000041534,ENSMUSG00000053773,ENSMUSG00000071648
## 9  ENSMUSG00000021804,ENSMUSG00000023978,ENSMUSG00000025386,ENSMUSG00000029491,ENSMUSG00000031293,ENSMUSG00000037446,ENSMUSG00000041534,ENSMUSG00000053773,ENSMUSG00000071648
## 10                                                          ENSMUSG00000023978,ENSMUSG00000029491,ENSMUSG00000031293,ENSMUSG00000037446,ENSMUSG00000048015,ENSMUSG00000071648
## 11                                                                             ENSMUSG00000023978,ENSMUSG00000028777,ENSMUSG00000029491,ENSMUSG00000037446,ENSMUSG00000071648
## 12                                                                             ENSMUSG00000023978,ENSMUSG00000028777,ENSMUSG00000029491,ENSMUSG00000037446,ENSMUSG00000071648
## 13                                                                             ENSMUSG00000023978,ENSMUSG00000028777,ENSMUSG00000029491,ENSMUSG00000037446,ENSMUSG00000071648
## 14                                                                             ENSMUSG00000023978,ENSMUSG00000028777,ENSMUSG00000029491,ENSMUSG00000037446,ENSMUSG00000071648

down_go$pvalue_plots$bpp_plot_over

## dlgn knockout translatome vs retina knockout translatome.
dlgnko_wrt_retko_translatome <- translatome_plotter(pair_mtrx,
                                                    x_axis="koret", y_axis="kodlgn")

## No fc/lfc was provided, defaulting to 10 fold.

dlgnko_wrt_retko_translatome$plot

up_go <- simple_gprofiler(sig_genes=dlgnko_wrt_retko_translatome$ups, species="mmusculus")

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

## GO search found 2 hits.

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

## KEGG search found 3 hits.

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

## REAC search found 2 hits.

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

## MI search found 0 hits.

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

## TF search found 2 hits.

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

## CORUM search found 2 hits.

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

## HP search found 0 hits.

down_go <- simple_gprofiler(sig_genes=dlgnko_wrt_retko_translatome$downs, species="mmusculus")

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

## GO search found 4 hits.

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

## KEGG search found 9 hits.

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

## REAC search found 2 hits.

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

## MI search found 0 hits.

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

## TF search found 0 hits.

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

## CORUM search found 0 hits.

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

## HP search found 0 hits.

down_go$go

##   query.number significant p.value term.size query.size overlap.size precision
## 1            1        TRUE 0.04480       319         12            4     0.333
## 2            1        TRUE 0.01250       194         14            4     0.286
## 3            1        TRUE 0.00675       166         14            4     0.286
## 4            1        TRUE 0.00161       116         14            4     0.286
##   recall    term.id domain subgraph.number
## 1  0.013 GO:0071496     BP               2
## 2  0.021 GO:0007249     BP               1
## 3  0.024 GO:0043122     BP               1
## 4  0.034 GO:0043123     BP               1
##                                                    term.name relative.depth
## 1                     cellular response to external stimulus              1
## 2                        I-kappaB kinase/NF-kappaB signaling              1
## 3          regulation of I-kappaB kinase/NF-kappaB signaling              1
## 4 positive regulation of I-kappaB kinase/NF-kappaB signaling              2
##                                                                  intersection
## 1 ENSMUSG00000026029,ENSMUSG00000026989,ENSMUSG00000027398,ENSMUSG00000032009
## 2 ENSMUSG00000026029,ENSMUSG00000027398,ENSMUSG00000032496,ENSMUSG00000035186
## 3 ENSMUSG00000026029,ENSMUSG00000027398,ENSMUSG00000032496,ENSMUSG00000035186
## 4 ENSMUSG00000026029,ENSMUSG00000027398,ENSMUSG00000032496,ENSMUSG00000035186

## Normal dlgn translatome vs knockout dlgn translatome.
dlgnnorm_wrt_dlgnko_translatome <- translatome_plotter(pair_mtrx,
                                                       x_axis="normdlgn",
                                                       y_axis="kodlgn")

## No fc/lfc was provided, defaulting to 10 fold.

dlgnnorm_wrt_dlgnko_translatome$plot

up_go <- simple_gprofiler(sig_genes=dlgnnorm_wrt_dlgnko_translatome$ups, species="mmusculus")

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

## GO search found 0 hits.

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

## KEGG search found 0 hits.

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

## REAC search found 1 hits.

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

## MI search found 0 hits.

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

## TF search found 0 hits.

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

## CORUM search found 0 hits.

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

## HP search found 0 hits.

down_go <- simple_gprofiler(sig_genes=dlgnnorm_wrt_dlgnko_translatome$downs, species="mmusculus")

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

## GO search found 3 hits.

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

## KEGG search found 0 hits.

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

## REAC search found 0 hits.

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

## MI search found 0 hits.

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

## TF search found 0 hits.

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

## CORUM search found 1 hits.

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

## HP search found 0 hits.

## dlgn knockout translatome vs. scn knockout translatome.
dlgnko_wrt_scnko_translatome <- translatome_plotter(pair_mtrx,
                                                    x_axis="kodlgn",
                                                    y_axis="koscn")

## No fc/lfc was provided, defaulting to 10 fold.

dlgnko_wrt_scnko_translatome$plot

## Warning: Removed 1 rows containing missing values (geom_point).

up_go <- simple_gprofiler(sig_genes=dlgnko_wrt_scnko_translatome$ups, species="mmusculus")

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

## GO search found 18 hits.

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

## KEGG search found 1 hits.

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

## REAC search found 1 hits.

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

## MI search found 0 hits.

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

## TF search found 0 hits.

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

## CORUM search found 1 hits.

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

## HP search found 0 hits.

down_go <- simple_gprofiler(sig_genes=dlgnko_wrt_scnko_translatome$downs, species="mmusculus")

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

## GO search found 2 hits.

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

## KEGG search found 0 hits.

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

## REAC search found 1 hits.

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

## MI search found 0 hits.

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

## TF search found 0 hits.

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

## CORUM search found 3 hits.

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

## HP search found 0 hits.

up_go$pvalue_plots$bpp_plot_over

4 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] 1877

ror_down <- ror[down_idx]
length(ror_down)

## [1] 476

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 1877 genes against mmusculus.

## GO search found 271 hits.

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

## KEGG search found 8 hits.

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

## REAC search found 11 hits.

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

## MI search found 0 hits.

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

## TF search found 376 hits.

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

## CORUM search found 1 hits.

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

## HP search found 5 hits.

## Writing data to: excel/20200115mm_ror_gpfoiler_up-v20200114.xlsx.

## Error in `levels<-`(`*tmp*`, value = as.character(levels)) : 
##   factor level [28] is duplicated

## Error in `levels<-`(`*tmp*`, value = as.character(levels)) : 
##   factor level [28] is duplicated

## Error in `levels<-`(`*tmp*`, value = as.character(levels)) : 
##   factor level [28] is duplicated

## Error in `levels<-`(`*tmp*`, value = as.character(levels)) : 
##   factor level [28] is duplicated

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

## Error in `levels<-`(`*tmp*`, value = as.character(levels)): factor level [28] is duplicated

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 476 genes against mmusculus.

## GO search found 82 hits.

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

## KEGG search found 7 hits.

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

## REAC search found 2 hits.

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

## MI search found 0 hits.

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

## TF search found 6 hits.

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

## CORUM search found 0 hits.

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

## HP search found 0 hits.

## Writing data to: excel/20200115mm_ror_gpfoiler_down-v20200114.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)