Hello,
I am trying to use a for loop in python to loop through a bam file and take out reads that map to certain chromosomes, for further downstream analysis. I am having problems in the beginning, seemingly easiest part of my code.
samfile = pysam.Samfile("mybam.bam", "rb")
positions=["7:151970856", "7:151970856"]
for p in positions:
start = p.split(":")
chr=start[0]
pos=start[1]
chromreads=[]
for read in samfile:
if int(chr) == int(read.tid):
chromreads.append(read)
print(len(chromreads))
Output I'm getting:
2525527
0
Output I expect:
2525527
2525527
It seems like it's not even looping through the second position. What am I doing wrong?
Thank you.
2 answers
Answered my own question from StackOverflow
You're iterating through the entire Samfile object, consuming it on the first pass of your inner for loop. When you arrive at the second iteration of your outer for loop, there is nothing left in the Samfile iterator and so it raises StopIteration which causes the inner loop to exit.
•
written
by
Matt Shirley
Log in to answer this question.
-
In pysam pileup, discrepancy between stepper all and stepper samtools reads
written by Michael •
-
Convert hg38 to a nested array for binary search
written by noahhelton98 •
-
Looping over list in R
written by rodolfo.peacewalker •
-
Use AWK to interrogate one file through another rather than for and while loops
written by s1667153 •
-
reads name for alt base pysam
written by novicebioinforesearcher •
-
Using pysam to parse bamfile and fetch based on bedfile, looking to speed up code.
written by Daniel
-
Problem with a "chr" variable in plotManhattan (Sushi R package)
written by yotiao •
-
obtaining a pileup with read names from a bam alignment
written by bernardo1963 •
-
How to generate vcf using samtools Version: 0.1.19-96b5f2294a?
written by lakhujanivijay
-
Program to split bam file according to size
written by kavitarege •
