Pysam fetch for records not exists in VCF files
2
0
Entering edit mode
7.5 years ago
always_learning ★ 1.1k
f1 = pysam.VariantFile("test.vcf.gz")

for f in f1.fetch(str(chr), int(start), int(end)):

            if f:
                    print "yes"  
            else :
                    print "hello"

In the above example Its not printing "hello" even for cases which doesn't have entry into VCF. I want to print program "hello" for cases where VariantFile doesn't have variants records.

Its their a way to do this ?

Thanks

pysam fetch • 4.0k views
ADD COMMENT
0
Entering edit mode

Try:

if not f1.fetch(str(chr), int(start), int(end)):
    print("hello")

instead of the for loop.

ADD REPLY
0
Entering edit mode

Thanks Ryan, Its not working. I tried this earlier as well.

ADD REPLY
0
Entering edit mode

What do you get if you print f for positions which are not in the vcf?

ADD REPLY
0
Entering edit mode

Its not going inside else condition.

ADD REPLY
1
Entering edit mode

There shouldn't be an else condition, since there shouldn't be a for loop in this case in question.

ADD REPLY
0
Entering edit mode

And looking at how f looks like doesn't depend on if else statement, you just need to now what fetch returns for a position you think isn't included in the vcf.

ADD REPLY
0
Entering edit mode
     if f1.fetch(str(chr), int(start), int(end)):
                print "hello"
     else:
                print "no"

Even this is not working as well.

ADD REPLY
0
Entering edit mode

What do you get if you print f1.fetch(str(chr), int(start), int(end)) for positions which are not in the vcf?

ADD REPLY
0
Entering edit mode
>>> print f1.fetch('2', 224629875, 224629876)
<pysam.ctabix.TabixIterator object at 0x7fd7e26101f0>
>>> print f1.fetch('2', 224629873, 224629874)
<pysam.ctabix.TabixIterator object at 0x7fd7e26101f0>

Reagrdless of absence or present its proting iterator object. I checked for object condition as well but no luck.

ADD REPLY
1
Entering edit mode

The TabixIterator is always True, regardless of precense or absence of that variant.

ADD REPLY
0
Entering edit mode

yes and Thats why its never going to else condition.

ADD REPLY
0
Entering edit mode

Exactly. What are you trying to achieve? This issue probably has a better solution.

ADD REPLY
0
Entering edit mode

you means other ways Using Pysam ?

ADD REPLY
0
Entering edit mode
7.5 years ago
always_learning ★ 1.1k

SInce I have a small files need to extracted from VCF file so i jus used hash and then search in hash so finally I was able to accomplished this task but stiil i have no idea why its not printing "hello" even for cases which doesn't have entry into VCF

ADD COMMENT
0
Entering edit mode
6.0 years ago
Vasisht ▴ 190

Pysam returns an iterator for the fetch call, so that condition will be always be true. You can query for an attribute of the VariantRecord (e.g. contig, pos etc.) to check if it's empty.

for f in f1.fetch(str(chr), int(start), int(end)):
      # This will return an empty list if a variant is not present at that location.
      if [r.contig for r in f]:
          print True
      else:
          print False
ADD COMMENT
0
Entering edit mode

I'd like to suggest to write python3-compatible code, especially when you are sharing examples. In this case you have to adapt the print statements to print(True) and print(False). As such, users who stumble upon your solution won't have issues running your code.

ADD REPLY

Login before adding your answer.

Traffic: 1963 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6