Some Pseudomonas RNAseq data: SK.

1 Introduction

This dataset contains multiple experiments.

2 Annotations

pa14_gff <- load_gff_annotations("reference/paeruginosa_pa14.gff", id_col="gene_id")

## Trying attempt: rtracklayer::import.gff3(gff, sequenceRegionsAsSeqinfo = TRUE)

## Had a successful gff import with rtracklayer::import.gff3(gff, sequenceRegionsAsSeqinfo = TRUE)

## Returning a df with 16 columns and 11946 rows.

rownames(pa14_gff) <- pa14_gff[["gene_id"]]
## The Alias column has PA14_00010

pa14_microbes <- as.data.frame(load_microbesonline_annotations("PA14"))

## Found 1 entry.

## Pseudomonas aeruginosa UCBPP-PA14Proteobacteria2006-11-22yes105972208963

## The species being downloaded is: Pseudomonas aeruginosa UCBPP-PA14

## Downloading: http://www.microbesonline.org/cgi-bin/genomeInfo.cgi?tId=208963;export=tab

## The sysName column has PA14_0010

pa14_annot <- merge(pa14_gff, pa14_microbes, by.x="Alias", by.y="sysName")
rownames(pa14_annot) <- pa14_annot[["gene_id"]]

pa14_go <- load_microbesonline_go(species="PA14")

## Found 1 entry.

## Pseudomonas aeruginosa UCBPP-PA14Proteobacteria2006-11-22yes105972208963

## The species being downloaded is: Pseudomonas aeruginosa UCBPP-PA14 and is being downloaded as 208963.tab.

3 Make the expressionset

Given the above annotations, now lets pull in the counts.

I am switching to the sheet all_samples_modified_gcd.xlsx for the moment because I added a space in the strain name for the gcd sampl

pa14_expt <- create_expt("sample_sheets/all_samples_modified_gcd.xlsx",
                         gene_info=pa14_annot, file_column="hisatcounttable")

## Reading the sample metadata.

## Did not find the condition column in the sample sheet.

## Filling it in as undefined.

## Did not find the batch column in the sample sheet.

## Filling it in as undefined.

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

## Matched 5972 annotations and counts.

## Bringing together the count matrix and gene information.

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

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

## The final expressionset has 5979 rows and 105 columns.

4 Quick global look at the data

While we are at it, lets drop the two sad samples.

pa14_libsize <- plot_libsize(pa14_expt)
pa14_libsize$plot

pa14_nonzero <- plot_nonzero(pa14_expt)
pa14_nonzero$plot

## Warning: ggrepel: 100 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

pa14_expt <- subset_expt(pa14_expt, nonzero=5500)

## The samples (and read coverage) removed when filtering 5500 non-zero genes are:

##  SM040  SM048 
##  43984 316632

## subset_expt(): There were 105, now there are 103 samples.

5 View all samples briefly

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

## Error in normalize_expt(sk_expt, transform = "log2", convert = "cpm", : object 'sk_expt' not found

plot_pca(sk_norm)$plot

## Error in plot_pca(sk_norm): object 'sk_norm' not found

6 Alginate!

initials_factor <- gsub(x=rownames(pData(pa14_expt)), pattern="^(..).*$", replacement="\\1")
pData(pa14_expt)[["initials"]] <- as.factor(initials_factor)
sk_expt <- subset_expt(pa14_expt, subset="initials=='SK'")

## subset_expt(): There were 103, now there are 45 samples.

alginate_expt <- subset_expt(sk_expt, subset="media!='LB'") %>%
  set_expt_batches(fact="bioreplicate") %>%
  set_expt_conditions(fact="media")

## subset_expt(): There were 45, now there are 36 samples.

rph_expt <- subset_expt(sk_expt, subset="media=='LB'") %>%
  set_expt_batches(fact="bioreplicate") %>%
  set_expt_conditions(fact="strains") %>%
  sanitize_expt_metadata()

## subset_expt(): There were 45, now there are 9 samples.

7 The alginate data

There are two media, an inducer, and multiple methods of reducing expression…

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

## Removing 18 low-count genes (5961 remaining).

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

plot_pca(alg_norm)$plot

## Warning: ggrepel: 1 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

conditions <- pData(alginate_expt)[["condition"]]
induce_bool <- grepl(pattern="IPTG", x=conditions)
induce_factor <- rep("not_induced", length(induce_bool))
induce_factor[induce_bool] <- "IPTG"

reduce_bool <- grepl(pattern="Eb", x=conditions)
reduce_factor <- rep("not_reduced", length(reduce_bool))
eb_samples <-  grepl(pattern="Eb", x=conditions)
reduce_factor[eb_samples] <- "Eb"
ebo_samples <- grepl(pattern="EbO", x=conditions)
reduce_factor[ebo_samples] <- "EbO"
ebs_samples <- grepl(pattern="EbS", x=conditions)
reduce_factor[ebs_samples] <- "EbS"

media_factor <- conditions
lb_bool <- grepl(pattern="LB", x=conditions)
media_factor[lb_bool] <- "LB"
media_factor[!lb_bool] <- "SCFM"

pData(alginate_expt)[["media_add"]] <- pData(alginate_expt)[["media"]]
pData(alginate_expt)[["media"]] <- media_factor
pData(alginate_expt)[["induced"]] <- as.factor(induce_factor)
pData(alginate_expt)[["reduced"]] <- as.factor(reduce_factor)

7.1 Compare media

alginate_media <- set_expt_conditions(alginate_expt, fact="media") %>%
  set_expt_batches(fact="bioreplicate")

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

## Removing 18 low-count genes (5961 remaining).

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

plot_pca(alginate_media_norm)$plot

## Warning: ggrepel: 6 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

7.2 Compare induce

alginate_induce <- set_expt_conditions(alginate_expt, fact="induced") %>%
  set_expt_batches(fact="bioreplicate")
alginate_induce_norm <- normalize_expt(alginate_induce, transform="log2",
                                       convert="cpm", norm="quant", filter=TRUE)

## Removing 18 low-count genes (5961 remaining).

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

plot_pca(alginate_induce_norm)$plot

## Warning: ggrepel: 4 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

alginate_reduce <- set_expt_conditions(alginate_expt, fact="reduced") %>%
  set_expt_batches(fact="bioreplicate")
alginate_reduce_norm <- normalize_expt(alginate_reduce, transform="log2",
                                       convert="cpm", norm="quant", filter=TRUE)

## Removing 18 low-count genes (5961 remaining).

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

plot_pca(alginate_reduce_norm)$plot

## Warning: ggrepel: 11 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

both_fact <- paste0(pData(alginate_expt)[["induced"]], "_",
                    pData(alginate_expt)[["reduced"]])
pData(alginate_expt)[["induce_reduce"]] <- both_fact

alginate_both <- set_expt_conditions(alginate_expt, fact="induce_reduce") %>%
  set_expt_batches(fact="bioreplicate")
alginate_both_norm <- normalize_expt(alginate_both, transform="log2", batch="svaseq",
                                     convert="cpm", norm="quant", filter=TRUE)

## Warning in normalize_expt(alginate_both, transform = "log2", batch = "svaseq", :
## Quantile normalization and sva do not always play well together.

## Removing 18 low-count genes (5961 remaining).

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

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

## Setting 60 low elements to zero.

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

plot_pca(alginate_both_norm)$plot

8 rph data

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

## Removing 71 low-count genes (5908 remaining).

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

plot_pca(rph_norm)$plot

rph_de <- all_pairwise(rph_expt, model_batch=TRUE, filter=TRUE)

## Plotting a PCA before surrogate/batch inclusion.

## Using limma's removeBatchEffect to visualize with(out) batch inclusion.

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

rph_tables <- combine_de_tables(
    rph_de,
    excel=glue::glue("excel/rph_tables-v{ver}.xlsx"))

## Deleting the file excel/rph_tables-v20220131.xlsx before writing the tables.

rph_sig <- extract_significant_genes(
    rph_tables,
    excel=glue::glue("excel/rph_sig-v{ver}.xlsx"))

## Deleting the file excel/rph_sig-v20220131.xlsx before writing the tables.

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