[howTo] Access GenomeArray() using a GenomeInterval() with a negative start

[shouldsee]

I was trying to iadd a GenomicArray() with a GenomicInterval() at negative start, and got the following error. How much effort would it take to make GenomicArray() support negative index?

arr = HTSeq.GenomicArray(chroms='auto')
iv = HTSeq.GenomicInterval('chr1',-10,50,'+')
arr[iv] += 1

IndexErrorTraceback (most recent call last)
<ipython-input-257-5253053be354> in <module>()
      1 arr = HTSeq.GenomicArray(chroms='auto')
      2 iv = HTSeq.GenomicInterval('chr1',-10,50,'+')
----> 3 arr[iv] += 1

python2/src/HTSeq/_HTSeq.pyx in HTSeq._HTSeq.GenomicArray.__getitem__()

python2/src/HTSeq/_HTSeq.pyx in HTSeq._HTSeq.ChromVector.__getitem__()

IndexError: start too small

[shouldsee]

NVM. I found a solution using HTSeq.StepVector.StepVector() instead. It may be benefical to subclass GenomicArray() to allow negative indexing.

iv = HTSeq.GenomicInterval('chr1',-5,30,'+')
vct = HTSeq.StepVector.StepVector.create(length=60,start_index=-10)
vct[iv.start:iv.end] +=1
list(vct)

[shouldsee]

Is HTSeq under active development at the moment?

[iosonofabio]

Yes, but we're not paid so PRs are welcome 😊

On Fri, May 3, 2019, at 03:05, shouldsee wrote:

Is HTSeq under active development at the moment?

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/simon-anders/htseq/issues/81#issuecomment-489043370, or mute the thread https://github.com/notifications/unsubscribe-auth/AAJFEABKLQVA3AUFGFH3DA3PTQE5FANCNFSM4HJ6A4QA.

[iosonofabio]

@shouldsee It seems like you found a workaround, would you like to try and see how to add negative indexing to GenomicArray or should I close this?

[shouldsee]

It was only for plotting purposes and the StepVector worked fine so I think we can just add this Tip somewhere: Use StepVector() if you are dealing with negative index.

I ended up plotting segments without a unified iterator so I don't really need this functionality anymore. But always good to save the snippet as an example.

[iosonofabio]

sounds good. where in the docs would you add it?

[shouldsee]

How about https://github.com/simon-anders/htseq/blob/master/python2/doc/genomic.rst

[shouldsee]

I will be PR in 5 mins

[iosonofabio]

Awesome thanks

[iosonofabio]

This got implemented meanwhile - I mentioned you in the history changes. Closing