How to use a fimo BED output as GRanges object in factorFootprints from ATACseqQC
1
1
Entering edit mode
pedro.rotcha ▴ 10
@pedrorotcha-24078
Last seen 3.3 years ago

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?

ATACseqQC • 1.1k views
ADD COMMENT
0
Entering edit mode

Hi pedro.rotcha, You did not supply the seqlenths for chrX and chrY. You may want to try this:

seqlev <- seqlevels(bindingSites)
seqinfo(bindingSites) <- seqinfo(genome)[seqlevels(bindingSites)]
factorFootprints(shiftedBamfile, pfm=GATA6[[1]],
             genome=genome, bindingSites=bindingSites,
             min.score="90%", seqlev=seqlev,
             upstream=100, downstream=100)

Let me known if you still have trouble in running the code.

jianhong.

ADD REPLY
0
Entering edit mode

Perfect! thanks so much

ADD REPLY
0
Entering edit mode
Julie Zhu ★ 4.3k
@julie-zhu-3596
Last seen 4 months ago
United States

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

ADD COMMENT
0
Entering edit mode

together with the line:

seqinfo(bindingSites) <- seqinfo(genome)[seqlevels(bindingSites)]

that Jianhong suggested this worked. Thanks

ADD REPLY

Login before adding your answer.

Traffic: 742 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6