How to read an mpileup file from Samtools into Python

Question

[M.O.L.S]

How do I read an mpileup file into Python?

I have tried several times, but it is just not working. I have a txt file instead on an mpileup file.

[Pierre Lindenbaum]

I have tried several times, but it is just not working. I have a txt file instead on an mpileup file.

that's not clear. Explain what you're doing. Show us the code.

Answer 1

[M.O.L.S]

import pysam

# Read in a BAM file
bamfile = pysam.AlignmentFile("outputFile.bam","rb")
for line in bamfile:
    print(line)

# Read in a SAM file 
samfile = pysam.AlignmentFile("outputFile.sam","r")
for line in samfile:
    print(line)

#Read in a CRAM file. 
CRAMfile = pysam.AlignmentFile("outputfile.cram", "rc")

#Create the iterable and print the first line 
it = iter(bamfile)
print(next(it))

#itt = iter(samfile)
print(next(itt))

# A different way to open the mpileup file.         
f = open("/Users/m.o.l.s/outputFile.mpileup","rt")

line_1 = f.readline()
print(line_1)

# Split up the strings based on tabs
individuals = line_1.split("\t")
print(individuals)

# The value we are interested in is in the 3rd index position
The_Coverage = individuals[3]
print(The_Coverage)

# If this number is less than 10, the row needs to be displayed
if The_Coverage <= "10":
    print(The_Coverage)
else:
    print("The Coverage is not less than or equal to 10")

[Carambakaracho]

FYI, it's considered good practice to close file handles when not needed anymore. It's not important in your code above, but just in case you plan getting serious with this

see first example from the API documentation

import pysam
samfile = pysam.AlignmentFile("ex1.bam", "rb")

for read in samfile.fetch('chr1', 100, 120):
     print read

samfile.close() # note the close operation

More reading on opening files without requiring explicit close operations (not sure how to apply this to the pysam API, though) https://docs.python.org/2.5/whatsnew/pep-343.html

[M.O.L.S]

I did not know that 3.1 years ago. Thanks !