error in evaluating the argument 'args' in selecting a method for function 'do.call'

[shangguandong1996]

Description of the bug

Hi, I am running hicar on my Arabidopsis thaliana hicar data.
But it seems that there is a error about ATACSeqQC

nextflow run /data5/sgd_data/biosoft/nf-core/nf-core-hicar-dev/workflow --input samplesheet.csv --outdir ../result --fasta ~/reference/genome/TAIR10/Athaliana_with_chr.fa --gtf ~/reference/annoation/Athaliana/Araport11/Araport11_GFF3_genes_transposons.201606.gtf --macs_gsize 1.1e8 -profile singularity --juicer_tools_jar /data5/sgd_data/biosoft/juicer_tools_2.13.06.jar --merge_map_py_source /data5/sgd_data/biosoft/ijuric/merge_map.py --feature_frag2bin_source /data5/sgd_data/biosoft/ijuric/feature_frag2bin.py --make_maps_runfile_source /data5/sgd_data/biosoft/ijuric/make_maps_runfile.py

Command used and terminal output

-[nf-core/hicar] Pipeline completed with errors-
Error executing process > 'NFCORE_HICAR:HICAR:ATAC_PEAK:ATACQC'

Caused by:
  Process `NFCORE_HICAR:HICAR:ATAC_PEAK:ATACQC` terminated with an error exit status (1)

Command executed:

  #!/usr/bin/env Rscript
  
  #######################################################################
  #######################################################################
  ## Created on Sept. 10, 2021 stats for ATAC reads
  ## Copyright (c) 2021 Jianhong Ou ([email protected])
  #######################################################################
  #######################################################################
  pkgs <- c("rtracklayer", "GenomicFeatures", "ATACseqQC")
  versions <- c("NFCORE_HICAR:HICAR:ATAC_PEAK:ATACQC:")
  for(pkg in pkgs){
      # load library
      library(pkg, character.only=TRUE)
      # parepare for versions.yml
      versions <- c(versions,
          paste0("    ", pkg, ": ", as.character(packageVersion(pkg))))
  }
  writeLines(versions, "versions.yml") # wirte versions.yml
  args <- strsplit("--pattern .merged.ATAC.bed.gz", "\\s+")[[1]]
  parse_args <- function(options, args){
      out <- lapply(options, function(.ele){
          if(any(.ele[-3] %in% args)){
              if(.ele[3]=="logical"){
                  TRUE
              }else{
                  id <- which(args %in% .ele[-3])[1]
                  x <- args[id+1]
                  mode(x) <- .ele[3]
                  x
              }
          }
      })
  }
  option_list <- list("pattern"=c("--pattern", "-p", "character"))
  opt <- parse_args(option_list, args)
  pattern <- opt[['pattern']] #"merged.ATAC.bed.gz"
  gtf <- "Araport11_GFF3_genes_transposons.201606.gtf"
  
  readsFiles <- dir(".", pattern) ## postfix from mergedreads.nf
  peaksFiles <- dir(".", "narrowPeak|broadPeak") ## output from macs2
  
  names(readsFiles) <- sub(pattern, "", readsFiles)
  names(peaksFiles) <- sub("_peaks.*Peak", "", peaksFiles)
  
  N <- intersect(names(readsFiles), names(peaksFiles))
  if(length(N)==0){
      if(length(peaksFiles)==1){ # for user defined peaks.
          peaksFiles <- rep(peaksFiles, length(readsFiles))
          names(peaksFiles) <- names(readsFiles)
          N <- intersect(names(readsFiles), names(peaksFiles))
      }
      stopifnot("no peak and signal pairs"=length(N)>0)
  }
  
  peaksFiles <- peaksFiles[N]
  readsFiles <- readsFiles[N]
  
  ## import reads
  readls <- lapply(readsFiles, function(f){
      reads <- read.table(f, colClasses=c("character", "integer", "NULL",
                              "NULL", "NULL", "character"))
      reads <- GRanges(reads[, 1],
                      IRanges(as.numeric(reads[, 2])+1, as.numeric(reads[, 2])+150),
                      strand = reads[, 3])
  })
  
  ## import peaks
  peakls <- lapply(peaksFiles, import)
  
  txdb <- makeTxDbFromGFF(gtf) #for TSS
  txs <- exons(txdb)
  
  stats <- mapply(peakls, readls, FUN = function(peaks, reads){
      ## calculate FRiP score (fraction of reads in peaks), must over 1%
      readsInPeaks <- countOverlaps(peaks, reads)
      FRiP <- 100*sum(readsInPeaks)/length(reads)
  
      reads <- as(reads, "GAlignments")
      ## calculate Transcription Start Site Enrichment Score
      tsse <- TSSEscore(reads, txs)
  
      ## Promoter/Transcript body (PT) score
      pt <- PTscore(reads, txs)
      pt <- pt[!is.na(pt$promoter) & !is.na(pt$transcriptBody) & !is.na(pt$PT_score)]
      pt <- pt[(pt$promoter>0 | pt$transcriptBody>0) & pt$PT_score!=0]
      promoterEnriched <- table(pt$PT_score>0)
      names(promoterEnriched) <-
          c("FALSE"="bodyEnrich", "TRUE"="promoterEnrich")[names(promoterEnriched)]
      promoterEnriched <-
          c(promoterEnriched,
              prop.test=prop.test(cbind(table(pt$PT_score>0), c(50, 50)))$p.value)
      list(tsse_FRiP = c(TSSEscore=tsse$TSSEscore, FRiP=FRiP, promoterEnriched),
          tsseValues = tsse$values)
  }, SIMPLIFY = FALSE)
  
  ## FRiP, TSSEscore table
  tsse_FRiP <- do.call(rbind, lapply(stats, function(.ele) .ele$tsse_FRiP))
  write.csv(tsse_FRiP, "TSSEscore_FRiP.csv")
  
  ## for TSSEscore to TSS plots
  tsse <- do.call(rbind, lapply(stats, function(.ele) .ele$tsseValues))
  if(ncol(tsse)==20){
      colnames(tsse) <- 100*(-9:10-.5)
      write.csv(tsse, "aggregateTSSEscoreToTSS.csv")
  }

Command exit status:
  1

Command output:
  (empty)

Command error:
  The following objects are masked from ‘package:stats’:
  
      IQR, mad, sd, var, xtabs
  
  The following objects are masked from ‘package:base’:
  
      anyDuplicated, append, as.data.frame, basename, cbind, colnames,
      dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,
      grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
      order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
      rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,
      union, unique, unsplit, which.max, which.min
  
  Loading required package: S4Vectors
  
  Attaching package: ‘S4Vectors’
  
  The following objects are masked from ‘package:base’:
  
      expand.grid, I, unname
  
  Loading required package: IRanges
  Loading required package: GenomeInfoDb
  Loading required package: AnnotationDbi
  Loading required package: Biobase
  Welcome to Bioconductor
  
      Vignettes contain introductory material; view with
      'browseVignettes()'. To cite Bioconductor, see
      'citation("Biobase")', and for packages 'citation("pkgname")'.
  
  Import genomic features from the file as a GRanges object ... OK
  Prepare the 'metadata' data frame ... OK
  Make the TxDb object ... OK
  Warning messages:
  1: In .get_cds_IDX(mcols0$type, mcols0$phase) :
    The "phase" metadata column contains non-NA values for features of type
    stop_codon, exon. This information was ignored.
  2: In .reject_transcripts(bad_tx, because) :
    The following transcripts were dropped because they have stop codons
    that cannot be mapped to an exon: AT1G07320.4, AT1G18180.2,
    AT1G30230.1, AT1G36380.1, AT1G52415.1, AT1G52940.1, AT2G18110.1,
    AT2G35130.1, AT2G35130.3, AT2G39050.1, AT3G13445.1, AT3G17660.1,
    AT3G17660.3, AT3G52070.2, AT3G54680.1, AT3G59450.2, AT4G13730.2,
    AT4G17730.1, AT4G39620.2, AT5G01520.2, AT5G22794.1, AT5G27710.2,
    AT5G45240.2
  Error in h(simpleError(msg, call)) : 
    error in evaluating the argument 'args' in selecting a method for function 'do.call': vector: cannot make a vector of mode 'Rle'.
  Calls: mapply ... emptyOpsReturnValue -> do.call -> .handleSimpleError -> h
  Execution halted

Work dir:
  /data5/sgd_data/newProject/202207/HiCAR_SIM3d_202207/rawdata/work/15/6993953699ef43f7501a50c93a6eb1

Tip: view the complete command output by changing to the process work dir and entering the command `cat .command.out`

Relevant files

No response

System information

No response

[shangguandong1996]

sorry, it is my mistake.
Because my inital fa is Chr1,Chr2, which is not UCSC name. So I have to change the fa into chr1, 2. But I do not change the gtf seqnames. Once I use the chrXX gtf, it run sucessfully.