How to read an mpileup file from Samtools into Python
1
1
Entering edit mode
3.3 years ago
M.O.L.S ▴ 40

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.

samtools mpileup python next-gen • 2.0k views
ADD COMMENT
0
Entering edit mode

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.

ADD REPLY
1
Entering edit mode
3.2 years ago
M.O.L.S ▴ 40
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")
ADD COMMENT
1
Entering edit mode

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

ADD REPLY
0
Entering edit mode

I did not know that 3.1 years ago. Thanks !

ADD REPLY

Login before adding your answer.

Traffic: 2503 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6