Compare polydT to ribozero library preparations.

1 Introduction

I aim to demonstrate the similarities/differences between a series of libraries which were prepared once using polydT beads and again using the ribozero kit. We have on hand 6 samples for which both methods were employed.

I therefore took a new download of our persistence shared sample sheet and culled it to just these 6 pairs: ‘compare_mRNA_ribozero.xlsx’.

sample_sheet <- "sample_sheets/compare_mRNA_ribozero.xlsx"
modified_sheet <- "sample_sheets/compare_mRNA_ribozero_modified.xlsx"

I should only need to run gather_metadata once in order to generate a new xlsx file which I will read in the future. There is one caveat, I need to either delete the first column of the original sample sheet or fill in some arguments. Because I am lazy and cannot remember the argument for setting the rownames… I deleted the first column.

modified <- gather_preprocessing_metadata(sample_sheet)

2 Grab the annotation data

hs_annot <- load_biomart_annotations(year = "2023")

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

specification <- make_rnaseq_spec()
basedir <- "preprocessing"
species <- "hg38_111"
modified <- gather_preprocessing_metadata(sample_sheet, specification = specification,
                                          species = species, feature_type = "gene")
head(modified[["new_meta"]])
ncRNA_gene_spec <- list(
  "hisat_observed_genes" = list(
    "file" = "{basedir}/{meta[['sampleid']]}/outputs/*hisat*_{species}/{species}_*ncRNA_gene*.count.xz",
    "column" = "ncRNA_genes_observed"),
  "hisat_sum_genes" = list(
    "file" = "{basedir}/{meta[['sampleid']]}/outputs/*hisat*_{species}/{species}_*ncRNA_gene*.count.xz",
    "column" = "ncRNA_gene_reads"),
  "hisat_count_table" = list(
    "file" = "{basedir}/{meta[['sampleid']]}/outputs/*hisat*_{species}/{species}_*ncRNA_gene*.count.xz",
    "column" = "ncRNA_gene_count_table"))
appended <- gather_preprocessing_metadata(modified[["new_meta"]], species = species,
                                          type = "genome", specification = ncRNA_gene_spec,
                                          feature_type = "ncRNA_gene",
                                          new_metadata = "sample_sheets/minimal_with_gene_nc.xlsx")
head(appended[["new_meta"]])
sno_spec <- list(
  "hisat_observed_genes" = list(
    "file" = "{basedir}/{meta[['sampleid']]}/outputs/*hisat*_{species}/{species}_*snoRNA*.count.xz",
    "column" = "snoRNA_observed"),
  "hisat_sum_genes" = list(
    "file" = "{basedir}/{meta[['sampleid']]}/outputs/*hisat*_{species}/{species}_*snoRNA*.count.xz",
    "column" = "snoRNA_reads"),
  "hisat_count_table" = list(
    "file" = "{basedir}/{meta[['sampleid']]}/outputs/*hisat*_{species}/{species}_*snoRNA*.count.xz",
    "column" = "snoRNA_count_table"))
appended <- gather_preprocessing_metadata(appended[["new_meta"]], species = species,
                                          type = "genome", specification = sno_spec,
                                          feature_type = "snoRNA",
                                          new_metadata = "sample_sheets/minimal_with_gene_nc_sno.xlsx")
head(appended[["new_meta"]])
exon_spec <- list(
  "hisat_observed_genes" = list(
    "file" = "{basedir}/{meta[['sampleid']]}/outputs/*hisat*_{species}/{species}_*exon*.count.xz",
    "column" = "exons_observed"),
  "hisat_sum_genes" = list(
    "file" = "{basedir}/{meta[['sampleid']]}/outputs/*hisat*_{species}/{species}_*exon*.count.xz",
    "column" = "exon_reads"),
  "hisat_count_table" = list(
    "file" = "{basedir}/{meta[['sampleid']]}/outputs/*hisat*_{species}/{species}_*exon*.count.xz",
    "column" = "exon_count_table"))
appended <- gather_preprocessing_metadata(appended[["new_meta"]], species = species,
                                          type = "genome", specification = exon_spec,
                                          feature_type = "exon",
                                          new_metadata = "sample_sheets/minimal_with_gene_nc_sno_exon.xlsx")
head(appended[["new_meta"]])
pseudogenic_transcript_spec <- list(
  "hisat_observed_genes" = list(
    "file" = "{basedir}/{meta[['sampleid']]}/outputs/*hisat*_{species}/{species}_*pseudogenic_transcript*.count.xz",
    "column" = "pseudogenic_transcripts_observed"),
  "hisat_sum_genes" = list(
    "file" = "{basedir}/{meta[['sampleid']]}/outputs/*hisat*_{species}/{species}_*pseudogenic_transcript*.count.xz",
    "column" = "pseudogenic_transcript_reads"),
  "hisat_count_table" = list(
    "file" = "{basedir}/{meta[['sampleid']]}/outputs/*hisat*_{species}/{species}_*pseudogenic_transcript*.count.xz",
    "column" = "pseudogenic_transcript_count_table"))
appended <- gather_preprocessing_metadata(appended[["new_meta"]], species = species,
                                          type = "genome", specification = pseudogenic_transcript_spec,
                                          feature_type = "pseudogenic_transcript",
                                          new_metadata = "sample_sheets/minimal_with_gene_nc_sno_exon_pseudo.xlsx")
head(appended[["new_meta"]])
lnc_RNA_spec <- list(
  "hisat_observed_genes" = list(
    "file" = "{basedir}/{meta[['sampleid']]}/outputs/*hisat*_{species}/{species}_*lnc_RNA*.count.xz",
    "column" = "lnc_RNAs_observed"),
  "hisat_sum_genes" = list(
    "file" = "{basedir}/{meta[['sampleid']]}/outputs/*hisat*_{species}/{species}_*lnc_RNA*.count.xz",
    "column" = "lnc_RNA_reads"),
  "hisat_count_table" = list(
    "file" = "{basedir}/{meta[['sampleid']]}/outputs/*hisat*_{species}/{species}_*lnc_RNA*.count.xz",
    "column" = "lnc_RNA_count_table"))
appended <- gather_preprocessing_metadata(appended[["new_meta"]], species = species,
                                          type = "genome", specification = lnc_RNA_spec,
                                          feature_type = "lnc_RNA", verbose = TRUE,
                                          new_metadata = "sample_sheets/minimal_with_gene_nc_sno_exon_pseudo_lnc.xlsx")
head(appended[["new_meta"]])
five_prime_UTR_spec <- list(
  "hisat_observed_genes" = list(
    "file" = "{basedir}/{meta[['sampleid']]}/outputs/*hisat*_{species}/{species}_*five_prime_UTR*.count.xz",
    "column" = "five_prime_UTRs_observed"),
  "hisat_sum_genes" = list(
    "file" = "{basedir}/{meta[['sampleid']]}/outputs/*hisat*_{species}/{species}_*five_prime_UTR*.count.xz",
    "column" = "five_prime_UTR_reads"),
  "hisat_count_table" = list(
    "file" = "{basedir}/{meta[['sampleid']]}/outputs/*hisat*_{species}/{species}_*five_prime_UTR*.count.xz",
    "column" = "five_prime_UTR_count_table"))
appended <- gather_preprocessing_metadata(appended[["new_meta"]], species = species,
                                          type = "genome", specification = five_prime_UTR_spec,
                                          feature_type = "five_prime_UTR",
                                          new_metadata = "sample_sheets/minimal_with_gene_nc_sno_exon_pseudo_lnc_fivep.xlsx")
head(appended[["new_meta"]])
three_prime_UTR_spec <- list(
  "hisat_observed_genes" = list(
    "file" = "{basedir}/{meta[['sampleid']]}/outputs/*hisat*_{species}/{species}_*three_prime_UTR*.count.xz",
    "column" = "three_prime_UTRs_observed"),
  "hisat_sum_genes" = list(
    "file" = "{basedir}/{meta[['sampleid']]}/outputs/*hisat*_{species}/{species}_*three_prime_UTR*.count.xz",
    "column" = "three_prime_UTR_reads"),
  "hisat_count_table" = list(
    "file" = "{basedir}/{meta[['sampleid']]}/outputs/*hisat*_{species}/{species}_*three_prime_UTR*.count.xz",
    "column" = "three_prime_UTR_count_table"))
appended <- gather_preprocessing_metadata(appended[["new_meta"]], species = species,
                                          type = "genome", specification = three_prime_UTR_spec,
                                          feature_type = "three_prime_UTR",
                                          new_metadata = "sample_sheets/minimal_with_gene_nc_sno_exon_pseudo_lnc_fivep_threep.xlsx")
head(appended[["new_meta"]])
pre_rRNA_spec <- list(
  "hisat_rrna_single_concordant" = list(
    "file" = "{basedir}/{meta[['sampleid']]}/outputs/*hisat*_{species}/hisat*_*rRNA*.stderr",
    "column" = "rRNA_reads_pre"),
  "hisat_rrna_percent_log" = list(
    "file" = "{basedir}/{meta[['sampleid']]}/outputs/*hisat*_{species}/hisat*_*rRNA*.stderr",
    "column" = "rRNA_pct_log_pre"))
appended <- gather_preprocessing_metadata(appended[["new_meta"]], species = "hglp",
                                          type = "rRNA", specification = pre_rRNA_spec,
                                          feature_type = "rRNA", verbose = verbose,
                                          new_metadata = "sample_sheets/minimal_with_gene_nc_sno_exon_pseudo_lnc_fivep_threep_rRNApre.xlsx")
head(appended[["new_meta"]])
post_rRNA_spec <- list(
  "hisat_observed_genes" = list(
    "file" = "{basedir}/{meta[['sampleid']]}/outputs/*hisat*_{species}/{species}_*rRNA*.count.xz",
    "column" = "rRNAs_observed"),
  "hisat_sum_genes" = list(
    "file" = "{basedir}/{meta[['sampleid']]}/outputs/*hisat*_{species}/{species}_*rRNA*.count.xz",
    "column" = "rRNA_reads"),
  "hisat_count_table" = list(
    "file" = "{basedir}/{meta[['sampleid']]}/outputs/*hisat*_{species}/{species}_*rRNA*.count.xz",
    "column" = "rRNA_count_table"))
appended <- gather_preprocessing_metadata(appended[["new_meta"]], species = species,
                                          type = "genome", specification = post_rRNA_spec,
                                          feature_type = "rRNA", verbose = TRUE,
                                          new_metadata = "sample_sheets/minimal_with_gene_nc_sno_exon_pseudo_lnc_fivep_threep_rRNApre_rRNApost.xlsx")
head(appended[["new_meta"]])

3 Make a summarizedExperiment

hs_se <- create_se(modified_sheet, gene_info = hs_annot,
                   file_column = "hisat_count_table") %>%
  set_se_conditions(fact = "library_type")

## Reading the sample metadata.

## Checking the state of the condition column.

## Checking the state of the batch column.

## Checking the condition factor.

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

## Defauting to the gene annotations.

## Matched 21557 annotations and counts.

## Bringing together the count matrix and gene information.

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

## The final summarized experiment has 21557 rows and 72 columns.

## The numbers of samples by condition are:

## 
## mRNA   RZ 
##    6    6

plot_libsize(hs_se)

## Library sizes of 12 samples, 
## ranging from 10,852,469 to 16,751,053.

plot_legend(hs_se)

## Warning in RColorBrewer::brewer.pal(12, "Dark2"): n too large, allowed maximum for palette Dark2 is 8
## Returning the palette you asked for with that many colors

## The colors used in the expressionset are: #7570B3, #1B9E77.

hs_se_ncrna <- create_se(appended[["new_meta"]], gene_info = hs_annot,
                         file_column = "ncRNA_gene_count_table")

## Error: object 'appended' not found

hs_se_ncrna <- set_se_conditions(hs_se_ncrna, fact = "library_type")

## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'colData': object 'hs_se_ncrna' not found

plot_libsize(hs_se_ncrna)

## Error in h(simpleError(msg, call)): error in evaluating the argument 'data' in selecting a method for function 'plot_libsize': object 'hs_se_ncrna' not found

hs_se_snorna <- create_se(appended[["new_meta"]], gene_info = hs_annot,
                         file_column = "snoRNA_count_table")

## Error: object 'appended' not found

hs_se_snorna <- set_se_conditions(hs_se_snorna, fact = "library_type")

## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'colData': object 'hs_se_snorna' not found

plot_libsize(hs_se_snorna)

## Error in h(simpleError(msg, call)): error in evaluating the argument 'data' in selecting a method for function 'plot_libsize': object 'hs_se_snorna' not found

hs_se_exons <- create_se(appended[["new_meta"]],
                         file_column = "exon_count_table")

## Error: object 'appended' not found

hs_se_exons <- set_se_conditions(hs_se_exons, fact = "library_type")

## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'colData': object 'hs_se_exons' not found

plot_libsize(hs_se_exons)

## Error in h(simpleError(msg, call)): error in evaluating the argument 'data' in selecting a method for function 'plot_libsize': object 'hs_se_exons' not found

hs_lnc <- create_se(appended[["new_meta"]],
                    file_column = "lnc_RNA_count_table")

## Error: object 'appended' not found

hs_lnc <- set_se_conditions(hs_lnc, fact = "library_type")

## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'colData': object 'hs_lnc' not found

plot_libsize(hs_lnc)

## Error in h(simpleError(msg, call)): error in evaluating the argument 'data' in selecting a method for function 'plot_libsize': object 'hs_lnc' not found

hs_pseudo <- create_se(appended[["new_meta"]],
                       file_column = "pseudogenic_transcript_count_table")

## Error: object 'appended' not found

hs_pseudo <- set_se_conditions(hs_pseudo, fact = "library_type")

## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'colData': object 'hs_pseudo' not found

plot_libsize(hs_pseudo)

## Error in h(simpleError(msg, call)): error in evaluating the argument 'data' in selecting a method for function 'plot_libsize': object 'hs_pseudo' not found

hs_fivep <- create_se(appended[["new_meta"]],
                       file_column = "five_prime_UTR_count_table")

## Error: object 'appended' not found

hs_fivep <- set_se_conditions(hs_fivep, fact = "library_type")

## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'colData': object 'hs_fivep' not found

plot_libsize(hs_fivep)

## Error in h(simpleError(msg, call)): error in evaluating the argument 'data' in selecting a method for function 'plot_libsize': object 'hs_fivep' not found

hs_threep <- create_se(appended[["new_meta"]],
                       file_column = "three_prime_UTR_count_table")

## Error: object 'appended' not found

hs_threep <- set_se_conditions(hs_threep, fact = "library_type")

## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'colData': object 'hs_threep' not found

plot_libsize(hs_threep)

## Error in h(simpleError(msg, call)): error in evaluating the argument 'data' in selecting a method for function 'plot_libsize': object 'hs_threep' not found

4 Play with the metadata

meta <- pData(hs_se)
meta[["nc_vs_gene"]] <- meta[["ncRNA_genes_observed"]] / meta[["hisat_observed_genes"]]

## Error in `[[<-`(`*tmp*`, "nc_vs_gene", value = numeric(0)): 0 elements in value to replace 12 elements

meta[["lnc_vs_gene"]] <- meta[["lnc_RNAs_observed"]] / meta[["hisat_observed_genes"]]

## Error in `[[<-`(`*tmp*`, "lnc_vs_gene", value = numeric(0)): 0 elements in value to replace 12 elements

meta[["sno_vs_gene"]] <- meta[["snoRNA_observed"]] / meta[["hisat_observed_genes"]]

## Error in `[[<-`(`*tmp*`, "sno_vs_gene", value = numeric(0)): 0 elements in value to replace 12 elements

meta[["exon_vs_gene"]] <- meta[["exons_observed"]] / meta[["hisat_observed_genes"]]

## Error in `[[<-`(`*tmp*`, "exon_vs_gene", value = numeric(0)): 0 elements in value to replace 12 elements

meta[["pseudo_vs_gene"]] <- meta[["pseudogenic_transcripts_observed"]] / meta[["hisat_observed_genes"]]

## Error in `[[<-`(`*tmp*`, "pseudo_vs_gene", value = numeric(0)): 0 elements in value to replace 12 elements

meta[["fivep_vs_gene"]] <- meta[["five_prime_UTRs_observed"]] / meta[["hisat_observed_genes"]]

## Error in `[[<-`(`*tmp*`, "fivep_vs_gene", value = numeric(0)): 0 elements in value to replace 12 elements

meta[["threep_vs_gene"]] <- meta[["three_prime_UTRs_observed"]] / meta[["hisat_observed_genes"]]

## Error in `[[<-`(`*tmp*`, "threep_vs_gene", value = numeric(0)): 0 elements in value to replace 12 elements

meta[["nc_vs_gene_reads"]] <- meta[["ncRNA_gene_reads"]] / meta[["hisat_sum_genes"]]

## Error in `[[<-`(`*tmp*`, "nc_vs_gene_reads", value = numeric(0)): 0 elements in value to replace 12 elements

meta[["lnc_vs_gene_reads"]] <- meta[["lnc_RNA_reads"]] / meta[["hisat_sum_genes"]]

## Error in `[[<-`(`*tmp*`, "lnc_vs_gene_reads", value = numeric(0)): 0 elements in value to replace 12 elements

meta[["sno_vs_gene_reads"]] <- meta[["snoRNA_reads"]] / meta[["hisat_sum_genes"]]

## Error in `[[<-`(`*tmp*`, "sno_vs_gene_reads", value = numeric(0)): 0 elements in value to replace 12 elements

meta[["exon_vs_gene_reads"]] <- meta[["exon_reads"]] / meta[["hisat_sum_genes"]]

## Error in `[[<-`(`*tmp*`, "exon_vs_gene_reads", value = numeric(0)): 0 elements in value to replace 12 elements

meta[["pseudo_vs_gene_reads"]] <- meta[["pseudogenic_transcript_reads"]] / meta[["hisat_sum_genes"]]

## Error in `[[<-`(`*tmp*`, "pseudo_vs_gene_reads", value = numeric(0)): 0 elements in value to replace 12 elements

meta[["fivep_vs_gene_reads"]] <- meta[["five_prime_UTR_reads"]] / meta[["hisat_sum_genes"]]

## Error in `[[<-`(`*tmp*`, "fivep_vs_gene_reads", value = numeric(0)): 0 elements in value to replace 12 elements

meta[["threep_vs_gene_reads"]] <- meta[["three_prime_UTR_reads"]] / meta[["hisat_sum_genes"]]

## Error in `[[<-`(`*tmp*`, "threep_vs_gene_reads", value = numeric(0)): 0 elements in value to replace 12 elements

pData(hs_se) <- meta
hs_se <- set_se_conditions(hs_se, fact = "library_type")

## The numbers of samples by condition are:

## 
## mRNA   RZ 
##    6    6

genes_obs <- plot_metadata_factors(hs_se, column = "hisat_observed_genes")
genes_obs

nc_vs_gene <- plot_metadata_factors(hs_se, column = "nc_vs_gene")
nc_vs_gene

## Error in `ggplot2::geom_violin()` at hpgltools/R/metadata.R:2018:5:
## ! Problem while computing aesthetics.
## ℹ Error occurred in the 1st layer.
## Caused by error in `.data[["nc_vs_gene"]]`:
## ! Column `nc_vs_gene` not found in `.data`.

lnc_vs_gene <- plot_metadata_factors(hs_se, column = "lnc_vs_gene")
lnc_vs_gene

## Error in `ggplot2::geom_violin()` at hpgltools/R/metadata.R:2018:5:
## ! Problem while computing aesthetics.
## ℹ Error occurred in the 1st layer.
## Caused by error in `.data[["lnc_vs_gene"]]`:
## ! Column `lnc_vs_gene` not found in `.data`.

sno_vs_gene <- plot_metadata_factors(hs_se, column = "sno_vs_gene")
sno_vs_gene

## Error in `ggplot2::geom_violin()` at hpgltools/R/metadata.R:2018:5:
## ! Problem while computing aesthetics.
## ℹ Error occurred in the 1st layer.
## Caused by error in `.data[["sno_vs_gene"]]`:
## ! Column `sno_vs_gene` not found in `.data`.

exon_vs_gene <- plot_metadata_factors(hs_se, column = "exon_vs_gene")
exon_vs_gene

## Error in `ggplot2::geom_violin()` at hpgltools/R/metadata.R:2018:5:
## ! Problem while computing aesthetics.
## ℹ Error occurred in the 1st layer.
## Caused by error in `.data[["exon_vs_gene"]]`:
## ! Column `exon_vs_gene` not found in `.data`.

fivep_vs_gene <- plot_metadata_factors(hs_se, column = "fivep_vs_gene")
fivep_vs_gene

## Error in `ggplot2::geom_violin()` at hpgltools/R/metadata.R:2018:5:
## ! Problem while computing aesthetics.
## ℹ Error occurred in the 1st layer.
## Caused by error in `.data[["fivep_vs_gene"]]`:
## ! Column `fivep_vs_gene` not found in `.data`.

threep_vs_gene <- plot_metadata_factors(hs_se, column = "threep_vs_gene")
threep_vs_gene

## Error in `ggplot2::geom_violin()` at hpgltools/R/metadata.R:2018:5:
## ! Problem while computing aesthetics.
## ℹ Error occurred in the 1st layer.
## Caused by error in `.data[["threep_vs_gene"]]`:
## ! Column `threep_vs_gene` not found in `.data`.

reads_obs <- plot_metadata_factors(hs_se, column = "hisat_sum_genes")
reads_obs

nc_vs_gene_reads <- plot_metadata_factors(hs_se, column = "nc_vs_gene_reads")
nc_vs_gene_reads

## Error in `ggplot2::geom_violin()` at hpgltools/R/metadata.R:2018:5:
## ! Problem while computing aesthetics.
## ℹ Error occurred in the 1st layer.
## Caused by error in `.data[["nc_vs_gene_reads"]]`:
## ! Column `nc_vs_gene_reads` not found in `.data`.

lnc_vs_gene_reads <- plot_metadata_factors(hs_se, column = "lnc_vs_gene_reads")
lnc_vs_gene_reads

## Error in `ggplot2::geom_violin()` at hpgltools/R/metadata.R:2018:5:
## ! Problem while computing aesthetics.
## ℹ Error occurred in the 1st layer.
## Caused by error in `.data[["lnc_vs_gene_reads"]]`:
## ! Column `lnc_vs_gene_reads` not found in `.data`.

sno_vs_gene_reads <- plot_metadata_factors(hs_se, column = "sno_vs_gene_reads")
sno_vs_gene_reads

## Error in `ggplot2::geom_violin()` at hpgltools/R/metadata.R:2018:5:
## ! Problem while computing aesthetics.
## ℹ Error occurred in the 1st layer.
## Caused by error in `.data[["sno_vs_gene_reads"]]`:
## ! Column `sno_vs_gene_reads` not found in `.data`.

exon_vs_gene_reads <- plot_metadata_factors(hs_se, column = "exon_vs_gene_reads")
exon_vs_gene_reads

## Error in `ggplot2::geom_violin()` at hpgltools/R/metadata.R:2018:5:
## ! Problem while computing aesthetics.
## ℹ Error occurred in the 1st layer.
## Caused by error in `.data[["exon_vs_gene_reads"]]`:
## ! Column `exon_vs_gene_reads` not found in `.data`.

fivep_vs_gene_reads <- plot_metadata_factors(hs_se, column = "fivep_vs_gene_reads")
fivep_vs_gene_reads

## Error in `ggplot2::geom_violin()` at hpgltools/R/metadata.R:2018:5:
## ! Problem while computing aesthetics.
## ℹ Error occurred in the 1st layer.
## Caused by error in `.data[["fivep_vs_gene_reads"]]`:
## ! Column `fivep_vs_gene_reads` not found in `.data`.

threep_vs_gene_reads <- plot_metadata_factors(hs_se, column = "threep_vs_gene_reads")
threep_vs_gene_reads

## Error in `ggplot2::geom_violin()` at hpgltools/R/metadata.R:2018:5:
## ! Problem while computing aesthetics.
## ℹ Error occurred in the 1st layer.
## Caused by error in `.data[["threep_vs_gene_reads"]]`:
## ! Column `threep_vs_gene_reads` not found in `.data`.

5 Examine kraken bacterial data

hs_kraken <- create_expt(modified_sheet, file_type = "kraken", file_column = "kraken_matrix") %>%
  set_expt_conditions(fact = "library_type")

## Reading the sample metadata.

## Checking the state of the condition column.

## Checking the state of the batch column.

## Checking the condition factor.

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

## Warning in first_rownames != current_rownames: longer object length is not a multiple of shorter object
## length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, : The file:
## /lab/home/trey/sshfs/scratch/atb/rnaseq/lpanamensis_persistence_2023/preprocessing/PRHU0038/outputs/06kraken_bacteria/kraken_report_matrix.tsv
## has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a multiple of shorter object
## length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, : The file:
## /lab/home/trey/sshfs/scratch/atb/rnaseq/lpanamensis_persistence_2023/preprocessing/PRHU0004/outputs/06kraken_bacteria/kraken_report_matrix.tsv
## has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a multiple of shorter object
## length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, : The file:
## /lab/home/trey/sshfs/scratch/atb/rnaseq/lpanamensis_persistence_2023/preprocessing/PRHU0031/outputs/06kraken_bacteria/kraken_report_matrix.tsv
## has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a multiple of shorter object
## length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, : The file:
## /lab/home/trey/sshfs/scratch/atb/rnaseq/lpanamensis_persistence_2023/preprocessing/PRHU0008/outputs/06kraken_bacteria/kraken_report_matrix.tsv
## has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a multiple of shorter object
## length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, : The file:
## /lab/home/trey/sshfs/scratch/atb/rnaseq/lpanamensis_persistence_2023/preprocessing/PRHU0037/outputs/06kraken_bacteria/kraken_report_matrix.tsv
## has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a multiple of shorter object
## length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, : The file:
## /lab/home/trey/sshfs/scratch/atb/rnaseq/lpanamensis_persistence_2023/preprocessing/PRHU0006/outputs/06kraken_bacteria/kraken_report_matrix.tsv
## has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a multiple of shorter object
## length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, : The file:
## /lab/home/trey/sshfs/scratch/atb/rnaseq/lpanamensis_persistence_2023/preprocessing/PRHU0035/outputs/06kraken_bacteria/kraken_report_matrix.tsv
## has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a multiple of shorter object
## length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, : The file:
## /lab/home/trey/sshfs/scratch/atb/rnaseq/lpanamensis_persistence_2023/preprocessing/PRHU0007/outputs/06kraken_bacteria/kraken_report_matrix.tsv
## has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a multiple of shorter object
## length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, : The file:
## /lab/home/trey/sshfs/scratch/atb/rnaseq/lpanamensis_persistence_2023/preprocessing/PRHU0036/outputs/06kraken_bacteria/kraken_report_matrix.tsv
## has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a multiple of shorter object
## length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, : The file:
## /lab/home/trey/sshfs/scratch/atb/rnaseq/lpanamensis_persistence_2023/preprocessing/PRHU0005/outputs/06kraken_bacteria/kraken_report_matrix.tsv
## has mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a multiple of shorter object
## length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, : The file:
## /lab/home/trey/sshfs/scratch/atb/rnaseq/lpanamensis_persistence_2023/preprocessing/PRHU0034/outputs/06kraken_bacteria/kraken_report_matrix.tsv
## has mismatched rownames.

## Warning in create_expt(modified_sheet, file_type = "kraken", file_column = "kraken_matrix"): There are
## some NAs in this data, the 'handle_nas' parameter may be required.

## Matched 175 annotations and counts.

## Bringing together the count matrix and gene information.

## Saving the expressionset to 'expt.rda'.

## The final expressionset has 175 features and 12 samples.

## The numbers of samples by condition are:

## 
## mRNA   RZ 
##    6    6

plot_libsize(hs_kraken)

## Library sizes of 12 samples, 
## ranging from 56,929 to 594,340.

hs_kraken_norm <- normalize_expt(hs_kraken, filter = TRUE, convert = "cpm", norm = "quant")

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

plot_pca(hs_kraken_norm)

## The result of performing a fast_svd dimension reduction.
## The x-axis is PC1 and the y-axis is PC2
## Colors are defined by mRNA, RZ
## Shapes are defined by undefined.

kraken_de <- all_pairwise(hs_kraken, filter = TRUE)

## 
## mRNA   RZ 
##    6    6 
## 
## undefined 
##        12

## Basic step 0/3: Normalizing data.

## Basic step 0/3: Converting data.

## Basic step 0/3: Transforming data.

## converting counts to integer mode

## gene-wise dispersion estimates

## mean-dispersion relationship

## -- note: fitType='parametric', but the dispersion trend was not well captured by the
##    function: y = a/x + b, and a local regression fit was automatically substituted.
##    specify fitType='local' or 'mean' to avoid this message next time.

## final dispersion estimates

## libsize was not specified, this parameter has profound effects on limma's result.

## Using the libsize from expt$best_libsize.

## Warning in createContrastL(objFlt$formula, objFlt$data, L): Contrasts with only a single non-zero term are
## already evaluated by default.

## libsize was not specified, this parameter has profound effects on limma's result.

## Using the libsize from expt$best_libsize.

kraken_table <- combine_de_tables(kraken_de, excel = "excel/kraken_compare.xlsx")

## Deleting the file excel/kraken_compare.xlsx before writing the tables.