use factorFootprints function in ATACseqQC to get DNA footprinting for certain gene region
2
0
Entering edit mode
@aiminatwork-12324
Last seen 5.4 years ago

I am using your ATACseqQC to generate DNA footprinting for certain gene region.

 

I noted that for the function of factorFootprints, there is an argument of “bindingSite” to define gene region, so I did

something like the following:

 

which <- GRanges(seqnames = “chr1", ranges = IRanges(1000,2000))

 

factorFootprints(bamfile, pfm=CTCF[[1]],genome=Hsapiens,bindingSites=which,min.score="95%", seqlev="chr1",upstream=100, downstream=100)

 

But I got the error like:

 

>  factorFootprints(bamfile, pfm=CTCF[[1]],
+                 genome=Hsapiens,bindingSites=which,
+                 min.score="95%", seqlev="chr1",
+                 upstream=100, downstream=100)
Error in factorFootprints(bamfile, pfm = CTCF[[1]], genome = Hsapiens,  : 
  all(!is.na(seqlengths(bindingSites))) is not TRUE

 

I am just wondering if you have an example to show how to use  the argument of “bindingSite” to define gene region, then to get DNA footprinting for the defined gene region.

 

Thank you,

 

Aimin
ATACseqQC • 1.8k views
ADD COMMENT •
0
Entering edit mode

Aimin,

The bindingSites you set might not contain any CTCF binding sites. Please try to run the analysis without setting the parameter bindingSites. By default, the program will perform genome-wide search of the CTCF binding sites. 

Best,

Julie

ADD REPLY •
0
Entering edit mode
Ou, Jianhong ★ 1.3k
@ou-jianhong-4539
Last seen 5 days ago
United States

Hi Aimin,

Thank you for trying ATACseqQC to plot the TF footprints. You can have a try following code:

which <- GRanges(seqnames = "chr1", ranges = IRanges(seq(10000, 20000, by=1000), width = 10))
seqinfo(which) <- seqinfo(Hsapiens)
factorFootprints(bamfile, pfm=CTCF[[1]],
                 genome=Hsapiens,bindingSites=which,
                 min.score="95%", seqlev="chr1",
                 upstream=100, downstream=100)

Let me know if you have any trouble in running this sample.

The idea is that the program will double check the bindingSites genome size.

 

Jianhong.
ADD COMMENT •
0
Entering edit mode
Thank you, Yes, I tried, but I got same error as before. Aimin > which <- GRanges(seqnames = "chr1", ranges = IRanges(seq(10000, 20000, by=1000), width = 10)) > which GRanges object with 11 ranges and 0 metadata columns: seqnames ranges strand <rle> <iranges> <rle> [1] chr1 10000-10009 * [2] chr1 11000-11009 * [3] chr1 12000-12009 * [4] chr1 13000-13009 * [5] chr1 14000-14009 * [6] chr1 15000-15009 * [7] chr1 16000-16009 * [8] chr1 17000-17009 * [9] chr1 18000-18009 * [10] chr1 19000-19009 * [11] chr1 20000-20009 * ------- seqinfo: 1 sequence from an unspecified genome; no seqlengths > seqinfo(which) <- seqinfo(Hsapiens) > factorFootprints(bamfile, pfm=CTCF[[1]], + genome=Hsapiens,bindingSites=which, + min.score="95%", seqlev="chr1", + upstream=100, downstream=100) Error in factorFootprints(bamfile, pfm = CTCF[[1]], genome = Hsapiens, : length(bindingSites$score) == length(bindingSites) is not TRUE Called from: factorFootprints(bamfile, pfm = CTCF[[1]], genome = Hsapiens, bindingSites = which, min.score = "95%", seqlev = "chr1", upstream = 100, downstream = 100) Browse[1]> Q > bamfile [1] "/home/aiminyan/R/x86_64-pc-linux-gnu-library/3.5/ATACseqQC/extdata/GL1.bam" > which <- GRanges(seqnames = "chr1", ranges = IRanges(seq(10000, 20000, by=1000), width = 10)) > seqinfo(which) <- seqinfo(Hsapiens) > factorFootprints(bamfile, pfm=CTCF[[1]], + genome=Hsapiens,bindingSites=which, + min.score="95%", seqlev="chr1", + upstream=100, downstream=100) Error in factorFootprints(bamfile, pfm = CTCF[[1]], genome = Hsapiens, : length(bindingSites$score) == length(bindingSites) is not TRUE Called from: factorFootprints(bamfile, pfm = CTCF[[1]], genome = Hsapiens, bindingSites = which, min.score = "95%", seqlev = "chr1", upstream = 100, downstream = 100) Browse[1]> On Thu, Jul 12, 2018 at 3:36 PM, Ou, Jianhong [bioc] < noreply@bioconductor.org> wrote: > Activity on a post you are following on support.bioconductor.org > > User Ou, Jianhong <https: support.bioconductor.org="" u="" 4539=""/> wrote Answer: > use factorFootprints function in ATACseqQC to get DNA footprinting for > certain gene region <https: support.bioconductor.org="" p="" 110957="" #110960="">: > > Hi Aimin, > > Thank you for trying ATACseqQC to plot the TF footprints. You can have a > try following code: > > which <- GRanges(seqnames = "chr1", ranges = IRanges(seq(10000, 20000, by=1000), width = 10)) > seqinfo(which) <- seqinfo(Hsapiens) > factorFootprints(bamfile, pfm=CTCF[[1]], > genome=Hsapiens,bindingSites=which, > min.score="95%", seqlev="chr1", > upstream=100, downstream=100) > > Let me know if you have any trouble in running this sample. > > Jianhong. > > ------------------------------ > > Post tags: ATACseqQC > > You may reply via email or visit https://support.bioconductor. > org/p/110957/#110960 >
ADD REPLY •
0
Entering edit mode

try to add scores for the binding sites like this 

which$score <- rep(1, length(which))

ADD REPLY •
0
Entering edit mode
Yes, I tried, but I got the following error: bamfile <- system.file("extdata", "GL1.bam",package="ATACseqQC") which <- GRanges(seqnames = "chr1", ranges = IRanges(seq(10000, 20000, by=1000), width = 10)) seqinfo(which) <- seqinfo(Hsapiens) which$score <- rep(1, length(which)) factorFootprints(bamfile, pfm=CTCF[[1]],genome=Hsapiens,bindingSites=which,min.score="95%", seqlev="chr1",upstream=100, downstream=100) Error in featureAlignedSignal(cvglists = cvglist, feature.gr = reCenterPeaks(mt, : Length of feature.gr less than 2. In addition: Warning message: In featureAlignedSignal(cvglists = cvglist, feature.gr = reCenterPeaks(mt, : feature.gr is set to the center of feature.gr On Thu, Jul 12, 2018 at 5:12 PM, Ou, Jianhong [bioc] < noreply@bioconductor.org> wrote: > Activity on a post you are following on support.bioconductor.org > > User Ou, Jianhong <https: support.bioconductor.org="" u="" 4539=""/> wrote Comment: > use factorFootprints function in ATACseqQC to get DNA footprinting for > certain gene region <https: support.bioconductor.org="" p="" 110957="" #110962="">: > > try to add scores for the binding sites like this > > which$score <- rep(1, length(which)) > > ------------------------------ > > Post tags: ATACseqQC > > You may reply via email or visit https://support.bioconductor. > org/p/110957/#110962 >
ADD REPLY •
0
Entering edit mode

Aimin,

 

Please use my following example, 

 

which <- GRanges(seqnames = "chr1", ranges = IRanges(seq(10000, 2000000, by=1000), width = 100), score = rep(1,1991))

seqlengths(which)[1] = 248956422

Best,

 

Julie

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

Aimin,

Score is required here. The following example should run fine.

which <- GRanges(seqnames = "chr1", ranges = IRanges(seq(10000, 2000000, by=1000), width = 100), score = rep(1,1991))

seqlengths(which)[1] = 248956422

bamfile <- system.file("extdata", "GL1.bam",

                            package="ATACseqQC")

     library(MotifDb)

     CTCF <- query(MotifDb, c("CTCF"))

     CTCF <- as.list(CTCF)

     library(BSgenome.Hsapiens.UCSC.hg19)

factorFootprints(bamfile, pfm=CTCF[[1]],genome=Hsapiens, min.score="95%", seqlev="chr1",upstream=100, downstream=100, bindingSites=which)

I suggest not set bindingSites unless you know the binding sites.

Best,

 

Julie
ADD COMMENT •
0
Entering edit mode
OK, thank you very much. it works now. Aimin On Thu, Jul 12, 2018 at 5:23 PM, Julie Zhu [bioc] <noreply@bioconductor.org> wrote: > Activity on a post you are following on support.bioconductor.org > > User Julie Zhu <https: support.bioconductor.org="" u="" 3596=""/> wrote Answer: > use factorFootprints function in ATACseqQC to get DNA footprinting for > certain gene region <https: support.bioconductor.org="" p="" 110957="" #110964="">: > > Aimin, > > *Score is required here. The following example should run fine.* > > which <- GRanges(seqnames = "chr1", ranges = IRanges(seq(10000, 2000000, > by=1000), width = 100), score = rep(1,1991)) > > seqlengths(which)[1] = 248956422 > > bamfile <- system.file("extdata", "GL1.bam", > > package="ATACseqQC") > > library(MotifDb) > > CTCF <- query(MotifDb, c("CTCF")) > > CTCF <- as.list(CTCF) > > library(BSgenome.Hsapiens.UCSC.hg19) > > factorFootprints(bamfile, pfm=CTCF[[1]],genome=Hsapiens, min.score="95%", > seqlev="chr1",upstream=100, downstream=100, bindingSites=which) > > *I suggest not set bindingSites unless you know the binding sites.* > > Best, > > > > Julie > > ------------------------------ > > Post tags: ATACseqQC > > You may reply via email or visit https://support.bioconductor. > org/p/110957/#110964 >
ADD REPLY •
0
Entering edit mode
Aimin, I suggest not set bindingSites unless you really know where the binding sites are. Best, Julie
ADD REPLY •

Login before adding your answer.

Similar Posts
Loading Similar Posts
Traffic: 753 users visited in the last hour
Content Search
Users
Tags
Badges
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