1#=========================================================================== 2# 3# PUBLIC DOMAIN NOTICE 4# National Center for Biotechnology Information 5# 6# This software/database is a "United States Government Work" under the 7# terms of the United States Copyright Act. It was written as part of 8# the author's official duties as a United States Government employee and 9# thus cannot be copyrighted. This software/database is freely available 10# to the public for use. The National Library of Medicine and the U.S. 11# Government have not placed any restriction on its use or reproduction. 12# 13# Although all reasonable efforts have been taken to ensure the accuracy 14# and reliability of the software and data, the NLM and the U.S. 15# Government do not and cannot warrant the performance or results that 16# may be obtained by using this software or data. The NLM and the U.S. 17# Government disclaim all warranties, express or implied, including 18# warranties of performance, merchantability or fitness for any particular 19# purpose. 20# 21# Please cite the author in any work or product based on this material. 22# 23#=========================================================================== 24# 25 26import sys 27import traceback 28 29from ngs import NGS 30from ngs.ErrorMsg import ErrorMsg 31from ngs.ReadCollection import ReadCollection 32from ngs.Reference import Reference 33from ngs.PileupEvent import PileupEvent 34 35 36def run(acc, refName, start, stop): 37 # open requested accession using SRA implementation of the API 38 with NGS.openReadCollection(acc) as run: 39 run_name = run.getName() 40 41 # get requested reference 42 with run.getReference(refName) as ref: 43 # start iterator on requested range 44 with ref.getPileupSlice(start-1, stop-start+1) as it: 45 i = 0 46 while it.nextPileup(): 47 qual = "" 48 base = "" 49 50 line = "{}\t{}\t{}\t{}".format( 51 it.getReferenceSpec(), 52 it.getReferencePosition()+1, 53 it.getReferenceBase(), 54 it.getPileupDepth(), 55 ) 56 while it.nextPileupEvent(): 57 e = it.getEventType() 58 if (e & PileupEvent.alignment_start ) != 0: 59 base = base + '^' 60 base = base + chr(it.getMappingQuality() + 33) 61 62 if ( e & PileupEvent.insertion ) != 0: 63 base = base + '+' 64 ibases= it.getInsertionBases() 65 c = len(ibases) 66 base = base + str(c) 67 68 if ( e & PileupEvent.alignment_minus_strand ) != 0: 69 base = base + ibases.lower() 70 else: 71 base = base + ibases 72 73 evt = e & 7 74 75 if ( e & PileupEvent.alignment_minus_strand ) != 0: 76 if evt == PileupEvent.deletion: 77 base = base + '<' 78 elif evt == PileupEvent.match: 79 base = base + ',' 80 elif evt == PileupEvent.mismatch: 81 base = base + str(it.getAlignmentBase()).lower() 82 else: 83 if evt == PileupEvent.deletion: 84 base = base + '>' 85 elif evt == PileupEvent.match: 86 base = base + '.' 87 elif evt == PileupEvent.mismatch: 88 base = base + str(it.getAlignmentBase()).upper() 89 90 if ( e & PileupEvent.alignment_stop ) != 0: 91 base = base + '$' 92 93 qual = qual + it.getAlignmentQuality() 94 95 i += 1 96 print ("{}\t{}\t{}".format(line, base, qual)) 97 print ("Read {} pileups for {}".format(i, run_name)) 98 99if len(sys.argv) != 5: 100 print ("Usage: PileupTest accession reference start stop\n") 101 exit(1) 102else: 103 try: 104 run(sys.argv[1], sys.argv[2], int(sys.argv[3]), int(sys.argv[4])) 105 except ErrorMsg as x: 106 print (x) 107 traceback.print_exc() 108 # x.printStackTrace - not implemented 109 exit(1) 110 except BaseException as x: 111 traceback.print_exc() 112 exit(1) 113