Repository navigation
invalid DNA letters #476
Description
Hi,
I used Bambu to quantify transcripts based on a modified reference genome (because I did a allelic specific analysis and modified the official assembly). I faced the problem as below:
Error: BiocParallel errors
1 remote errors, element index: 1
0 unevaluated and other errors
first remote error:
Error in value[3L]: record 686186 (16:79451803-79451804) contains invalid DNA letters
file: .../hap2.fa
Execution halted
I have gone through the websites and someone mentioned it being due to R object GRanges. Would you please give a clue how to bypass or solve this issue?
Best regards
Sang
Activity
cying111 commentedon Apr 3, 2025
Hi @Geil625 ,
Thanks for reporting this to us. it looks like the processing of the fa file is a bit problematic here.
For me to understand a bit more of the context here, could you share your command on how you use bambu, more specifically how the fa file is provided, as a string variable or as a FaFile?
Also, which OS system are you running bambu on?
Looking forward to your reply
Thanks
Ying
Geil625 commentedon Apr 3, 2025
Hi,
I have tried to use branch HongYhong_fix as a similar issue was reported in #402. Fortunately, my issue could be resolved but also a bunch of warnings were reported as:
WARNING - If you change the number of cores (ncore) between Bambu runs and there is no progress please restart your R session to resolve the issue that originates from the XGboost package.
--- Start generating read class files ---
Warning: unable to access index for repository https://bioconductor.org/packages/3.20/data/annotation/src/contrib:
cannot open URL 'https://bioconductor.org/packages/3.20/data/annotation/src/contrib/PACKAGES'
[13:14:52] WARNING: src/learner.cc:553:
If you are loading a serialized model (like pickle in Python, RDS in R) generated by
older XGBoost, please export the model by calling
Booster.save_modelfrom that versionfirst, then load it back in current version. See:
https://xgboost.readthedocs.io/en/latest/tutorials/saving_model.html
for more details about differences between saving model and serializing.
[13:14:52] WARNING: src/learner.cc:553:
If you are loading a serialized model (like pickle in Python, RDS in R) generated by
older XGBoost, please export the model by calling
Booster.save_modelfrom that versionfirst, then load it back in current version. See:
https://xgboost.readthedocs.io/en/latest/tutorials/saving_model.html
for more details about differences between saving model and serializing.
Detected 2 warnings across the samples during read class construction. Access warnings with metadata(bambuOutput)$warnings
--- Start extending annotations ---
Using a novel discovery rate (NDR) of: 0.229
--- Start isoform quantification ---
--- Finished running Bambu ---
Warning message:
There was 1 warning in
filter().ℹ In argument:
newGeneId == min(newGeneId).Caused by warning in
min():! no non-missing arguments to min; returning Inf
I wonder if these warnings would impact the quantification result or not?
Acutually I was doing allelic specific analysis using two haplotype-aware sequences with following codes:
rm(list=ls())
library(bambu)
options(BIOCONDUCTOR_CONFIG_FILE = "bambu/config.yaml")
gtf.file <- 'modified.gtf'
annotations <- prepareAnnotations(gtf.file)
args <- commandArgs(T)
work_dir <- args[1]
ppn <- args[2]
res <- bambu(reads = paste0(work_dir,'/hap1.bam'), annotations = annotations, genome = paste0(work_dir,'/hap1.fa'),ncore = ppn)
writeBambuOutput(res, path = paste0(work_dir,'/hap1_output'))
res <- bambu(reads = paste0(work_dir,'/hap2.bam'), annotations = annotations, genome = paste0(work_dir,'/hap2.fa'),ncore = ppn)
writeBambuOutput(res, path = paste0(work_dir,'/hap2_output'))
The code was running well for hap1 result but reporting the error I displayed at the very beginning for hap2. However, what is weird is if I marked hap1 code and only running hap2 code as:
# res <- bambu(reads = paste0(work_dir,'/hap1.bam'), annotations = annotations, genome =paste0(work_dir,'/hap1.fa'),ncore = ppn)
# writeBambuOutput(res, path = paste0(work_dir,'/hap1_output'))
There was no error reported but just the warnings I showed above! Now I have to separately run hap1 and hap2 in different scripts .
By the way, I compared the devel and HongYhong_fix branches results regarding transcripts quantification, and found there is a bit difference. I wonder the core code or algorithm between the two branches is same?
Cheers,
Sang
Geil625 commentedon Apr 3, 2025
And also ppn <- args[2] seems made ppn as a character, so should I made it as an integer using as.integer(ppn)?
I am running Bambu on Linux system.
cying111 commentedon Apr 3, 2025
Hi,
I am glad that it works out when using HongYhong_fix branch. You can ignore the warnings for now and they won't affect the quantification results.
Regarding the need to run the two haplotypes separately, it could be that the annotation object changes after running Bambu for hap1.bam. A possible solution is to initialize the annotations beforehand:
bambuAnnotations <- prepareAnnotations(gtf.file)and later use bambuAnnotations when calling bambu.As for the difference in quantification results between HongYhong_fix and devel, how large is the discrepancy? The devel branch has had some recent bug fixes that may not be present in HongYhong_fix, which could explain the variation.
Let me know what you find!
Thanks
Ying
Geil625 commentedon Apr 3, 2025
Thanks for the quick turnaround. Using one sample from my data as an example, the Pearson correlation between CPM of official annotated transcripts quantified by the two branches is 0.998. For the novel transcripts, the number differred (1779 vs. 1753) and could not be compared.
The package version for devel and HongYhong_fix are respectively 3.4.0 and 3.3.5.
By the way, the ncore used in the main function must be numeric or character is also OK?
Geil625 commentedon Apr 3, 2025
Unfortunately, using HongYhong_fix just could addresss one sample. For all other samples the error is also reported like:
Error: BiocParallel errors
1 remote errors, element index: 1
0 unevaluated and other errors
first remote error:
Error in value[3L]: record 686186 (16:79451803-79451804) contains invalid DNA letters
file: .../hap2.fa
Execution halted
I have gone through the exact letters in chr16: 79451803-79451804 in hap2.fa and found no problem (all letters are A C G T a c g t).