[2] | 1 | #!/usr/bin/env python |
---|
| 2 | |
---|
| 3 | """ |
---|
| 4 | Adapted from bx/scripts/qv_to_bqv.py |
---|
| 5 | |
---|
| 6 | Convert a qual (qv) file to several BinnedArray files for fast seek. |
---|
| 7 | This script takes approximately 4 seconds per 1 million base pairs. |
---|
| 8 | |
---|
| 9 | The input format is fasta style quality -- fasta headers followed by |
---|
| 10 | whitespace separated integers. |
---|
| 11 | |
---|
| 12 | usage: %prog qual_file output_file |
---|
| 13 | """ |
---|
| 14 | |
---|
| 15 | import pkg_resources |
---|
| 16 | pkg_resources.require( "bx-python" ) |
---|
| 17 | pkg_resources.require( "numpy" ) |
---|
| 18 | import string |
---|
| 19 | import psyco_full |
---|
| 20 | import sys, re, os, tempfile |
---|
| 21 | from bx.binned_array import BinnedArrayWriter |
---|
| 22 | from bx.cookbook import * |
---|
| 23 | import fileinput |
---|
| 24 | |
---|
| 25 | def load_scores_ba_dir( dir ): |
---|
| 26 | """ |
---|
| 27 | Return a dict-like object (keyed by chromosome) that returns |
---|
| 28 | FileBinnedArray objects created from "key.ba" files in `dir` |
---|
| 29 | """ |
---|
| 30 | return FileBinnedArrayDir( dir ) |
---|
| 31 | |
---|
| 32 | def main(): |
---|
| 33 | args = sys.argv[1:] |
---|
| 34 | try: |
---|
| 35 | qual_file_dir = args[0] |
---|
| 36 | #mydir="/home/gua110/Desktop/chimp_quality_scores/chr22.qa" |
---|
| 37 | mydir="/home/gua110/Desktop/rhesus_quality_scores/rheMac2.qual.qv" |
---|
| 38 | qual_file_dir = mydir.replace(mydir.split("/")[-1], "") |
---|
| 39 | output_file = args[ 1 ] |
---|
| 40 | fo = open(output_file,"w") |
---|
| 41 | except: |
---|
| 42 | print "usage: qual_file output_file" |
---|
| 43 | sys.exit() |
---|
| 44 | |
---|
| 45 | tmpfile = tempfile.NamedTemporaryFile() |
---|
| 46 | cmdline = "ls " + qual_file_dir + "*.qa | cat >> " + tmpfile.name |
---|
| 47 | os.system (cmdline) |
---|
| 48 | for qual_file in tmpfile.readlines(): |
---|
| 49 | qual = fileinput.FileInput( qual_file.strip() ) |
---|
| 50 | outfile = None |
---|
| 51 | outbin = None |
---|
| 52 | base_count = 0 |
---|
| 53 | mega_count = 0 |
---|
| 54 | |
---|
| 55 | for line in qual: |
---|
| 56 | line = line.rstrip("\r\n") |
---|
| 57 | if line.startswith(">"): |
---|
| 58 | # close old |
---|
| 59 | if outbin and outfile: |
---|
| 60 | print "\nFinished region " + region + " at " + str(base_count) + " base pairs." |
---|
| 61 | outbin.finish() |
---|
| 62 | outfile.close() |
---|
| 63 | # start new file |
---|
| 64 | region = line.lstrip(">") |
---|
| 65 | #outfname = output_file + "." + region + ".bqv" #CHANGED |
---|
| 66 | outfname = qual_file.strip() + ".bqv" |
---|
| 67 | print >>fo, "Writing region " + region + " to file " + outfname |
---|
| 68 | outfile = open( outfname , "wb") |
---|
| 69 | outbin = BinnedArrayWriter(outfile, typecode='b', default=0) |
---|
| 70 | base_count = 0 |
---|
| 71 | mega_count = 0 |
---|
| 72 | else: |
---|
| 73 | if outfile and outbin: |
---|
| 74 | nums = line.split() |
---|
| 75 | for val in nums: |
---|
| 76 | outval = int(val) |
---|
| 77 | assert outval <= 255 and outval >= 0 |
---|
| 78 | outbin.write(outval) |
---|
| 79 | base_count += 1 |
---|
| 80 | if (mega_count * 1000000) <= base_count: |
---|
| 81 | sys.stdout.write(str(mega_count)+" ") |
---|
| 82 | sys.stdout.flush() |
---|
| 83 | mega_count = base_count // 1000000 + 1 |
---|
| 84 | if outbin and outfile: |
---|
| 85 | print "\nFinished region " + region + " at " + str(base_count) + " base pairs." |
---|
| 86 | outbin.finish() |
---|
| 87 | outfile.close() |
---|
| 88 | |
---|
| 89 | if __name__ == "__main__": |
---|
| 90 | main() |
---|