On Tue, Mar 23, 2010 at 5:24 PM, Chen-Yi Chen <[email protected]> wrote:

>
>
> ----- Original Message -----
> From: Michael Lawrence <[email protected]>
> Date: Tuesday, March 23, 2010 4:57 pm
> Subject: Re: [Bioc-sig-seq] Questions about coverage visualization in wig
> file output
> To: Chen-Yi Chen <[email protected]>
> Cc: [email protected]
>
> > On Tue, Mar 23, 2010 at 4:46 PM, Chen-Yi Chen <[email protected]>
> > wrote:
> > > Hi all,
> > > Recently we found the greatness of this bioconductor ht-seq
> > pacakage and
> > > we would like to do some analysis in this package. We have
> > aligned our
> > > ChIP and input data by Bowtie aligner and try to output wig file
> > for the
> > > coverage, so we are able to use UCSC genome browser to visualize
> > them. By
> > > searching the documentation and work flow online, I found this great
> > > thread in bioc-sig-sequencing archive:
> > >
> > >
> > > https://stat.ethz.ch/pipermail/bioc-sig-sequencing/2009-
> > September/000632.html>
> > > We basically followed the work flow in the above thread, with some
> > > slightly modifications:
> > > ChIPread = readAligned("./", pattern="0_Rad_H3K4_ChIP.out",
> > type="Bowtie")> ChIP = as(ChIPread, "GenomeData")
> > > inputread = readAligned('./', pattern='0_Rad_input.out',
> > type='Bowtie')> input = as(inputread, "GenomeData")
> > > gdlist = GenomeDataList(list(input=input, ChIP=ChIP))
> > > input1 = extendReads(gdlist$input$chr1, seqLen=200)
> > > ChIP1 = extendReads(gdlist$ChIP$chr1, seqLen=200)
> > > library(BSgenome.Hsapiens.UCSC.hg18)
> > > human.chromlens = seqlengths(Hsapiens)
> > > input1.cov = coverage(input1, width=human.chromlens['chr1'])
> > > ChIP1.cov = coverage(ChIP1, width=human.chromlens['chr1'])
> > > inputrange1 = IRanges(start=start(input1.cov), end=end(input1.cov))
> > > score = runValue(input1.cov)
> > > chr1 = rep("chr1", length(inputrange1))#
> > > export(GenomicData(inputrange1, score, chrom=chr1),
> > "input_test.wig")>
> > > We successfully got a wig file output, however, when we try to
> > visualize> it on UCSC genome browser, the peaks won't display
> > correctly.>
> >
> > Could you please give a little more detail? A screenshot perhaps?
> I attached a UCSC genome browser screen shot.
>

Thanks, but it's not clear what's wrong there. You're looking at an entire
chromosome in heatmap mode. Could you maybe describe in words what's wrong?



> >
> >
> > > I took a glance on the output wig file, it has the format like the
> > > following:
> > > track name="R Track" type=wiggle_0
> > > variableStep chrom=chr1 span=1
> > > 251     7
> > > 395     7
> > > 14041   3
> > > 78204   1
> > > ...
> > > variableStep chrom=chr1 span=2
> > > 12      2
> > > 212     9
> > > 346     7
> > > ...
> > > variableStep chrom=chr1 span=3
> > > 343     8
> > > 404     7
> > > ...
> > > ...
> > > fixedStep chrom=chr1 start=117727703 step=129030611 span=1990
> > > 0
> > > 0
> > > variableStep chrom=chr1 span=1991
> > > ...
> > >
> > > Apparently the output wig file format we got is different than
> > what's> mentioned in the thread, and I believed those different
> > variable steps and
> > > fixed steps are confusing the UCSC genome browser.
> > >
> >
> > I would hope not. It looks like valid WIG. It's compressing the
> > data using
> > multiple variableStep and fixedStep blocks.
> >
> I believed so too, but I've also emailed the supports in UCSC genome
> browser. They said that usually WIG file only has very few number of
> "variable steps", and the fixedSteps in our output WIG file didn't make
> sense.
>
>
This is probably because people rarely try to encode coverage in WIG, given
bedGraph. Could you be a bit more specific about how the fixedSteps don't
make sense?

Thanks,
Michael


> >
> > > Is there a way to output the wig file in the format without all
> > those> steps, for example,
> > > track name="R Track" type=wiggle_0
> > > chrI    0       2       169
> > > chrI    2       4       176
> > > chrI    4       5       178
> > > chrI    5       7       179
> > > ...
> > > which mentioned in the thread?
> > >
> > >
> > This is no longer valid WIG. It's now called bedGraph. So you can
> > export(track, "input_test.bedGraph").
> Thanks Michael!
> The coverage graph displays pretty well in bedGraph format!
>
> Thanks,
> -Charlie-
>
> >
> > But I'd like to figure out the problem with the WIG, if possible.
> >
> > Thanks,
> > Michael
> >
> >
> >
> >
> > > Thanks for your help!
> > >
> > > -Charlie-
> > >
> > > _______________________________________________
> > > Bioc-sig-sequencing mailing list
> > > [email protected]
> > > https://stat.ethz.ch/mailman/listinfo/bioc-sig-sequencing
> > >
> >
>

        [[alternative HTML version deleted]]

_______________________________________________
Bioc-sig-sequencing mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/bioc-sig-sequencing

Reply via email to