Seqanswers Leaderboard Ad

Collapse

Announcement

Collapse
No announcement yet.
X
 
  • Filter
  • Time
  • Show
Clear All
new posts

  • bedGraphToBigWig does not work with single base coverage

    Hi,

    I've identified a "bug" in UCSC's bedGraphToBigWig.
    It does not work with single base coverage.

    If I have a file where the start and end are always equal (i.e. single base coverage), the bigWig file generated by bedGraphToBigWig is empty.
    Does anyone have a workaround, or am I simply not respecting the format for bedgraph files? If I want to include the count just for base 1, how am I supposed to write it?
    1 1 count
    or
    1 2 count
    The 1st option appears correct but bedGraphToBigWig will not take it.
    It would be a big job to convert all my bedgraph files to the format bedGraphToBigWig appears to want.

    Here are 2 example bedgraph files: test1.bedgraph with the start and the end equal to 1, and test2.bedgraph where the start and end are 1 and 2, respectively. The bigWig files generated by bedGraphToBigWig is empty in the first case, and with 1 base covered in the second case.

    blancha@sw-2r12-n36 ~]$ more test1.bedgraph
    chr1 1 1 1
    [blancha@sw-2r12-n36 ~]$ more test2.bedgraph
    chr1 1 2 1
    [blancha@sw-2r12-n36 ~]$ bedGraphToBigWig test1.bedgraph UCSCChromInfo.txt test1.bw
    [blancha@sw-2r12-n36 ~]$ bedGraphToBigWig test2.bedgraph UCSCChromInfo.txt test2.bw
    [blancha@sw-2r12-n36 ~]$ bigWigInfo test1.bw
    version: 4
    isCompressed: yes
    isSwapped: 0
    primaryDataSize: 30
    primaryIndexSize: 6,204
    zoomLevels: 1
    chromCount: 1
    basesCovered: 0
    mean: -nan
    min: 1.000000
    max: 1.000000
    std: -nan
    [blancha@sw-2r12-n36 ~]$ bigWigInfo test2.bw
    version: 4
    isCompressed: yes
    isSwapped: 0
    primaryDataSize: 32
    primaryIndexSize: 6,204
    zoomLevels: 1
    chromCount: 1
    basesCovered: 1
    mean: 1.000000
    min: 1.000000
    max: 1.000000
    std: 0.000000

    Thank you for your help.
    Last edited by blancha; 05-28-2014, 01:51 PM.

  • #2
    I think your bedgraph should be the second way, 1 2 not 1 1, because bedgraph files are "half open".

    Comment


    • #3
      Great. Thanks BioComputer.
      That is exactly the formal definition I was looking for.

      I did waste half a day figuring out the "bug", though.
      Thankfully, you had a quick answer to my question, as I would never have figured out the reason behind the "bug" on my own.
      I guess it's more of a "feature".
      Now, I just need to decide if I should write a Python script to correct the coordinates, or whether if I run bismark_methylation_extractor again with the --zero_based parameter I will get the desired format.

      I used bismark_methylation_extractor to generate the BED files. The output was a bit confusing since, by default, the bedgraph file outputted has 0-based start coords and 1-based end coords, whereas the other BED files generated have 1-based genomic coords.

      To be fair, it is in the documentation, although I'm still unsure if adding the --zero-based parameter will result in 0-based start and end coordinates, or just 0-based start coordinates, and 1-based end coordinates. I imagine it is the latter.
      Bioinformatics is so much fun.

      Thanks again!

      Comment


      • #4
        So you just have to add 1 to each value in the end position column? I'm still learning this stuff so there's probably a better way to do this, but this command below worked when I tested it on a small file. It assumes tab-separated columns, should add 1 to each value in the 3rd column, and should print out any non-position lines (ie. any lines that don't start with "chr" like a file header) as is. Test first and use at your own risk

        Code:
        awk '/^chr/ {FS="\t";$3=$3+1;print} !/^chr/{print}' inputFile > outputFile
        Or instead of adding 1 to the end position do you need to substract 1 from each start position?

        Code:
        awk '/^chr/ {FS="\t";$2=$2-1;print} !/^chr/{print}' inputFile > outputFile
        Here's another version that should do the same thing, maybe it'll be faster?

        Code:
        awk 'BEGIN { OFS="\t"}{if (/^chr/) $3=$3+1; print}' inputFile > outputFile
        Last edited by biocomputer; 05-28-2014, 05:35 PM.

        Comment


        • #5
          If you are going to use UCSC browser to display the files consider this: http://genome.ucsc.edu/FAQ/FAQtracks.html#tracks1

          Comment

          Latest Articles

          Collapse

          • seqadmin
            Recent Developments in Metagenomics
            by seqadmin





            Metagenomics has improved the way researchers study microorganisms across diverse environments. Historically, studying microorganisms relied on culturing them in the lab, a method that limits the investigation of many species since most are unculturable1. Metagenomics overcomes these issues by allowing the study of microorganisms regardless of their ability to be cultured or the environments they inhabit. Over time, the field has evolved, especially with the advent...
            09-23-2024, 06:35 AM
          • seqadmin
            Understanding Genetic Influence on Infectious Disease
            by seqadmin




            During the COVID-19 pandemic, scientists observed that while some individuals experienced severe illness when infected with SARS-CoV-2, others were barely affected. These disparities left researchers and clinicians wondering what causes the wide variations in response to viral infections and what role genetics plays.

            Jean-Laurent Casanova, M.D., Ph.D., Professor at Rockefeller University, is a leading expert in this crossover between genetics and infectious...
            09-09-2024, 10:59 AM

          ad_right_rmr

          Collapse

          News

          Collapse

          Topics Statistics Last Post
          Started by seqadmin, 10-02-2024, 04:51 AM
          0 responses
          13 views
          0 likes
          Last Post seqadmin  
          Started by seqadmin, 10-01-2024, 07:10 AM
          0 responses
          21 views
          0 likes
          Last Post seqadmin  
          Started by seqadmin, 09-30-2024, 08:33 AM
          0 responses
          25 views
          0 likes
          Last Post seqadmin  
          Started by seqadmin, 09-26-2024, 12:57 PM
          0 responses
          18 views
          0 likes
          Last Post seqadmin  
          Working...
          X