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
            Current Approaches to Protein Sequencing
            by seqadmin


            Proteins are often described as the workhorses of the cell, and identifying their sequences is key to understanding their role in biological processes and disease. Currently, the most common technique used to determine protein sequences is mass spectrometry. While still a valuable tool, mass spectrometry faces several limitations and requires a highly experienced scientist familiar with the equipment to operate it. Additionally, other proteomic methods, like affinity assays, are constrained...
            04-04-2024, 04:25 PM
          • seqadmin
            Strategies for Sequencing Challenging Samples
            by seqadmin


            Despite advancements in sequencing platforms and related sample preparation technologies, certain sample types continue to present significant challenges that can compromise sequencing results. Pedro Echave, Senior Manager of the Global Business Segment at Revvity, explained that the success of a sequencing experiment ultimately depends on the amount and integrity of the nucleic acid template (RNA or DNA) obtained from a sample. “The better the quality of the nucleic acid isolated...
            03-22-2024, 06:39 AM

          ad_right_rmr

          Collapse

          News

          Collapse

          Topics Statistics Last Post
          Started by seqadmin, 04-11-2024, 12:08 PM
          0 responses
          30 views
          0 likes
          Last Post seqadmin  
          Started by seqadmin, 04-10-2024, 10:19 PM
          0 responses
          32 views
          0 likes
          Last Post seqadmin  
          Started by seqadmin, 04-10-2024, 09:21 AM
          0 responses
          28 views
          0 likes
          Last Post seqadmin  
          Started by seqadmin, 04-04-2024, 09:00 AM
          0 responses
          52 views
          0 likes
          Last Post seqadmin  
          Working...
          X