Dear all,

I have indexed bam file, and I want to print all alignments on chromosome "2" using pysam.

My code is 

    import pysam
    bam = pysam.AlignmentFile("Aligned.sortedByCoord_rep1.out.bam", "rb")
    for line in bam.fetch("2"):
         print line

What I get is:

    D00733:162:CADM2ANXX:2:1311:7466:63232	163	18	3328	255	1S47M2S	18	3419	47	TTGCTTAGTGTCCGAAATACCATCCTCAAGGCTAAGAACTAAATCGATTA	array('B', [33, 33, 33, 33, 33, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38])	[('NH', 1), ('HI', 1), ('AS', 70), ('nM', 10)]
    D00733:162:CADM2ANXX:2:1215:11810:89291	163	18	3356	255	44M6S	18	3441	44	GGCTAAGAACTAAATCGATTATTCTGGCTCGTAACGCATATAATATGGCA	array('B', [33, 33, 33, 32, 33, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 35, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 37, 38, 38, 38, 38, 38, 38, 33])	[('NH', 1), ('HI', 1), ('AS', 70), ('nM', 10)]

So there are 2 questions:
 1. Why pysam does not preserve format of bam file?
 2. Why does it print out chromosome 18 instead of 2?

I would appreciate if somebody points out on what am I doing wrong.
Thanks
Hi, you are not wrong, you are doing it correctly :)

The AlignedSegment object is independent of a SAM-file. The AlignedSegment contains only an index of the chromosomal identifier (tid), and, thus, it does not contain the chromosome number information. It means that the number you see as 18 in the field 3 of the line is the number corresponding to chromosome 2 in the sequence dictionary. You can check that by doing the following:

    import pysam
    bam = pysam.AlignmentFile("Aligned.sortedByCoord_rep1.out.bam", "rb")
    for line in bam.fetch("2"):
        print line
        print line.reference_id
        print bam.get_reference_name(line.reference_id)


As you might have notice, I have added two lines to your previous code:

 1.`line.reference_id` contains the index of the reference sequence in the sequence dictionary 
 
2.`bam.get_reference_name(line.reference_id)` will return the reference name corresponding to numerical tid. In your case chromosome 2.








