Minimap2 -ax map-pb doesn't output tlen field

Question

I have used minimap2 to map some pacbio reads to a reference genome. I would like to know the "insert size" (true length of sequence) relative to the reference. More specifically, I want to know the first base position in the reference where the read maps and the last base position in the reference where the mapping ends.

As I understand it I need the TLEN field from bam files (sam format field 9). However, all the alignment files from minimap2 have 0 in the TLEN field.

Does anyone know how to get pacbio alignments with the TLEN field? I am also happy using bwa or another aligner.

The command I used was minimap2 -ax map-pb {fasta} {reads} | samtools view -h -b - | samtools sort - > myfile.bam.

Answer 1

The TLEN for single-end reads is defined as 0 in the SAM specification, which is why this happens. You can still get the end position of the alignment by parsing the CIGAR string. An example using the pysam package in python is below:

import pysam

bam = pysam.AlignmentFile("myfile.bam")
for b in bam.fetch():
    print("This alignment spans {}:{}-{}".format(b.reference_name, b.reference_start+1, b.reference_end))

The reference_end value is automatically calculated from the CIGAR string.