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