coverage histogram causes malformedBedLineError

[radaniba]

Hello

I am trying to create a histogram using pybedttols and coverage wrapper but when I look at the result (the header) it looks fine except the option hist cause the MalformedBedLineError since it is a not a valid bed file

here is the error message

>>> for position in coverage:
...     if(position[0]=="all"):
...             print position
...
...
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "cbedtools.pyx", line 680, in pybedtools.cbedtools.IntervalFile.__next__ (pybedtools/cbedtools.cpp:8685)
pybedtools.cbedtools.MalformedBedLineError: malformed line: ['all', '0', '981', '30513', '0.0321502']

my command to get the coverage was :

coverage = alignment.coverage(regions, hist=True)

Any idea on how to get only the all lines ?

Thanks

[daler]

Yep, this is happening because the all lines have one less field than the rest of the file. For performance reasons, BedTools (and pybedtools) check the first valid BED line and then assume that all lines in the file have the same number of fields.

Running any BedTools operations on that coverage file with the "all" histogram lines at the end will make BedTools complain, which is then passed to pybedtools to raise an exception.

On the command line, you could just grep out the all lines. In pybedtools, you could just use typical file parsing stuff.

Here's a solution that takes the results from a BedTool.coverage(regions, hist=True) call and returns two values: 1) a BedTool object representing the non-"all", valid BED lines and 2) a separate pandas.DataFrame of the all lines.

https://gist.github.com/daler/99cdf714d438e161dad6

[radaniba]

Thank you @daler !

I will try your solution and try to capture the all lines into a dataframe it is indeed a good solution

Thanks
Rad

[adomingues]

Hi @daler
I've just hit up a little snag with this issue, that is the output of coverageBed not being a bed - which is fine. To cut a long story short, I wanted to convert this output into a bed file to use in other operations (closestBed). So I was calculating coverage:

import pybedtools
in_reads = pybedtools.example_bedtool('x.bam')
strand = '+'
three_bg = in_reads.genome_coverage(dz=True, strand=strand, **{'3': True})

And created a function, cov2bed6 that with the help of extend_fields and show_value would change this:

chr2L 4999310 1

to this:

chr2L 4999310 4999311 . 1 +

However, since three_bg is not a proper bed it whole thing just crashed and burn. Do you happen to know of a handy solution to convert the coverage output into bed? Otherwise I will just wrangle something with Pandas.

Out of interest, the cov2bed6 looks like this:

def cov2bed6(bg, strand):
    """
    Converts single nucleotide coverage (bedTool):
    chr2L   4999307 1
    chr2L   4999308 1
    chr2L   4999310 1

    to bed6 (bedTool) in which the score is put in 5th field. The name is set to ".". If no strand is provided - not recommend! - "." is the default.
    Returns a bedTool obj
    """
    assert strand is "+" or strand is "-" or strand is "."
    bead4 = bg.each(add_end).saveas()
    stranded = bead4.each(add_strand, strand).saveas()
    bed6 = stranded.each(update_score).saveas()
    return bed6

Cheers.

[daler]

Yep, in this case I would switch to pandas for the manipulation and conversion into a BED-like format.

[adomingues]

Just in case someone needs a similar solution, here is my function to convert coverage (chr position score) into bed4:

def coverage_to_bedTool(cov):
    """
    Converts single nucleotide coverage (bedTool):
    chr2L   4999307 1
    chr2L   4999308 1
    chr2L   4999310 1

    to bed4 (bedTool). This step is needed because bedTools expects col2 < col3 which is not always the case since col3 is the coverage. The end coordinate is position (col2) + 1. See:
    https://github.com/daler/pybedtools/issues/145
    https://github.com/daler/pybedtools/issues/110
    https://github.com/daler/pybedtools/issues/112

    Returns a bedTool obj
    Requires pandas.
    """
    df = pd.read_table(cov.fn, names=['chrom', 'start', 'coverage'])
    df['end'] = df['start'] + 1
    df = df[['chrom', 'start', 'end', 'coverage']]
    bed4 = pybedtools.BedTool.from_dataframe(df)
    return bed4

Btw, and out of curiosity, even though cov.head() returns a "malformedBedLineError" message, cov.tail() works fine. Just use my previous example and see the difference between three_bg.tail() and three_bg.head().

[daler]

Thanks.

The different head() and tail() behavior is from head doing the easy thing -- "for i in bedtool; print(i)" -- while tail breaks the last bytes of the file off into chunks to find the newlines to avoid having to iterate over the entire file.

After thinking about it, it probably makes sense to have head work on the file directly just like tail does so that debugging the MalformedBedLineError messages would be easier. I'll open an issue for that.

And seeing the issues referenced in your docstring reminds me that it hasn't been just one person hitting the assumption that coverage output is a BED file. So maybe it's worth adding something like you've shown as a helper function. Does this make sense as a method on a BedTool object? Or a function in the helpers module?

[adomingues]

Does this make sense as a method on a BedTool object? Or a function in the helpers module?

I would probably make it a helpers function because it is very case-specific for some of the coverage function output, and it could keep the class as is. If anyone needs to change format, it would simply call the function and all would be good. But you know more about the specific of the BedTool objects.

[daler]

It might be convenient to include in a "pipeline" like:

x = (
    BedTool('a.bam')
    .genome_coverage(dz=True)
    .coverage_to_bed()
    .intersect('b.bed')
)

In which case being a method would be the way to go. I was just checking out helpers.py and it doesn't make sense with the other kinds of functions there.

[adomingues]

You are probably right. I was thinking more along lines of functions like pybedtools.featurefuncs.three_prime.

Thank you once again for being so damn fast fixing things.