Trying to use a fimo output file with a list of motifs for a TF found at that TFs ChIPseq peaks. Bed file example here Did

> library(ATACseqQC)
library(TxDb.Mmusculus.UCSC.mm10.knownGene)
library(Rsamtools)
library(BSgenome.Mmusculus.UCSC.mm10)
library(MotifDb)
library(regioneR)
 file<-("fimo.bed")
seqlev = paste0("chr", c(1:19, "X", "Y"))
genome=Mmusculus
outPath <- "splited"
shiftedBamfile <- file.path(outPath, "shifted.bam")
bindingSites <- toGRanges(file)
bindingSites$name<-NULL
seqlengths(bindingSites) = c(195471971,130694993,122082543,120129022,120421639,124902244,104043685,98207768,94987271,90702639,61431566,182113224,160039680,156508116,151834684,149736546,145441459,129401213,124595110)
 GATA6 <- query(MotifDb, c("GATA6"))
GATA6 <- as.list(GATA6)
 factorFootprints(shiftedBamfile, pfm=GATA6[[1]],
                 genome=genome, bindingSites=bindingSites,
                 min.score="90%", seqlev=seqlev,
                 upstream=100, downstream=100)

factorFootprints outputs

> Error in asMethod(object) :
  cannot create a GRanges object from a Seqinfo object with NA seqlengths

This is how bindingSites looks like:

> GRanges object with 26898 ranges and 1 metadata column:
          seqnames              ranges strand |     score
             <Rle>           <IRanges>  <Rle> | <numeric>
      [1]     chr1     3446064-3446079      + |      36.4
      [2]     chr1     3446360-3446375      - |      32.8
      [3]     chr1     3446020-3446035      + |      30.4
      [4]     chr1     3949250-3949265      - |      50.9
      [5]     chr1     3949463-3949478      + |      43.9
      ...      ...                 ...    ... .       ...
  [26894]     chr9 123939661-123939676      - |      31.7
  [26895]     chr9 123939543-123939558      - |      31.2
  [26896]     chr9 123939519-123939534      - |      30.9
  [26897]     chr9 123958011-123958026      - |      47.6
  [26898]     chr9 123958371-123958386      - |      33.4
  -------seqinfo: 19 sequences from an unspecified genome

so I guess the question is how to actually generate the bindingSites GRange object?

Hi Pedro,

Please try to set seqlev using the following command. seqlev <- seqlevels(mt)

If this does not solve the issue, please share your bam file as well. Thanks!

Best, Julie