hiv-zipcode-tools/process-bam.py at master · KiddLab/hiv-zipcode-tools · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
#!/usr/bin/env python2
import sys
import pysam

from optparse import OptionParser

###############################################################################
USAGE = """
python process-bam.py   --bam <path for bam file from insertion mapping>

Will output hit set to file in same dir.


"""
parser = OptionParser(USAGE)
parser.add_option('--bam',dest='bamFile', help = 'bam file')

(options, args) = parser.parse_args()


if options.bamFile is None:
    parser.error('bamFile not given')
###############################################################################


print options.bamFile

of = options.bamFile + '.align-parse.txt'
print 'writing output to',of
outFile = open(of,'w')


nl = ['#zipCode','readName','insChrom','insBP','insStrand','insMapQ','shearChrom','shearBP','ShearReadOrient','shearMapQ']
nl = '\t'.join(nl) + '\n'
outFile.write(nl)

bamFile = pysam.AlignmentFile(options.bamFile, "rb")
# get r1 and r2, primary...
r1 = None
r2 = None
for read in bamFile:
    if read.is_secondary is True:
        continue
    if read.is_read1 is True:
        if r1 == None:
            r1 = read
        else:
            print 'found r1 when have r1...???'
            print r1
            print read
            sys.exit()
    if read.is_read2 is True:
        if r2 == None:
            r2 = read
        else:
            print 'found r2 when have r2...???'
            print r2
            print read
            sys.exit()
    if r1 == None or r2 == None:
        continue # get next read
    else:
        # do the analysis
        name1 = r1.query_name
        name2 = r2.query_name
        if name1 != name2:
            print 'names do not match'
            print r1
            print r2
            sys.exit()

        # r1 should be to the genome
        r1Chrom = ''
        r1Pos = ''
        r1Dir = ''
        r2Chrom = ''
        r2Pos = '' # shear point
        r2Dir = ''

        qName = r1.query_name
        zc = qName.split(':')[0]

        if r1.is_unmapped is True:
            r1Chrom = 'NA'
            r1Pos = '.'
            r1Dir = '.'
            insPoint = '.'
            insOrient = '.'
            mapQ = -1
        else:
            r1Chrom = r1.reference_name
            refPos = r1.get_reference_positions()
            r1Start = refPos[0] +1 # will work in 1 based...
            r1End = refPos[-1] + 1 # will work in 1 based...
            mapQ = r1.mapping_quality
            if r1.is_reverse is True:
                r1Dir = '-'  # reverse strand, so in positive dir   # this is the dir of HIV insertion relative to genome
                insPoint = r1End
                insOrient = 'fwd'
            else:
                r1Dir = '+'  # positive strand, so in reverse dir
                insPoint = r1Start
                insOrient = 'rev'

        # now do to R2...
        if r2.is_unmapped is True:
            r2Chrom = 'NA'
            r2Pos = '.'
            r2Dir = '.'
            shearPoint = 0
            mapQR2 = -1
            shearOrient = '.'
        else:
            r2Chrom = r2.reference_name
            refPos = r2.get_reference_positions()
            r2Start = refPos[0] +1 # will work in 1 based...
            r2End = refPos[-1] + 1 # will work in 1 based...
            mapQR2 = r1.mapping_quality

            if r2.is_reverse is True:  # this is true shear orientation
                r2Dir = '-'  # reverse strand, so in positive dir
                shearOrient = 'rev'
                shearPoint = r2End
            else:
                r2Dir = '+'  # positive strand, so in reverse dir
                shearOrient = 'fwd'
                shearPoint = r2Start


#        print 'info!!!'
#        print r1
#        print r1Chrom,r1Start,r1End,r1Dir,insPoint,insOrient,mapQ

#        print 'r2!!!!'
#        print r2
#        print r2Chrom,r2Start,r2End,r2Dir,shearPoint,r2Dir,mapQR2

        nl = [zc,qName,r1Chrom,insPoint,insOrient,mapQ,r2Chrom,shearPoint,shearOrient,mapQR2]
        nl = [str(j) for j in nl ]
        nl = '\t'.join(nl) + '\n'
        outFile.write(nl)

       ###### end do the analysis
        r1 = None
        r2 = None