key not found in last pysam.PileupColumn
- 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.