0

我正在尝试使用包 bambu 来量化 bam 文件中的基因计数。我正在使用我大学的 HPC，所以我编写了一个 R 脚本和一个批处理提交文件来启动它。

当脚本运行 bambu 函数时，它会给出以下错误：

Start generating read class files
  |                                                                      |   0%[W::hts_idx_load2] The index file is older than the data file: ./results/minimap2/KD_R1.sorted.bam.bai
[W::hts_idx_load2] The index file is older than the data file: ./results/minimap2/KD_R3.sorted.bam.bai
[W::hts_idx_load2] The index file is older than the data file: ./results/minimap2/WT_R1.sorted.bam.bai
[W::hts_idx_load2] The index file is older than the data file: ./results/minimap2/WT_R2.sorted.bam.bai
  |==================                                                    |  25%
Error: BiocParallel errors
  element index: 1, 2, 3
  first error: cannot open the connection
In addition: Warning message:
stop worker failed:
  attempt to select less than one element in OneIndex 
Execution halted

所以看起来 BiocParallel 不高兴并且无法打开某个连接，但我不知道如何解决这个问题？

这是我的 R 脚本：

#Bambu R script

#load libraries
library(Rsamtools)
library(bambu)

#Creating files
bamFiles<- Rsamtools::BamFileList(c("./results/minimap2/KD_R1.sorted.bam","./results/minimap2/KD_R2.sorted.bam","./results/minimap2/KD_R3.sorted.bam","./results/minimap2/WT_R1.sorted.bam","./results/minimap2/WT_R2.sorted.bam","./results/minimap2/WT_R3.sorted.bam"))
annotation<-prepareAnnotations("./ref_data/Homo_sapiens.GRCh38.104.chr.gtf")
fa.file<-"./ref_data/Homo_sapiens.GRCh38.dna.primary_assembly.fa"

#Running bambu
se<- bambu(reads=bamFiles, annotations=annotation, genome=fa.file,ncore=4)
se
seGene<- transcriptToGeneExpression(se)

#Saving files
save.file<-tempfile(fileext=".gtf")
writeToGTF(rowRanges(se),file=save.file)
save.dir <- tempdir()
writeBambuOutput(se,path=save.fir,prefix="Nanopore_")
writeBambuOutput(seGene,path=save.fir,prefix="Nanopore_")

如果您对为什么会发生这种情况有任何想法，那将非常有帮助！谢谢

4

1 回答 1

0

我认为@Chris 有一个很好的观点。在引擎盖下，bambu 似乎是htslib根据这些警告运行的。虽然它们可能确实只是警告，但我想知道如果您以交互方式运行它，结果会是什么样子。

这个问题现在很难回答，因为它缺少一些信息（文件是什么样的，一个最小的可重现示例等）。但与此同时，这里有一些可能有用的问题来解决这个问题：

  1. 长什么样bamFiles？它是否有正确数量的读取记录？所有这些文件都有非零读取记录吗？有可疑的小吗？
  2. baivs文件上的时间戳是什么bam（例如ls -lh /results/minimap2/）？它们是关于你所期望的还是很不稳定？它们中的任何一个（比如说，./results/minimap2/WT_R2.sorted.bam.bai）都小得奇怪吗？
  3. 当您以交互方式运行它时会发生什么？它在哪里失败？你说是在bambu()通话中，但你怎么知道的？
  4. 当你运行时会发生bambu()什么ncores=1？

这似乎很可能是由于文件问题造成的，并且只有在这biocParallel一步错误才会冒泡到顶部。许多实用程序都有一个令人讨厌的习惯，即乐于接受空文件，但当被要求对空文件执行某些操作时，只会在没有信息性错误消息的情况下令人困惑地失败。

您也可以考虑向开发人员提出问题。

（为什么警告可能只是一个问题：索引文件有时具有类似于以编程方式生成和索引的非常小的对齐文件的时间戳，其中索引步骤几乎是即时的。）

于 2021-08-09T23:40:34.973 回答