Question: How to read an mpileup file from Samtools into Python
1
gravatar for M.O.L.S
17 months ago by
M.O.L.S10
M.O.L.S10 wrote:

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.

mpileup samtools next-gen python • 1.0k views
ADD COMMENTlink modified 17 months ago • written 17 months ago by M.O.L.S10

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 REPLYlink written 17 months ago by Pierre Lindenbaum127k
1
gravatar for M.O.L.S
17 months ago by
M.O.L.S10
M.O.L.S10 wrote:
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 COMMENTlink modified 17 months ago by finswimmer13k • written 17 months ago by M.O.L.S10
1

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 REPLYlink modified 17 months ago • written 17 months ago by Carambakaracho2.0k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1120 users visited in the last hour