DESeq2 plot functions not found with bcbioRNASeq

[VictoriaPatten]

Hi, I am using bcbioRNASeq and working from the linked workflow paper: Steinbaugh MJ, Pantano L, Kirchner RD, Barrera V, Chapman BA, Piper ME, Mistry M, Khetani RS, Rutherford KD, Hoffman O, Hutchinson JN, Ho Sui SJ. (2018).

I have gone through all of the QC steps and created the plots no problem, but now when I move on to the Differential Expression Analysis using DESeq2 I keep running into trouble. I keep getting errors like Error in plotMeanAverage(res) : could not find function "plotMeanAverage" and the same happens for the plotVolcano(res). The alphaSummary function was also not found.

Please can you advise on how to proceed as I can't seem to create any DGE plots.

Is there updated documentation on using bcbioRNASeq and doing DGE analysis, or is the 2018 article above the only one?

Thanks.

[mjsteinbaugh]

Hi @VictoriaPatten , I need to revise the F1000 workflow paper a bit, as we've changed a couple of the functions slightly:

• Ensure our DESeqAnalysis package is installed and loaded.
• Use plotMA instead of plotMeanAverage.
• alphaSummary, which works on DESeqDataSet input, is now defined in DESeqAnalysis instead of bcbioRNASeq.

Current DE template is here, for reference: https://github.com/hbc/bcbioRNASeq/blob/master/inst/rmarkdown/templates/02-differential-expression/skeleton/skeleton.Rmd

Recommended installation methods: either use our bioconda recipe or install in R with:

if (!requireNamespace("BiocManager", quietly = TRUE)) {
    install.packages("BiocManager")
}
install.packages(
    pkgs = "bcbioRNASeq",
    repos = c(
        "https://r.acidgenomics.com",
        BiocManager::repositories()
    ),
    dependencies = TRUE
)

[VictoriaPatten]

Thanks @mjsteinbaugh, this is very helpful. I'll install DESeqAnalysis package and hopefully all will come right. Thanks.

[VictoriaPatten]

@mjsteinbaugh I've followed the template you linked above and all the plots worked beautifully until the very end where we get to 'TopTables'. Here I keep getting the error:

I ran rowData on the dds object and this is the output:

Is the issue that I have Ensemble ID's not gene names? How can I get this TopTables function to work?

Thanks very much Victoria

[VictoriaPatten]

@mjsteinbaugh and @roryk Do you have any advice on my above problem? I've followed the template and run it as a script in R and have also used the markdown template in Rstudio and keep getting the same error (attached above) when I reach the Top Tables section.

Can you advise? Thanks. Victoria

[mjsteinbaugh]

@VictoriaPatten I can reproduce this with a minimal DESeq2 example, and will post a bug fix to DESeqAnalysis later today. Will be in touch with an update!

[VictoriaPatten]

@mjsteinbaugh Thanks so much. I can share my bcb.rds file through googledrive, but obviously cannot upload it here on GitHub due to ethics and confidentiality of patient data etc. Can I use mike@steinbaugh.com or do you have a different preferred email?

[mjsteinbaugh]

Not necessary, I think I made a minimal test dataset that covers your issue. Thanks though!

[mjsteinbaugh]

OK try installing DESeqAnalysis 0.6.1 and see if that resolves your issue:

if (!requireNamespace("BiocManager", quietly = TRUE)) {
    install.packages("BiocManager")
}
install.packages(
    pkgs = "DESeqAnalysis",
    repos = c(
        "https://r.acidgenomics.com",
        BiocManager::repositories()
    ),
    dependencies = TRUE
)

[VictoriaPatten]

@mjsteinbaugh I had already installed this. I've just reinstalled and run it again now to check and I am still getting the exact same error.

[mjsteinbaugh]

What does packageVersion("DESeqAnalysis") return? Make sure to rerun in a clean R session. That error message was only defined in the preivous version.

[VictoriaPatten]

@mjsteinbaugh Okay I think I see the issue. I only have R/4.1.1 at the moment and DESeqAnalysis will only load V0.5.0. I will get the system admin to upgrade to R/4.2.0 and re install. Hopefully that will do the trick! Will be in touch if I hit any more stumbling blocks. Thank you. P.S. I am very new to bioinformatics, so your help is greatly appreciated.

[mjsteinbaugh]

Ah OK yeah that's it, try upgrading to R 4.2 / Bioconductor 3.15. There are a number of changes in Bioconductor 3.15 that make it difficult to continue providing compatiblity support with the R 4.1 release series.

PS You can check your Bioconductor version with:

BiocManager::version()

[VictoriaPatten]

Hi @mjsteinbaugh, the version upgrade seemed to help but now I get a new error after the Top Table section, right at the end of the template. Lines 368-372

→ Exporting <DESeqDataSet> to 'data'.
→ Exporting <SummarizedExperiment> to '/home/hocine/my-projects/victoria/rnaseq/results/2022-05-19/differential-expression/data'.
→ Exporting assays "counts", "normalized", "fpkm" to '/home/hocine/my-projects/victoria/rnaseq/results/2022-05-19/differential-expression/data'.
Quitting from lines 368-372 (bcb_template.Rmd)
Error in `fileExt()`:
! unused argument (pattern = compressExtPattern)
Backtrace:
  1. BiocIO::export(object = deseq, con = params[["results_dir"]])
  2. DESeqAnalysis::export(object = deseq, con = params[["results_dir"]])
  3. DESeqAnalysis .local(object, con, format, ...)
  5. DESeqAnalysis::export(...)
  6. DESeqAnalysis .local(object, con, format, ...)
     ...
 10. AcidExperiment:::.exportAssays(...)
 12. base::lapply(...)
 13. AcidExperiment FUN(X[[i]], ...)
 15. pipette::export(...)
 16. pipette .local(object, con, format, ...

I don't know how to fix this error as it seems it is related to the source code? Thanks

[mjsteinbaugh]

It sounds like you may need to upgrade some other packages. What returns when you run this?

BiocManager::valid(site_repository = "https://r.acidgenomics.com")

[VictoriaPatten]

@mjsteinbaugh that returns: [1] TRUE

[mjsteinbaugh]

OK great, what do these return on your machine?

vapply(
    X = c(
        "AcidBase",
        "AcidExperiment",
        "AcidGenerics",
        "DESeqAnalysis",
        "basejump",
        "bcbioBase",
        "bcbioRNASeq",
        "goalie",
        "pipette",
        "syntactic"
    ),
    FUN = function(x) {
        as.character(packageVersion(x))
    },
    FUN.VALUE = character(1L)
)

print(AcidBase::fileExt)

[VictoriaPatten]

> vapply(
+     X = c(
+         "AcidBase",
+         "AcidExperiment",
+         "AcidGenerics",
+         "DESeqAnalysis",
+         "basejump",
+         "bcbioBase",
+         "bcbioRNASeq",
+         "goalie",
+         "pipette",
+         "syntactic"
+     ),
+     FUN = function(x) {
+         as.character(packageVersion(x))
+     },
+     FUN.VALUE = character(1L)
+ )
      AcidBase AcidExperiment   AcidGenerics  DESeqAnalysis       basejump 
       "0.6.1"        "0.4.0"        "0.6.2"        "0.6.1"       "0.16.0" 
     bcbioBase    bcbioRNASeq         goalie        pipette      syntactic 
       "0.8.0"        "0.5.0"        "0.6.3"        "0.9.1"        "0.6.0" 
> 
> print(AcidBase::fileExt)
function (path) 
{
    pos <- regexpr(pattern = extPattern, text = path, ignore.case = TRUE)
    ifelse(test = pos > -1L, yes = substring(text = path, first = pos + 
        1L), no = NA_character_)
}
<bytecode: 0x55a2745474c0>
<environment: namespace:AcidBase>

[mjsteinbaugh]

Ah OK I need to push a minor fix to r-acidbase. Thanks, this is helpful!

[mjsteinbaugh]

This should be resolved with AcidBase v0.6.2. Try updating with:

if (![requireNamespace](https://rdrr.io/r/base/ns-load.html)("BiocManager", quietly = TRUE)) {
    [install.packages](https://rdrr.io/r/utils/install.packages.html)("BiocManager")
}
[install.packages](https://rdrr.io/r/utils/install.packages.html)(
    pkgs = "AcidBase",
    repos = [c](https://rdrr.io/r/base/c.html)(
        "https://r.acidgenomics.com",
        BiocManager::[repositories](https://rdrr.io/pkg/BiocManager/man/repositories.html)()
    ),
    dependencies = TRUE
)

edit: this got messed up from copy and paste from pkgdown website, see below

[VictoriaPatten]

When I try to run this in R/4.2.0 I get:

> if (![requireNamespace](https://rdrr.io/r/base/ns-load.html)("BiocManager", quietly = TRUE)) {
Error: unexpected '[' in "if (!["
>     [install.packages](https://rdrr.io/r/utils/install.packages.html)("BiocManager")
Error: unexpected '[' in "    ["
> }
Error: unexpected '}' in "}"
> [install.packages](https://rdrr.io/r/utils/install.packages.html)(
Error: unexpected '[' in "["
>     pkgs = "AcidBase",
Error: unexpected ',' in "    pkgs = "AcidBase","
>     repos = [c](https://rdrr.io/r/base/c.html)(
Error: unexpected '[' in "    repos = ["
>         "https://r.acidgenomics.com",
Error: unexpected ',' in "        "https://r.acidgenomics.com","
>         BiocManager::[repositories](https://rdrr.io/pkg/BiocManager/man/repositories.html)()
Error: unexpected '[' in "        BiocManager::["
>     ),
Error: unexpected ')' in "    )"
>     dependencies = TRUE
> )

[mjsteinbaugh]

Sorry about that, here's the correct code:

if (!requireNamespace("BiocManager", quietly = TRUE)) {
    install.packages("BiocManager")
}
install.packages(
    pkgs = "AcidBase",
    repos = c(
        "https://r.acidgenomics.com",
        BiocManager::repositories()
    ),
    dependencies = TRUE
)

[VictoriaPatten]

@mjsteinbaugh Thanks, this seemed to work. When looking at the table output though, some files seem to have the ENS ID's but the top tables (Up-, down-, all-, both.csv) do not have the Ensemble gene ID's. Is this something that can be fixed? Currently I don't see how to identify which genes are differentially expressed in these tables as the rows don't link to specific genes. Also, can row names be added to that tables?

Many thanks. Victoria

[mjsteinbaugh]

Ah good point, I'll work on fixing that

[VictoriaPatten]

Hi @mjsteinbaugh, I know this is a lot for you to look into, but while you are fixing the issue above with the ENS ID's, I've also run into an error with the Functional Analysis template (on RStudio).

Quitting from lines 96-136 (functional_analysis.Rmd) 
Error in `str_match()`:
! could not find function "str_match"
Execution halted

And when my colleague ran the template on his machine he got a bit further (somehow?) but then got the error:

Quitting from lines 177-196 (functional_analysis.Rmd) 
Error in `alert()`:
! could not find function "alert"
Execution halted

Then another thing, because I had done the QC steps from your F1000 paper, not using the rmakdown template, I decided to now go back and re-do the QC using the template, and I'm afraid, more errors:

Quitting from lines 223-228 (qc.Rmd) 
Error:
! Assert failure.
[1] hasLength(keep) is not TRUE.
Cause: `keep` has length 0.
Backtrace:
     ▆
  1. ├─BiocIO::export(...)
  2. └─AcidExperiment::export(...)
  3.   └─AcidExperiment .local(object, con, format, ...)
  4.     └─AcidExperiment:::.exportRowData(...)
  5.       ├─AcidGenerics::atomize(data)
  6.       └─pipette::atomize(data)
  7.         └─pipette .local(object, ...)
  8.           ├─AcidGenerics::atomize(object)
  9.           └─pipette::atomize(object)
 10.             └─pipette .local(object, ...)
 11.               └─goalie::assert(hasLength(keep))
 12.                 └─AcidCLI stop(...)
 13.                   └─cli::cli_abort(x, call = NULL)
 14.                     └─rlang::abort(message, ..., call = call, use_cli_format = TRUE)

Execution halted

I'm terribly sorry for taking up so much of your time with these errors. Are you able to fix them? Thank you so much for your help.

[mjsteinbaugh]

Hi @VictoriaPatten I'm working on some package updates today, will send an update soon

[mjsteinbaugh]

PS alert is from AcidCLI, and str_match is from stringr

[VictoriaPatten]

Hi @mjsteinbaugh have you managed to do those package updates to add ENS ID's to the top tables and also fix the issues in the QC template?

Also, I have both stringr and AcidCLI packages installed, but I am still getting this error with the functional analysis template.

Quitting from lines 96-136 (functional_analysis.Rmd) 
Error in `str_match()`:
! could not find function "str_match"
Execution halted

Hope to hear from you soon. Thanks

[mjsteinbaugh]

Hi @VictoriaPatten , I'm working on an update that I'll push soon!

[VictoriaPatten]

Hi @mjsteinbaugh, have you been able to finish the updates you were working on?

We ran the functional analysis again and commented out the alert function and now there is another error related to the organism.

The organism of dds is Homo Sapiens (with capital S). In line 297 for the RankedList we get the error:

! Assert failure.
[1] identical(organism(dds), "Homo sapiens") is not TRUE.

I don't know which script generated the Homo Sapiens (with capital S), how can I fix this?

Thanks so much Victoria

[mjsteinbaugh]

Hi @VictoriaPatten , I'm working on finishing a couple of updates to bcbioRNASeq, and will try to get those out today if possible. I'll post a link to an updated R Markdown template that you can try out with your dataset.

[VictoriaPatten]

Hi @mjsteinbaugh sorry to keep bothering you about this, I'm just on quite a tight time schedule as I have to submit my thesis in 6 weeks time so I really need to get this analysis completed. Have you managed to do the R Markdown template updates?

[mjsteinbaugh]

Hi @VictoriaPatten apologies for the delay, as I'm on vacation this week. I'll work on getting this fixed ASAP and ping you back!

[VictoriaPatten]

@mjsteinbaugh any update? Are you still on vacation? I'm still running into issues with the functional analysis template and the Ensemble IDs vs Entrez ID for the GSEA and KEGG analysis, and also still the issue with 'Homo sapiens'

Quitting from lines 297-324 (functional_analysis.Rmd) 

Error:
! Assert failure.
[1] identical(organism(dds), "Homo sapiens") is not TRUE.
Backtrace:
    ▆
 1. ├─AcidGSEA::RankedList(object = deseq, value = "stat", keyType = "geneName")
 2. └─AcidGSEA::RankedList(object = deseq, value = "stat", keyType = "geneName")
 3.   └─AcidGSEA .local(object, ...)
 4.     └─goalie::assert(identical(organism(dds), "Homo sapiens"))
 5.       └─AcidCLI stop(...)
 6.         └─cli::cli_abort(x, call = NULL)
 7.           └─rlang::abort(message, ..., call = call, use_cli_format = TRUE)
Warning message:
ggrepel: 1 unlabeled data points (too many overlaps). Consider increasing max.overlaps 
Execution halted

Hope to hear from you soon. Victoria

[mjsteinbaugh]

Hi @VictoriaPatten just got back, I'll work on package updates this week. Apologies for the delay!

[mjsteinbaugh]

Hi @VictoriaPatten , quick followup -- I'm blocking out some time tomorrow and Wednesday to to finish this update.

[VictoriaPatten]

@mjsteinbaugh Thank you, thank you, thank you!! I am beyond grateful.

[VictoriaPatten]

Hi @mjsteinbaugh its been two weeks... were you able to finish the update? I am still getting errors with the functional analysis with ensemble and entrezID issues. I'm a bit desperate with this now.

[mjsteinbaugh]

Apologies for the delay, working on it this morning. I'm free to do a Zoom chat this week as well, which may speed up the troubleshooting

[mjsteinbaugh]

OK I can confirm that the template is working with my example dataset. I'm working on pushing a couple of updates to AcidGSEA and bcbioRNASeq that should resolve your issues.

[VictoriaPatten]

Thats great! I'll run my dataset through and hopefully it will work too. Does your dataset use Ensembl ID's? Perhaps if I cant get things working by the end of the week we could do that zoom chat on Friday.

[VictoriaPatten]

Hi @mjsteinbaugh I am still getting this error:

Quitting from lines 374-378 (functional_analysis.Rmd) 
Error:
! Assert failure.
[1] isSubset(x = c("geneId", "entrezId"), y = colnames(mcols(object))) is not
TRUE.
Cause: `c("geneId", "entrezId")` has elements not in `colnames(mcols(object))`:
entrezId
Backtrace:
    ▆
 1. ├─AcidGenerics::Ensembl2Entrez(object = dds, format = "1:1")
 2. └─AcidExperiment::Ensembl2Entrez(object = dds, format = "1:1")
 3.   ├─AcidGenerics::Ensembl2Entrez(rowRanges(object), ...)
 4.   └─AcidGenomes::Ensembl2Entrez(rowRanges(object), ...)
 5.     └─AcidGenomes .local(object, ...)
 6.       └─goalie::assert(isSubset(x = c("geneId", "entrezId"), y = colnames(mcols(object))))
 7.         └─AcidCLI stop(...)
 8.           └─cli::cli_abort(x, call = NULL)
 9.             └─rlang::abort(message, ..., call = call, use_cli_format = TRUE)
Warning messages:
1: ggrepel: 1 unlabeled data points (too many overlaps). Consider increasing max.overlaps 
2: ggrepel: 1 unlabeled data points (too many overlaps). Consider increasing max.overlaps 
Execution halted

Does this have something to do with the AcidGSEA update you were going to push? How should I proceed?

[mjsteinbaugh]

OK I think the issue with your dataset is that you don't have enough gene metadata in the rowRanges / rowData of your original bcbioRNASeq object to correctly run the clusterProfiler enrichment analysis template with Ensembl-to-Entrez identifier interconversion. This requires setting the organism and ensemblRelease arguments during the main bcbioRNASeq() generator call. Refer to the bcbioRNASeq documentation for details.

Don't worry, you don't need to re-run the analysis from scratch. We just need to tweak your DESeqAnalysis object to contain the sufficient Entrez identifier metadata so you can successfully run Ensembl-to-Entrez conversion.

Assuming you have a DESeqAnalysis object, generated by our differential expression R Markdown template.

## DESeqAnalysis object
## https://r.acidgenomics.com/packages/deseqanalysis/
## Generated with code like:
## > deseq <- DESeqAnalysis::DESeqAnalysis(
## >   data = <DESeqDataSet>,
## >   transform = <DESeqTransform>,
## >   results = <list of unshrunken results>,
## >   lfcShrink = <list of shrunken results>
## )

If you need to assign rowRanges, and your data was processed using Ensembl as the reference genome, here's the way. This gets automatically generated by bcbioRNASeq if you use organism during the main bcbioRNASeq() call.

suppressPackageStartupMessages({
    library(SummarizedExperiment)
    library(AcidGenomes)
})
## NOTE Update the Ensembl release to match the version used in your analysis (e.g. 106).
rowRanges <- makeGRangesFromEnsembl(
    organism = "Homo sapiens",
    release = 106L
)
c("geneId", "entrezId") %in% colnames(mcols(rowRanges))
## If this step errors, you likely have a mismatch with your reference genome.
stopifnot(all(names(rowRanges) %in% rownames(deseq@data)))
## Reorder the GenomicRanges to match your DESeqDataSet.
rowRanges <- rowRanges[rownames(deseq@data)]
## Slot the updated ranges back into your DESeqDataSet and DESeqTransform objects.
rowRanges(deseq@data) <- rowRanges
rowRanges(deseq@transform) <- rowRanges
validObject(deseq)

Now your object should contain the desired geneId and entrezId columns in rowRanges metadata required for these GSEA steps. Give that a shot!

[mjsteinbaugh]

PS a draft update of the functional analysis template is available here, for reference:

https://github.com/acidgenomics/r-bcbiornaseq/blob/develop/inst/rmarkdown/templates/03-functional-analysis/skeleton/skeleton.Rmd

[VictoriaPatten]

@mjsteinbaugh I couldnt get the tweaks to work so I ended up re-running the entire analysis from scratch creating the original bcb.rds object setting organism and ensemblRelease. I then ran this object through the DGE template to get a new deseq.rds object and ran this through the original Functional Analysis template.

When keytype = geneName, this error:

preparing geneSet collections...
--> Expected input gene ID: ENSG00000130270,ENSG00000188488,ENSG00000133794,ENSG00000070770,ENSG00000158955,ENSG00000135624
Quitting from lines 329-363 (functional_analysis.Rmd) 

Error in `check_gene_id()`:
! --> No gene can be mapped....
Backtrace:
 1. clusterProfiler::gseGO(...)
 2. DOSE:::GSEA_internal(...)
 3. DOSE .GSEA(...)
 4. DOSE:::check_gene_id(geneList, geneSets)
Warning messages:
1: ggrepel: 1 unlabeled data points (too many overlaps). Consider increasing max.overlaps 
2: ggrepel: 1 unlabeled data points (too many overlaps). Consider increasing max.overlaps 
Execution halted

I also tried keytype = ensemblId, or entrezId, or geneId, but all gave erorrs from line 299 keytype = entrezId.... Error: Failed to map at least 50% of gene identifiers. keytype = ensemblId...Error in as.vector(x, mode = "character"): no method for coercing this S4 class to a vector keytype = geneId ...Error in as.vector(x, mode = "character"): no method for coercing this S4 class to a vector

[mjsteinbaugh]

Thanks I'm taking a look. Can you post or email a link to the deseq object used? That will make it much easier to debug. This may be an issue with clusterProfiler and DOSE.

[mjsteinbaugh]

It's this part, correct?

gseaGo <- clusterProfiler::gseGO(
    geneList = rankedListEnsembl,
    ont = params[["go_class"]],
    OrgDb = get(orgDb, envir = asNamespace(orgDb)),
    keyType = "ENSEMBL",
    pvalueCutoff = 0.05
)

Can you print what rankedListEnsembl looks like on your machine?

[mjsteinbaugh]

Here's a minimal example that should work:

suppressPackageStartupMessages({
    library(DESeqAnalysis) # 0.6.4
    library(AcidGSEA) # 0.8.4
})
data(deseq, package = "DESeqAnalysis")
ensembl <- RankedList(
    object = deseq,
    value = "stat",
    keyType = "ensemblId"
)
print(ensembl)
## RankedList of length 2
## names(2): condition_B_vs_A treatment_D_vs_C
print(head(ensembl[[1L]]))
# ENSG00000004848 ENSG00000008118 ENSG00000015133 ENSG00000011451 ENSG00000005100
#        4.681582        4.532673        4.309957        4.048786        3.993137
# ENSG00000007545
#        3.916430
entrez <- RankedList(
    object = deseq,
    value = "stat",
    keyType = "entrezId"
)
print(entrez)
## RankedList of length 2
## names(2): condition_B_vs_A treatment_D_vs_C
print(head(entrez[[1L]]))
##   170302    57172   440193    58525    56919    57585
## 4.681582 4.532673 4.309957 4.048786 3.993137 3.916430

[VictoriaPatten]

@mjsteinbaugh I sent you a googledrive link to my deseq.rds object (mike@steinbaugh.com ?). Will test out the above shortly.

[mjsteinbaugh]

Perfect thanks! I'll try debugging with your dataset later today

[VictoriaPatten]

@mjsteinbaugh not sure if what I'm doing is correct. I tried to run the above minimal example you sent added in to the template at the start of the GO GSEA Analysis section, and get this error:

ℹ 97.6% of genes mapped to primary seqnames (57323 / 58735).

ℹ 97.6% of genes mapped to primary seqnames (57323 / 58735).
ℹ Mapping 43.5% of genes to `entrezId` (24936 / 57323).

Quitting from lines 299-347 (functional_analysis.Rmd) 

Error:
! BiocParallel errors
  1 remote errors, element index: 1
  0 unevaluated and other errors
  first remote error:
Error: Failed to map at least 50% of gene identifiers.
Backtrace:
  1. rmarkdown::render(...)
  2. knitr::knit(knit_input, knit_output, envir = envir, quiet = quiet)
  3. knitr:::process_file(text, output)
  6. knitr:::process_group.block(group)
  7. knitr:::call_block(x)
     ...
 12. knitr (local) evaluate(...)
 13. evaluate::evaluate(...)
 14. evaluate:::evaluate_call(...)
 19. evaluate:::eval_with_user_handlers(expr, envir, enclos, user_handlers)
 23. base::eval(expr, envir, enclos)
Warning messages:
1: ggrepel: 1 unlabeled data points (too many overlaps). Consider increasing max.overlaps 
2: ggrepel: 1 unlabeled data points (too many overlaps). Consider increasing max.overlaps 
Execution halted

The Package versions are all up to date