JuliaPy / JuliaPy/PyCall.jl

key not found in last pysam.PileupColumn

Open
#168 24 comments 0 reactions 0 assignees View on GitHub
Dominant language
Julia
Stars
1.5k
Forks
186
PR merge metrics
No merged PRs in 30d

Description

Dear PyCall developper,

I'm new at programming in julia and I'm trying to use the pysam python library using PyCall.
I work with the following test file describing the alignment of two sequencing reads on a reference:

`````` bash
bli@sei-lot:/tmp$ cat test.sam
@HD VN:1.0 SO:unsorted
@SQ SN:test_ref LN:17637
SRR1524970.144283 16 test_ref 1706 255 25M * 0 0 TGCTGATGAAGCAGAACAACTTTAA ]YG[^baaaa^W`ab]]````aaba AS:i:0XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:25 YT:Z:UU
SRR1524970.316478 16 test_ref 1706 255 24M * 0 0 TGCTGATGAAGCAGAACAACTTTA `\X_`aaaaaY]``b_aa_aaaaa AS:i:0XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:24 YT:Z:UU
``````

This file is then sorted and indexed using [samtools](http://www.htslib.org/):

``` bash
bli@sei-lot:/tmp$ samtools sort test.sam test_sorted
bli@sei-lot:/tmp$ samtools index test_sorted.bam
bli@sei-lot:/tmp$ ls test*
test.jl test.py test.sam test.sam.old test_sorted.bam test_sorted.bam.bai
```

And I want to count reads overlapping each reference position.

Here is how I do it in python:

``` bash
bli@sei-lot:/tmp$ cat test.py
```

``` python
#!/usr/bin/env python

import sys
import pysam

def main():
samfile = pysam.Samfile(sys.argv[1], "rb")
for (ref, reflen) in zip(samfile.references, samfile.lengths):
print "%s: %d" % (ref, reflen)
for pile in samfile.pileup(ref):
ref_pos = pile.pos
nb_reads = 0
print ref_pos,
for pileup_read in pile.pileups:
nb_reads += 1
print nb_reads,
print

main()
```

``` bash
bli@sei-lot:/tmp$ python test.py test_sorted.bam
test_ref: 17637
1705 1 2
1706 1 2
1707 1 2
1708 1 2
1709 1 2
1710 1 2
1711 1 2
1712 1 2
1713 1 2
1714 1 2
1715 1 2
1716 1 2
1717 1 2
1718 1 2
1719 1 2
1720 1 2
1721 1 2
1722 1 2
1723 1 2
1724 1 2
1725 1 2
1726 1 2
1727 1 2
1728 1 2
1729 1
```

And here is how I tried to translate the code in julia:

``` bash
bli@sei-lot:/tmp$ cat test.jl
```

``` julia
#!/usr/bin/env julia

using PyCall
@pyimport pysam

function main()
samfile = pysam.Samfile(ARGS[1], "rb")
for (ref, reflen) in zip(samfile["references"], samfile["lengths"])
println("$ref: $reflen")
for pile in pycall(samfile["pileup"], PyAny, ref)
ref_pos = convert(PyAny, pile["pos"])
nb_reads = 0
print("$ref_pos")
for pileup_read in pile["pileups"]
nb_reads += 1
print(" $nb_reads")
end
println()
end
end
end

main()
```

``` bash
bli@sei-lot:/tmp$ julia test.jl test_sorted.bam
test_ref: 17637
1705 1 2
1706 1 2
1707 1 2
1708 1 2
1709 1 2
1710 1 2
1711 1 2
1712 1 2
1713 1 2
1714 1 2
1715 1 2
1716 1 2
1717 1 2
1718 1 2
1719 1 2
1720 1 2
1721 1 2
1722 1 2
1723 1 2
1724 1 2
1725 1 2
1726 1 2
1727 1 2
1728 1 2
1729ERROR: key not found: "pileups"
in getindex at /home/bli/.julia/v0.3/PyCall/src/PyCall.jl:255
in main at /tmp/test.jl:14
in include at ./boot.jl:245
in include_from_node1 at loading.jl:128
in process_options at ./client.jl:285
in _start at ./client.jl:354
while loading /tmp/test.jl, in expression starting on line 23
```

It seems that the last [PileupColumn](http://pysam.readthedocs.org/en/latest/api.html#pysam.PileupColumn) is" anormal" when accessed from within julia.

Do you have a clue of what may be causing this problem ?

I take the opportunity to ask another question: is there a difference between the `object[:attribute]` and `object["attribute"]` syntaxes?

Contributor guide

No contributing guide indexed for this repository

Assessment

This issue has not been assessed yet.

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.