Extracting insert size from .bam using pysam

[Tom.Booker]

Hi there, this is my first time posting. I am fairly new to bioinformatics so I may make silly errors here and there.

Anyway, I have been trying to extract the mate-pair distances from an indexed bam file of paired end reads. I have been using pysam to iterate through a a bam using the pysam.Samfile().fetch() methods. By manually inspecting my output I know that I have the reads that I am interested in. The problem is that pretty much all the reads extracted using pysam have in incorrect (or possibly missing!) value in the insert size (isize) field which I have confirmed by extracting a few corresponding reads from the sam file. All fields match up but for the overwhelming majority of reads the insert size field have an exact value of 101, which made me think that there could be another value taking the place of isize in the pysam.AlignedRead object.

Does anyone with experience using pysam know what might be going on?

Here is the code I've been using, if it helps:

sam_file = pysam.Samfile(FILE_NAME, "rb")
output_file = open("mate_pair_distances.tsv","w")
iter= sam_file.fetch("REGION_OF_INTEREST",1,500)
for x in iter:
line = (str(x))
cols = line.rstrip("\n").split("\t")
print(cols)
if cols[2] == "REGION_OF_INTEREST":
mate_dist=int(cols[8])
position = int(cols[3])

[dpryan]

Why are you converting the convenient to use AlignedRead class to a string? Just:

sam_file = pysam.Samfile(FILE_NAME, "rb")
output_file = open("mate_pair_distances.tsv","w")
for alignment in  sam_file.fetch("REGION_OF_INTEREST",1,500) :
    isize=alignment.tlen
    ...

[dpryan]

Note that you might want to toss an abs() in there, size the insert size will be negative for one of the mates.

[Tom.Booker]

Thanks! That's cleared it right up.