Skip to content

Python: when c2 <= c1 output is not flipped back #119

@lldelisle

Description

@lldelisle

Describe the bug
When fetching the contacts between 'chr1' 'chrM', it gives the same result as 'chrM' 'chr1' without any notice to the user.

To Reproduce

import numpy as np
import hicstraw
hic_file =  'ENCFF080DPJ.hic'
chrom1 = 'chr1'
chrom2 = 'chr2'
result = hicstraw.straw('observed', 'NONE', hic_file, chrom1, chrom2, 'BP', 1000000)
for i in range(10):
     print("{0}\t{1}\t{2}".format(result[i].binX, result[i].binY, result[i].counts))
result = hicstraw.straw('observed', 'NONE', hic_file, chrom2, chrom1, 'BP', 1000000)
for i in range(10):
     print("{0}\t{1}\t{2}".format(result[i].binX, result[i].binY, result[i].counts))

Expected behavior
I would expect to have the same result but with column 1 and 2 shifted.
I got exactly the same result.

The problem is that the chromosomes are flipped here:

} else { // flip
this->c1 = c02;
this->c2 = c01;
this->numBins1 = static_cast<int32_t>(chrom2.length / resolution);
this->numBins2 = static_cast<int32_t>(chrom1.length / resolution);
}
but this is information is not stored.
Therefore, the results are not flipped back in:
if (!isnan(c) && !isinf(c)){
contactRecord record = contactRecord();
record.binX = static_cast<int32_t>(x);
record.binY = static_cast<int32_t>(y);
record.counts = c;
records.push_back(record);
}

and the matrix is not transposed in:
return py::array(py::cast(finalMatrix));

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions