diff --git a/atomicModel.py b/atomicModel.py
index 78fc9f0432fd8c996b7b4b09002005cf0ef4b536..0b3fc2c997e4aa41e238e7bd8f0cd15242ebf608 100644
--- a/atomicModel.py
+++ b/atomicModel.py
@@ -378,6 +378,7 @@ class atomicModel():
last = s
strandSegments = [(i+1,j) for i,j in zip([lo-1]+split,split+[hi])]
+
## Find zIdx for each nucleotide
zIdxs = []
for l,h in strandSegments:
@@ -391,25 +392,34 @@ class atomicModel():
assert(len(zIdxs) == numNt)
- baseIds = range(numNt)
if isFwd:
baseToSeg = baseToSegFwd[hid]
else:
- baseIds = reversed(baseIds)
baseToSeg = baseToSegRev[hid]
zIdxs = reversed(zIdxs)
+ ## Find strandOccupancy index for each nucleotide
+ strandOccIdx = []
+ for l,h in strandSegments:
+ for i in range(l,h+1):
+ nts = 1
+ for ins in strand.insertionsOnStrand(i,i):
+ nts += ins.length()
+ baseToSeg[i] = []
+ for j in range(nts):
+ strandOccIdx.append(i)
+ if not isFwd:
+ strandOccIdx = reversed(strandOccIdx)
+
## Add nucleotides
- for i,s,zIdx in zip(baseIds,seqs,zIdxs):
+ for s,zIdx,idx in zip(seqs,zIdxs,strandOccIdx):
p = zIdxToPos(zIdx)
a = zIdxToAngle(zIdx)
- lastNt = seg.addNucleotide(hid, i, zIdx, s, isFwd,
+ lastNt = seg.addNucleotide(hid, idx, zIdx, s, isFwd,
p, a, lastNt)
- # baseToSeg[i+1] = [lastNt]
- assert(zIdx not in baseToSeg)
- baseToSeg[zIdx] = [lastNt]
+ baseToSeg[idx].append(lastNt)
## 2) Find basepairs and stacks
for hid in range(numHID):
@@ -418,7 +428,6 @@ class atomicModel():
continue
ends1,ends2 = self._helixStrandsToEnds(helixStrands)
-
strandOccupancies = [ [x for i in range(0,len(e),2)
for x in range(e[i],e[i+1]+1)]
for e in (ends1,ends2) ]
@@ -435,7 +444,7 @@ class atomicModel():
reversed(baseToSegRev[hid][i])):
self._pairBases(nt1,nt2)
- if prevBasepair is not None and prevBasepair[0] == i-1:
+ if prevBasepair is not None and i - prevBasepair[0] <= 3: # TODO: find better way of locating stacks
atomicModel._stackBases( prevBasepair[1], nt1)
atomicModel._stackBases( nt2, prevBasepair[2])
# else:
@@ -450,16 +459,6 @@ class atomicModel():
below.addNodeAbove(above, sep)
above.addNodeBelow(below, sep)
- def _addCrossover(self, hid1, hid2, xo):
- zid1, zid2, isFwd1, isFwd2 = xo
- node1 = self.helices[hid1].nodes[zid1]
- node2 = self.helices[hid2].nodes[zid2]
-
- ## TODO add polarity
- polarity = 0
- node1.addXover(node2, isFwd1)
- node2.addXover(node1, isFwd2)
-
def _getNeighborHelixDict(part):
props, origins = part.helixPropertiesAndOrigins()
neighborHelices = dict()
@@ -860,7 +859,8 @@ class atomicModel():
enmTemplate = enmTemplateHC
else:
raise Exception("Lattice type '%s' not supported" % self.latticeType)
-
+ noStackPrime = 0
+ noBasepair = 0
with open("%s.exb" % prefix,'w') as fh:
# natoms=0
for seg in self.segments:
@@ -875,12 +875,17 @@ class atomicModel():
if nt2 is not None and nt2.firstAtomicIndex > nt1.firstAtomicIndex:
other.append((nt2,'paircross'))
+ else:
+ noBasepair += 1
+
nt2 = nt1.stack3prime
if nt2 is not None:
other.append((nt2,'stack'))
nt2 = nt2.basepair
if nt2 is not None and nt2.firstAtomicIndex > nt1.firstAtomicIndex:
other.append((nt2,'cross'))
+ else:
+ noStackPrime += 1
a1 = nt1.atoms(transform=False)
for nt2,key in other:
@@ -888,10 +893,11 @@ class atomicModel():
for n1, n2, d in enmTemplate[key]:
d = float(d)
a2 = nt2.atoms(transform=False)
- try:
- i,j = [a.index(n)-1 for a,n in zip((a1,a2),(n1,n2))]
- except:
- continue
+ i,j = [a.index(n)-1 for a,n in zip((a1,a2),(n1,n2))]
+ # try:
+ # i,j = [a.index(n)-1 for a,n in zip((a1,a2),(n1,n2))]
+ # except:
+ # continue
k = 0.1
if self.latticeType == 'honeycomb':
@@ -905,6 +911,8 @@ class atomicModel():
fh.write("bond %d %d %f %.2f\n" % (i,j,k,d))
+ print("NO STACKS found for : %d" % noStackPrime)
+ print("NO BPs: %d" % noBasepair)
props, origins = self.part.helixPropertiesAndOrigins()
origins = [np.array(o) for o in origins]
@@ -956,7 +964,7 @@ class atomicModel():
if nt.basepair is None: continue
hid = nt.prop['helixId']
if hid not in xo:
- print("WARNING: writeENM: helix %d not in xo dict" % hid)
+ # print("WARNING: writeENM: helix %d not in xo dict" % hid)
continue
zidx = nt.prop['zIdx']
isFwd = nt.prop['isFwd']
@@ -964,7 +972,7 @@ class atomicModel():
for hid2 in neighborHelices[hid]:
# if hid2 not in xo[hid]:
if hid2 not in helixIds:
- print("WARNING: writeENM: helix %d not in helixIds" % hid2) # xo[%d] dict" % (hid2,hid))
+ # print("WARNING: writeENM: helix %d not in helixIds" % hid2) # xo[%d] dict" % (hid2,hid))
continue
dist = min([np.abs(z-zidx) for z in xo[hid][hid2]])
key = (hid2,zidx,isFwd)
diff --git a/run.py b/run.py
index 8e58251f00445b3dd87f69acb781abc23e002e28..5f2c71be78a77b6fef5ea617bd72cdf415e060c8 100644
--- a/run.py
+++ b/run.py
@@ -95,7 +95,7 @@ def simulateDesign(inFile, outPrefix, numCoarsestSteps=5000000, gpu=0, shiftEule
finerModel = beadModelTwist(part, twistPersistenceLength=twistLp,
maxBpsPerDNode=1, maxNtsPerSNode=1)
- finestModel = atomicModel(part, ntScale=0.5)
+ finestModel = atomicModel(part, ntScale=0.25)
coarsestModelNoLongBonds._removeIntrahelicalConnectionsBeyond(100)
coarsestModelNoLongBonds._removeCrossoversBeyond(100)