-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathBIRDbath1.py
executable file
·199 lines (184 loc) · 7.17 KB
/
BIRDbath1.py
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
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
#!/usr/bin/env python
#=========================================================================
# This is OPEN SOURCE SOFTWARE governed by the Gnu General Public
# License (GPL) version 3, as described at www.opensource.org.
# Copyright (C)2018 William H. Majoros ([email protected])
#=========================================================================
from __future__ import (absolute_import, division, print_function,
unicode_literals, generators, nested_scopes, with_statement)
from builtins import (bytes, dict, int, list, object, range, str, ascii,
chr, hex, input, next, oct, open, pow, round, super, filter, map, zip)
# The above imports should allow this program to run in both Python 2 and
# Python 3. You might need to update your version of module "future".
import sys
import os
import math
import ProgramName
from Rex import Rex
rex=Rex()
import TempFilename
import getopt
from StanParser import StanParser
from PooledParser import PooledParser
from Stan import Stan
DEBUG=False
WARMUP=1000
ALPHA=0.05
#POP_CONC=1500 ### THIS NEEDS TO BE ESTIMATED FROM EACH DATA SET!
STDERR=TempFilename.generate(".stderr")
INPUT_FILE=TempFilename.generate(".staninputs")
INIT_FILE=TempFilename.generate(".staninit")
OUTPUT_TEMP=TempFilename.generate(".stanoutputs")
def printFields(fields,hFile):
numFields=len(fields)
for i in range(7,numFields):
print(i-6,"=",fields[i],sep="",end="",file=hFile)
if(i<numFields-1): print("\t",end="",file=hFile)
print(file=hFile)
def getFieldIndex(label,fields):
numFields=len(fields)
index=None
for i in range(7,numFields):
if(fields[i]==label): index=i
return index
def writeToFile(fields,OUT):
numFields=len(fields)
for i in range(7,numFields):
print(fields[i],end="",file=OUT)
if(i<numFields-1): print("\t",end="",file=OUT)
print(file=OUT)
def writeInitializationFile(stan,variant,filename):
OUT=open(filename,"wt")
print("theta <- 1",file=OUT)
freqs=variant.getFreqs()
numPools=variant.numPools()
stan.writeOneDimArray("p",freqs,numPools,OUT)
maxRnaReps=variant.getMaxRnaReps()
#qiInit=[freqs]*maxRnaReps
qiInit=makeQiArray(freqs,maxRnaReps)
#print("pinit=",freqs)
#print("qiInit=",qiInit)
#print("numPools=",numPools,", maxRnaReps=",maxRnaReps)
stan.writeTwoDimArray("qi",qiInit,numPools,maxRnaReps,OUT)
print("c <- 100",file=OUT)
print("s <- 1",file=OUT)
OUT.close()
def makeQiArray(freqs,maxRnaReps):
array=[]
for f in freqs:
array.append([f]*maxRnaReps)
return array
def writeInputsFile(stan,variant,filename):
OUT=open(filename,"wt")
numPools=variant.numPools()
print("N_POOLS <-",numPools,file=OUT)
maxDnaReps=variant.getMaxDnaReps()
maxRnaReps=variant.getMaxRnaReps()
print("MAX_DNA <- ",maxDnaReps,file=OUT)
print("MAX_RNA <- ",maxRnaReps,file=OUT)
print("pop_conc <- ",POP_CONC,file=OUT)
freqs=variant.getFreqs()
stan.writeOneDimArray("pop_freq",freqs,numPools,OUT)
dnaReps=variant.getDnaReps()
rnaReps=variant.getRnaReps()
stan.writeOneDimArray("N_DNA",dnaReps,numPools,OUT)
stan.writeOneDimArray("N_RNA",rnaReps,numPools,OUT)
dnaAltCounts=[[rep.alt for rep in pool.DNA] for pool in variant.pools]
dnaRefCounts=[[rep.ref for rep in pool.DNA] for pool in variant.pools]
rnaAltCounts=[[rep.alt for rep in pool.RNA] for pool in variant.pools]
rnaRefCounts=[[rep.ref for rep in pool.RNA] for pool in variant.pools]
dnaAltCounts=expandArray(dnaAltCounts,maxDnaReps)
dnaRefCounts=expandArray(dnaRefCounts,maxDnaReps)
rnaAltCounts=expandArray(rnaAltCounts,maxRnaReps)
rnaRefCounts=expandArray(rnaRefCounts,maxRnaReps)
stan.writeTwoDimArray("a",dnaAltCounts,numPools,maxDnaReps,OUT)
stan.writeTwoDimArray("b",dnaRefCounts,numPools,maxDnaReps,OUT)
stan.writeTwoDimArray("k",rnaAltCounts,numPools,maxRnaReps,OUT)
stan.writeTwoDimArray("m",rnaRefCounts,numPools,maxRnaReps,OUT)
OUT.close()
def expandArray(array,desiredSize):
newArray=[]
numPools=len(array)
for i in range(numPools):
oldPool=array[i]
newPool=[x for x in oldPool]
oldPoolSize=len(oldPool)
for j in range(oldPoolSize,desiredSize):
newPool.append(0)
newArray.append(newPool)
return newArray
def runVariant(stan,variant,numSamples,outfile):
# Write inputs file for STAN
writeInputsFile(stan,variant,INPUT_FILE)
writeInitializationFile(stan,variant,INIT_FILE)
# Run STAN model
cmd=stan.getCmd(WARMUP,numSamples,INPUT_FILE,OUTPUT_TEMP,STDERR,INIT_FILE)
if(DEBUG):
print(cmd)
exit()
os.system(cmd)
# Parse MCMC output
parser=StanParser(OUTPUT_TEMP)
thetas=parser.getVariable("theta")
return (thetas,parser)
def summarize(parser,thetas,ID,minRight):
(median,CI_left,CI_right)=parser.getMedianAndCI(1.0-ALPHA,"theta")
maxLeft=1.0/minRight
leftP=parser.getLeftTail("theta",maxLeft)
rightP=parser.getRightTail("theta",minRight)
Preg=leftP if leftP>rightP else rightP
print(ID,round(median,3),round(CI_left,3),round(CI_right,3),
round(Preg,3),sep="\t")
#=========================================================================
# main()
#=========================================================================
(options,args)=getopt.getopt(sys.argv[1:],"s:t:")
if(len(args)!=7):
exit(ProgramName.get()+" [-s stanfile] [-t thetafile] <model> <min-effect> <beta-concentration-parm> <input.essex> <output.txt> <#MCMC-samples> <firstVariant-lastVariant>\n -s = save raw STAN file\n -t = save theta samples\n variant range is zero-based and inclusive\n min-effect (lambda) must be >= 1\n")
(model,minEffect,POP_CONC,inFile,outfile,numSamples,numVariants)=args
stanFile=None
thetaFile=None
for pair in options:
(key,value)=pair
if(key=="-s"): stanFile=value
if(key=="-t"): thetaFile=value
if(not rex.find(r"(\d+)-(\d+)",numVariants)):
exit(numVariants+": specify range of variants: first-last")
firstIndex=int(rex[1])
lastIndex=int(rex[2])
minEffect=float(minEffect)
POP_CONC=float(POP_CONC)
if(minEffect<1): raise Exception("Min-effect must be >= 1")
THETA=None
if(thetaFile is not None): THETA=open(thetaFile,"wt")
stan=Stan(model)
# Process all input lines, each line = one variant (one MCMC run)
thetaIndex=None
variantIndex=0
pooledParser=PooledParser(inFile)
while(True):
variant=pooledParser.nextVariant()
if(variant is None): break
# Check whether this variant is in the range to be processed
if(variantIndex<firstIndex):
variantIndex+=1
continue
elif(variantIndex>lastIndex): break
keep=variant.dropHomozygousPools()
if(not keep): continue
(thetas,stanParser)=runVariant(stan,variant,numSamples,outfile)
if(thetas is None): continue
summarize(stanParser,thetas,variant.ID,minEffect)
variantIndex+=1
if(THETA is not None):
for i in range(len(thetas)):
print(thetas[i],file=THETA,end="")
if(i<len(thetas)): print("\t",file=THETA,end="")
print(file=THETA)
os.remove(STDERR)
os.remove(INPUT_FILE)
if(stanFile is None):
os.remove(OUTPUT_TEMP)
else:
os.system("cp "+OUTPUT_TEMP+" "+stanFile)
if(THETA is not None): THETA.close()