# Strawberry Sort: clean version
**Note: this note is to be maintained by bj only, please do not muddle with it. It will contain the current "best version" of the full RGxHapil strawberry pipeline.**
Index:
[toc]
---
## Step 1: DBG and anchor runs
The idea behind this step is two-fold:
- Create an "inventory" of the 31-mer spectrum, that will be used on subsequent steps to define what a "unique anchor" is. The unique anchors are unitig regions that have all their 31-mers between MIN_UKMER_FREQ and MAX_UKMER_FREQ (i.e a "single peak").
- To make the graph simpler to start with, we create unitigs from a
```python=
import sys
sys.path.append('/hpc-home/ggarcia/git_sources/bsg/build/SDGpython')
import SDGpython as SDG
import os
MIN_UKMER_FREQ=57
MAX_UKMER_FREQ=135
K=31
## SDG-DBG assembly k=63
## Min coverage is important to reduce the noice
GK=63
MIN_CVG=30
NUM_BATCHES=8
ws=SDG.WorkSpace()
peds=ws.add_paired_reads_datastore('./data/strbrg_pe.prseq', 'pe')
gm=SDG.GraphMaker(ws.sdg)
gm.new_graph_from_paired_datastore(peds,K,MIN_CVG,NUM_BATCHES); #this is so the initial kmer count includes all 31-mers
print("First graph done")
kc = ws.add_kmer_counter('pek31', K)
kc.add_count('pe', peds)
print("kmers counted")
gm.new_graph_from_paired_datastore(peds,GK,MIN_CVG,NUM_BATCHES); #this is the real DBG k
print("second graph done")
ws.sdg.write_to_gfa1("./00-sdg_dbg_graph_MC%s.gfa" %MIN_CVG)
gc=SDG.GraphContigger(ws)
gc.clip_tips(2*GK-1,3)
print("Remaining nodes: %s" %(len(ws.sdg.get_all_nodeviews())))
kc.update_graph_counts()
ws.dump("./00-sdg_dbg_graph_MC%s.sdgws" %MIN_CVG)
## Check step1
passed=True
ccc=zip(kc.count_spectra('pe',present_in_graph=False),kc.count_spectra('pe',present_in_graph=True))
for i,x in list(enumerate(ccc))[MIN_UKMER_FREQ:MAX_UKMER_FREQ+1]:
if x[0]!=x[1]:
print(x[0]-x[1],"Kmers are missing at f=",i," -> ",x)
passed=False
if passed: print("Check passed!!!!")
else: print("Check FAILED!!!!")
## Make unique kmer runs
seqs = gc.contig_reduction_to_unique_kmers(MIN_UKMER_FREQ,MAX_UKMER_FREQ)
print(len(seqs),"sequences for anchoring")
with open("./pseudoLong2.fasta","w") as f:
for i, s in enumerate(seqs):
f.write(">seq_%s\n%s\n"%(i,s))
```
---
## Step 2: workspace, KmerCounters, and nanopore mapping
k=31 index, every unique mapping is kept.
---
## Step 3: filtering nanopore reads by F1 pattern
```python=
from collections import Counter
def best_patterns(patterns,verbose=False):
if len(patterns)<10:
if verbose: print('less than 10 matches, giving up!')
return []
all_1s='1'*len(patterns[0])
pc=Counter(patterns)
t=len(patterns)
if pc[all_1s]>=t*.75:
if verbose: print('Too many 1s, giving up!')
return []
t-=pc[all_1s]
pc[all_1s]=0
hc=sum([c for p,c in pc.items() if p[-1]=='1']) #count patterns with hapil
if hc>t*.7:
if verbose: print('Too many patterns with hapil, using full pattern')
fpc=pc
ft=t
else:
if verbose: print('Discarding',hc,'/',t,'patterns with hapil')
fpc=Counter()
ft=0
for p,c in pc.items():
if p[-1]!='1':
fpc[p]=c
ft+=c
w=fpc.most_common(1)[0][0]
if verbose: print("most common pattern:",fpc.most_common(1),'/',ft)
if fpc[w]>.60*ft:
if w[-1]=='1': return [w]
else: return [w,w[:-1]+'1']
else:
if verbose: print('Winner not 60% of filtered matches, giving up')
return []
def filter_read(rid,lrr,cf,verbose=False):
patterns=[cf.get_pattern(ws.sdg.get_nodeview(m.node)) for m in lrr.read_perfect_matches[rid]]
bp=best_patterns(patterns,verbose)
if verbose: print("winning patterns:",bp)
fm=[]
for i,p in enumerate(patterns):
if p in bp:
fm.append(lrr.read_perfect_matches[rid][i])
return fm
cf3=SDG.CountFilter('rgxh31','sdg',100,kc.list_names()[1:],[2 for x in kc.list_names()[1:-1]+[30]])
srpm2=[filter_read(rid,lrr,cf3) for rid in range(len(lrr.read_perfect_matches))]
for i in range(len(srpm2)):
lrr.read_perfect_matches[i]=srpm2[i]
lrr.dump('lrr_after_segfilterv2.dump')
```
---
## Step 4: multi-component thread dg partial orders
Load previos status:
```python=
import SDGpython as SDG
ws=SDG.WorkSpace()
ws.sdg.load_from_fasta('pseudoLong2.fasta')
lords = ws.add_long_reads_datastore('data/strbrg_nanopore.loseq','nanopore')
lrr=SDG.LongReadsRecruiter(ws.sdg,lords,31)
lrr.load('lords_kmerruns.mapped.data')
lrrf=SDG.LongReadsRecruiter(ws.sdg,lords,31)
lrrf.load('lrr_after_segfilterv2.dump')
```
```python=
def extend_solution(tgs, min_support=2,min_shared=10,min_new=10):
finished=False
for step in range(10000):
print("\nStarting expansion step",step)
s=sorted(tgs.sort_graph().items(), key=lambda x: x[1])
if s:
print("Current sorted size:",s[-1][1]-s[0][1])
else:
print("No sorted size")
# print("Current sorted size:",s[-1][1]-s[0][1])
added=0
for nodes in [[x[0] for x in s[:30]],[x[0] for x in s[-30:]]]:
rids=set()
for nid in nodes:
rids.update(tgs.rids_from_node(trg_nt.get_nodeview(nid)))
while rids:
to_add=0
largest_new=min_new-1
to_add_old=0
for rid in rids:
if rid in tgs.used_reads: continue
old,new=tgs.evaluate_read(rid)
if old>=min_shared and new>largest_new:
to_add=rid
largest_new=new
to_add_old=old
if to_add:
print("adding read %d with %d new and %d old nodes" % (to_add,largest_new,to_add_old))
rno=list(tgs.evaluate_read_nodeorder(to_add))
print(rno)
print("%d ordered, %d unordered and %d unused nodes" %(rno[0], rno[1], rno[2]))
if rno[1]:
print("skipping rid with unordered nodes")
rids.remove(to_add)
else:
tgs.add_read(to_add,min_support)
added+=1
break
else:
print("can't find any read to add!!!")
break
if not added:
print("can't find any read to add on both directions!!!")
finished=True
break
if finished: break
class CCSorter():
"""Takes an initial node and builds all TGSs needed to sort most long reads"""
def __init__(self,trg_nt,founding_node,min_nodes_per_tracked_read=30):
self.sorters=[]
self.trg_nt=trg_nt
self.founding_node=founding_node
self.master_tgs=SDG.TheGreedySorter(trg_nt,founding_node)
#for every read, check how many nodes it's got, and keep track of every read with more than X nodes
self.read_nodes={}
for rid,ends in self.master_tgs.read_ends.items():
rnodes=set()
nv=self.trg_nt.get_nodeview(ends[0])
while True:
try:
rnodes.add(abs(nv.node_id()))
nv=[x for x in nv.next() if x.support().id==rid][0].node()
except:
break
if len(rnodes)>=min_nodes_per_tracked_read:
self.read_nodes[rid]=rnodes
def read_status(self,rid):
used=unused=0
sortersets=[set(list(s.used_nodes)) for s in self.sorters]
for nid in self.read_nodes[rid]:
for s in sortersets:
if nid in s:
used+=1
break
else:
unused+=1
return used,unused
def run_new_sorter(self):
possible_starts=[]
for rid in self.read_nodes.keys():
used,unused=self.read_status(rid)
if used<3:
possible_starts.append((rid,used,unused))
possible_starts.sort(key=lambda x:(-x[1],x[2]),reverse=True)
if len(possible_starts)<10:
print("less than 10 possible starts, aborting")
return
self.sorters.append(SDG.TheGreedySorter(self.trg_nt,self.founding_node))
self.sorters[-1].start_from_read(possible_starts[0][0])
extend_solution(self.sorters[-1])
def read_usage(self):
possible_starts=[]
used_reads=unused_reads=0
for rid in self.read_nodes.keys():
used,unused=self.read_status(rid)
if used>10: used_reads+=1
else: unused_reads+=1
return used_reads,unused_reads
def make_all_orders(self):
old_order_count=-1
while len(self.sorters)>old_order_count:
old_order_count=len(self.sorters)
self.run_new_sorter()
```
```python=
lrrf.thread_reads(end_matches=1)
trg_nt=lrrf.dg_from_threads(False,True)
ccs100=trg_nt.get_all_connected_components(100)
multisorters=[]
for x in ccs100:
multisorters.append(CCSorter(trg_nt,list(x)[0]))
multisorters[-1].make_all_orders()
```
---
## Step 5: total order
```python=
lrr.thread_reads(end_matches=1)
trg_ntu=lrr.dg_from_threads(False,True,20)
tgsu=SDG.TheGreedySorter(trg_ntu)
rtgsu=SDG.TheGreedySorter(trg_ntu)
```
**TODO**
- Project read links across the order.
- Find unordered nodes already included in the order and pop them.
- For each read mapped to the node, walk fw and bw until hitting a node included in the order, and check the consistency.
```python=
from datetime import datetime
def write_connected(dg,filename):
dg.write_to_gfa1(filename,selected_nodes=[x.node_id() for x in dg.get_all_nodeviews(include_disconnected=False)])
def printlog(*args, **kwargs):
print(datetime.now().strftime('%H:%M:%S'),*args, **kwargs)
def local_happiness_assembler(tgs, rtgs, reads_tg, starting_rid, min_support=2, max_unhappy=0, max_disconnected_perc=.4,min_shared=10,min_new=10,max_unhappy_nodes_in_new_read=10,max_steps=10000000):
tgs.start_from_read(starting_rid,min_support)
for stepnumber in range(max_steps):
printlog("step #",stepnumber,"starting")
write_connected(tgs.dg,"step%d_a-starting.gfa"%(stepnumber))
# remove unhappy nodes from order
printlog("sorting graph")
order=tgs.sort_graph()
printlog("full order hapiness assessment")
nbh=[(x[0],x[1],SDG.assess_node_happiness(x[0],order,reads_tg)) for x in sorted(order.items(),key=lambda x:x[1])] # (nid, pos, (happy, unhappy, disconnected))
#first remove all nodes that are too disconnected: they won't get any better
printlog("removing disconnected")
unhappy_nodes=[]
removed_disconnected=removed_unhappy=0
for x in nbh:
if x[2][0]==0 or 1.0*x[2][2]/sum(x[2])>max_disconnected_perc:
#pop_node_from_all(tgs.dg,)
#print(x, "should go!")
tgs.remove_node(x[0])
removed_disconnected+=1
elif x[2][1]>max_unhappy:
unhappy_nodes.append(x)
printlog("removing unhappy")
#now remove any unhappy nodes
unhappy_nodes.sort(key=lambda x:x[2][1])
for x in unhappy_nodes:
h=SDG.assess_node_happiness(x[0],tgs.sort_graph(),reads_tg)
if h[1]>max_unhappy:
tgs.remove_node(x[0])
removed_unhappy+=1
printlog("%d disconnected and %d unhappy nodes removed from order"%(removed_disconnected,removed_unhappy))
write_connected(tgs.dg,"step%d_b-cleaned.gfa"%(stepnumber))
###TODO: THIS IS OLD CODE!
added=0
s=sorted(tgs.sort_graph().items(),key=lambda x:x[1])
for nodes in [[x[0] for x in s[:30]],[x[0] for x in s[-30:]]]:
rids=set()
for nid in nodes:
rids.update(tgs.rids_from_node(reads_tg.get_nodeview(nid)))
printlog("choosing among %d candidate reads" % len(rids))
while rids:
to_add=0
largest_new=min_new-1
to_add_old=0
for rid in rids:
if rid in tgs.used_reads: continue
old,new=tgs.evaluate_read(rid)
if old>=min_shared and new>largest_new:
to_add=rid
largest_new=new
to_add_old=old
if to_add:
printlog("evaluating read %d with %d new and %d old nodes" % (to_add,largest_new,to_add_old))
rnv=reads_tg.get_nodeview(tgs.read_ends[to_add][0])
unhappy=[]
while True:
if rnv.node_id() in tgs.sort_graph() or -rnv.node_id() in tgs.sort_graph():
h=SDG.assess_node_happiness(rnv.node_id(), tgs.sort_graph(), reads_tg)
if h[1]>max_unhappy: unhappy.append(rnv.node_id())
rnvn=[x for x in rnv.next() if x.support().id==to_add]
if not rnvn:break
rnv=rnvn[0].node()
if len(unhappy)>max_unhappy_nodes_in_new_read:
print("skipping rid with unordered nodes")
rids.remove(to_add)
else:
printlog("adding read")
tgs.add_read(to_add,min_support)
write_connected(tgs.dg,"step%d_c-read_added%d.gfa"%(stepnumber,added))
for u in unhappy:
tgs.pop_node(u,to_add)
added+=1
write_connected(tgs.dg,"step%d_c-read_added%d_popped.gfa"%(stepnumber,added))
break
else:
print("can't find any read to add!!!")
break
if added==0:
print("no read was added")
break
# find a read to add that makes the order grow and has few unhappy nodes
printlog("resorting order")
s=sorted(tgs.sort_graph().items(),key=lambda x:x[1])
printlog("order is now",s[-1][1],"bp long")
# if no read found: STOP
# remove unhappy nodes from read until all read is happy
# add read to order
```
```python=
local_happiness_assembler(tgsu,rtgsu,trg_ntu,11676965,max_steps=30)
```
TODO: before adding a read, measure its' node's future happiness, remove future unhappy nodes, but only use reads that have few of them
TODO: Choose reads not by "nodes shared/new" but by growing the order in one direction.
Future hapiness: a node with good connections fw/bw to the order (in the raw graph)
Keep a "nodes used per read" count for all reads in the TGS. Future happy nodes appear in reads that have a high "nodes used per read" count.
In TheGreedySorter:
- NodesUsedPerRead: Keep a uint64_t vector, one entry per read, count how many nodes of that read are included in the order.
- If a node "belongs" in the order, all the reads mapping to the node in the raw graph will have HIGH NodesUsedPerRead values.
- If a node "belongs" in a parallel order, all the reads mapping to the node will have low or null NodesUsedPerRead values.
- Nodes towards the end of the order will have some reads that span in the direction the order has not yet expanded to, so their values will always be low. We need to judge a node belonging only on reads that go towards the order (or, allow for some of the reads to have a small value, which would not hurt anyway).
To add a read to the order, pick one where all its nodes have a high belonging value.
```
c=Counter()
for nv in tgsu.dg.get_all_nodeviews(include_disconnected=False):
nv=trg_ntu.get_nodeview(nv.node_id())
c.update(list(set([x.support().id for x in nv.next()+nv.prev()])))
```
now the counter has a list of reads and how many nodes they share with the order.
if node should go in the order, all its reads should have a high shared value with the order.
node future happines: how many reads the node shares with other nodes in the haplotype (ish).
read future happiness: sum of all its node's happiness.
you want to add reads that have a certain amount of shared nodes (i.e. 10/20), and really high read future happiness.
```python=
def reverse_order(dg,order):
end=max([x[1]+dg.get_nodeview(x[0]).size() for x in order.items()])
reverse_order={}
for x in order.items():
reverse_order[-x[0]]=end-x[1]-dg.get_nodeview(x[0]).size()
return reverse_order
def score_read_to_add(tgs,rid,reverse=False,):
"""gives a score to add the read fw, or bw if reversed"""
order=tgs.sort_graph()
if reverse:
order=reverse_order(trg_ntu,order)
rtn=tgs.thread_nodes(rid)
#first figure out what the relative position of the read is
fw=bw=0
for nid in rtn:
if nid in order: fw+=1
elif -nid in order: bw+=1
#print("read",rid,"shared nodes with order:",fw,"fw,",bw,"bw")
if bw>fw:
#print("reversing read")
rtn=[-x for x in rtn[::-1]]
rid=-rid
bw,fw=fw,bw
if fw<10 or bw>2:
#print("read does not have 10+ nodes FW and <3 BW, scoring None")
return None
tnodes_in_order=[x for x in enumerate(rtn) if x[1] in order]
last_order_nodes=[x[0] for x in sorted(order.items(),key=lambda x:x[1])[-30:]]
if tnodes_in_order[-1][1] not in last_order_nodes: return None
#TODO: validate the order and such?
#nodes in order, position in read of last_node_in_order, nodes after
return (len(tnodes_in_order),tnodes_in_order[-1][0]+1,len(rtn)-tnodes_in_order[-1][0])
def find_reads_to_add(tgs,reverse=False,last_nodes=30, min_hit_perc=.30, min_new=20):
rids=set()
order=sorted(tgs.sort_graph().items(),key=lambda x:x[1])
if reverse:
nodes=[x[0] for x in order[:last_nodes]]
else:
nodes=[x[0] for x in order[-last_nodes:]]
for nid in nodes:
for rid in tgs.rids_from_node(trg_ntu.get_nodeview(nid)):
rids.add(rid)
rta=[]
for rid in rids:
s=score_read_to_add(tgs,rid,reverse)
if s is None: continue
if s[0]/s[1]>=min_hit_perc and s[2]>=min_new:
rta.append((rid,s))
#print(rid,s,'%.2f'%(100.0*s[0]/s[1]))
#else: print(rid,s)
return sorted(rta,key=lambda x:-x[1][2])
from datetime import datetime
def write_connected(dg,filename):
dg.write_to_gfa1(filename,selected_nodes=[x.node_id() for x in dg.get_all_nodeviews(include_disconnected=False)])
def printlog(*args, **kwargs):
print(datetime.now().strftime('%H:%M:%S'),*args, **kwargs)
def local_happiness_assembler(tgs, rtgs, reads_tg, starting_rid,
min_support=2, max_unhappy=0, max_disconnected_perc=.4, min_shared=10, min_new=10, max_unhappy_nodes_in_new_read=10, max_steps=10000000):
tgs.start_from_read(starting_rid,min_support)
for stepnumber in range(max_steps):
printlog("step #",stepnumber,"starting")
write_connected(tgs.dg,"step%d_a-starting.gfa"%(stepnumber))
# remove unhappy nodes from order
printlog("sorting graph")
order=tgs.sort_graph()
printlog("full order hapiness assessment")
nbh=[(x[0],x[1],SDG.assess_node_happiness(x[0],order,reads_tg)) for x in sorted(order.items(),key=lambda x:x[1])] # (nid, pos, (happy, unhappy, disconnected))
#first remove all nodes that are too disconnected: they won't get any better
printlog("removing disconnected")
unhappy_nodes=[]
removed_disconnected=removed_unhappy=0
for x in nbh:
if x[2][0]==0 or 1.0*x[2][2]/sum(x[2])>max_disconnected_perc:
#pop_node_from_all(tgs.dg,)
#print(x, "should go!")
tgs.remove_node(x[0])
removed_disconnected+=1
elif x[2][1]>max_unhappy:
unhappy_nodes.append(x)
printlog("removing unhappy")
#now remove any unhappy nodes
unhappy_nodes.sort(key=lambda x:x[2][1])
for x in unhappy_nodes:
h=SDG.assess_node_happiness(x[0],tgs.sort_graph(),reads_tg)
if h[1]>max_unhappy:
tgs.remove_node(x[0])
removed_unhappy+=1
printlog("%d disconnected and %d unhappy nodes removed from order"%(removed_disconnected,removed_unhappy))
write_connected(tgs.dg,"step%d_b-cleaned.gfa"%(stepnumber))
added=0
for rids in [find_reads_to_add(tgs),find_reads_to_add(tgs,True)]:
while rids:
to_add=rids[0][0]
if to_add in tgs.used_reads:
print("skipping already-used read")
rids=rids[1:]
continue
printlog("evaluating read",rids[0])
rnv=reads_tg.get_nodeview(tgs.read_ends[to_add][0])
unhappy=[]
while True:
if rnv.node_id() in tgs.sort_graph() or -rnv.node_id() in tgs.sort_graph():
h=SDG.assess_node_happiness(rnv.node_id(), tgs.sort_graph(), reads_tg)
if h[1]>max_unhappy: unhappy.append(rnv.node_id())
rnvn=[x for x in rnv.next() if x.support().id==to_add]
if not rnvn:break
rnv=rnvn[0].node()
if len(unhappy)>max_unhappy_nodes_in_new_read:
print("skipping rid with unordered nodes")
rids=rids[1:]
else:
printlog("adding read")
tgs.add_read(to_add,min_support)
write_connected(tgs.dg,"step%d_c-read_added%d.gfa"%(stepnumber,added))
for u in unhappy:
tgs.pop_node(u,to_add)
write_connected(tgs.dg,"step%d_c-read_added%d_popped.gfa"%(stepnumber,added))
added+=1
break
if added==0:
print("no read was added")
break
# find a read to add that makes the order grow and has few unhappy nodes
printlog("resorting order")
s=sorted(tgs.sort_graph().items(),key=lambda x:x[1])
printlog("order is now",s[-1][1],"bp long")
# if no read found: STOP
# remove unhappy nodes from read until all read is happy
# add read to order
```
Strategy #1: Partial orders from Longer reads
1 - Take all reads that have more than x nodes.
2 - Filter their nodes to a curated list of "happy nodes in the read".
use this in a single read:
```python=
# def sort_a_read(tgs,rid,reads_tg,max_unhappy=1,max_disconnected_perc=.4):
# tgs.start_from_read(rid,0)
# printlog("sorting graph")
# order=tgs.sort_graph()
# printlog("full order hapiness assessment")
# nbh=[(x[0],x[1],SDG.assess_node_happiness(x[0],order,reads_tg)) for x in sorted(order.items(),key=lambda x:x[1])] # (nid, pos, (happy, unhappy, disconnected))
# #first remove all nodes that are too disconnected: they won't get any better
# printlog("removing disconnected")
# unhappy_nodes=[]
# removed_disconnected=removed_unhappy=0
# for x in nbh:
# if x[2][0]==0 or 1.0*x[2][2]/sum(x[2])>max_disconnected_perc:
# #pop_node_from_all(tgs.dg,)
# #print(x, "should go!")
# tgs.remove_node(x[0])
# removed_disconnected+=1
# elif x[2][1]>max_unhappy:
# unhappy_nodes.append(x)
# printlog("removing unhappy")
# #now remove any unhappy nodes
# unhappy_nodes.sort(key=lambda x:x[2][1])
# for x in unhappy_nodes:
# h=SDG.assess_node_happiness(x[0],tgs.sort_graph(),reads_tg)
# if h[1]>max_unhappy:
# tgs.remove_node(x[0])
# removed_unhappy+=1
# printlog("%d disconnected and %d unhappy nodes removed from order"%(removed_disconnected,removed_unhappy))
# return tgs.sort_graph()
class HappySorter():
def __init__(self,sorter,all_threads_graph):
self.sorted_graph=SDG.DistanceGraph(all_threads_graph.sdg)
self.used_reads=set()
self.sorter=sorter
self.all_threads_graph=all_threads_graph
def thread_size(self,rid):
ends=self.sorter.get_thread_ends(rid)
e0ix=[x.support().index for x in self.all_threads_graph.get_nodeview(ends[0]).next() if x.support().id==rid][0]
e1ix=[x.support().index for x in self.all_threads_graph.get_nodeview(ends[1]).next() if x.support().id==rid][0]
#print(e0ix,e1ix)
return max([e0ix,e1ix])+1
# def sort_a_read(self,rid,max_unhappy=1,max_disconnected_perc=.4):
# self.sorter.start_from_read(rid,0)
# printlog("sorting graph")
# order=self.sorter.sort_graph()
# printlog("full order hapiness assessment")
# nbh=[(x[0],x[1],SDG.assess_node_happiness(x[0],order,self.all_threads_graph)) for x in sorted(order.items(),key=lambda x:x[1])] # (nid, pos, (happy, unhappy, disconnected))
# #first remove all nodes that are too disconnected: they won't get any better
# printlog("removing disconnected")
# unhappy_nodes=[]
# removed_disconnected=removed_unhappy=0
# for x in nbh:
# if x[2][0]==0 or 1.0*x[2][2]/sum(x[2])>max_disconnected_perc:
# #pop_node_from_all(tgs.dg,)
# #print(x, "should go!")
# self.sorter.remove_node(x[0])
# removed_disconnected+=1
# elif x[2][1]>max_unhappy:
# unhappy_nodes.append(x)
# printlog("removing unhappy")
# #now remove any unhappy nodes
# unhappy_nodes.sort(key=lambda x:x[2][1])
# for x in unhappy_nodes:
# h=SDG.assess_node_happiness(x[0],self.sorter.sort_graph(), self.all_threads_graph)
# if h[1]>max_unhappy:
# self.sorter.remove_node(x[0])
# removed_unhappy+=1
# printlog("%d disconnected and %d unhappy nodes removed from order"%(removed_disconnected,removed_unhappy))
# return self.sorter.sort_graph()
def start_from_read(self,rid):
self.sorted_graph=SDG.DistanceGraph(self.all_threads_graph.sdg)
self.used_reads=set()
self.add_read(rid)
#def add_read(self,rid):
# self.used_reads.add(rid)
# for x in sorted(self.sort_a_read(rid).items(),key=lambda x:x[1])[:-1]:
# nid=x[0]
# l=self.sorter.dg.get_nodeview(nid).next()[0]
# self.sorted_graph.add_link(-nid,l.node().node_id(),l.distance(),l.support())
def add_read(self,rid):
self.used_reads.add(rid)
hr=SDG.make_thread_happy(lrr.read_threads[abs(rid)],self.all_threads_graph)
for i in range(len(hr)-1):
self.sorted_graph.add_link(-hr[i].node,hr[i+1].node,hr[i+1].start-hr[i].end,SDG.Support(id=rid))
def find_longest_threads_per_node(self,include_used=False,top=2):
c=Counter()
for nv1 in self.sorted_graph.get_all_nodeviews(False,False):
nv=self.all_threads_graph.get_nodeview(nv1.node_id())
ml=0
mlrid=0
for l in nv.prev()+nv.next():
if include_used or l.support().id not in self.used_reads:
if len(lrr.read_threads[l.support().id])>ml:
mlrid=l.support().id
ml=len(lrr.read_threads[l.support().id])
if ml: c[mlrid]+=1
return c
def sorted_order(self):
cc=self.sorted_graph.get_connected_component( self.sorted_graph.get_all_nodeviews(include_disconnected=False)[0].node_id(), signed_nodes=True)
#print(cc)
order = sorted(SDG.sort_cc(self.sorted_graph,cc).items(),key=lambda x:x[1])
#print(order)
if len(self.sorted_graph.get_all_nodeviews(include_disconnected=False))!=len(order):
print("WARNING: it appears there's more than one component in HappySorter's sorted_graph")
return order
def pick_for_extension(self,happy_in_order=4,end_size=10,fw=True):
"""picks a read fw/bw for extension, these are the furthest-reaching reads that have at least X happy nodes among the last Y nodes of the current order"""
order=self.sorted_order()
if fw: nids=[x[0] for x in order[-end_size:]]
else: nids=[x[0] for x in order[:end_size]][::-1]
rnids=[-x for x in nids]
#print(len(nids))
#print(nids)
c=Counter()
last_link={}
for nid in nids:
for l in self.all_threads_graph.get_nodeview(nid).next():
tid=l.support().id
if tid in self.used_reads: continue
tidx=l.support().index
if tid not in last_link or last_link[tid]<tidx: last_link[tid]=tidx
c[tid]+=1
#print(c.most_common(10))
c2=Counter()
for rid,oldc in c.items():
if oldc<happy_in_order: continue
hn=[x.node for x in SDG.make_thread_happy(lrr.read_threads[abs(rid)],self.all_threads_graph)]
c2[rid]= max([ len(set(hn).intersection(set(nids))), len(set(hn).intersection(set(rnids)))])
fts=[(x,self.thread_size(x)-last_link[x]) for x in c2.keys() if c2[x]>=happy_in_order]
return sorted(fts,key=lambda x:x[1],reverse=True)
#for x in :
# print(x)
# for nv1 in self.sorted_graph.get_all_nodeviews(False,False):
# nv=self.all_threads_graph.get_nodeview(nv1.node_id())
# ml=0
# mlrid=0
def write_connected(self,filename):
self.sorted_graph.write_to_gfa1(filename,selected_nodes=[x.node_id() for x in self.sorted_graph.get_all_nodeviews(include_disconnected=False)])
def extension_cycles(self,reads_per_end=1,cycles=10,happy_in_order=4,end_size=10):
print("Running %d extension cycles with the farthest reaching %d reads"%(cycles,reads_per_end))
order=self.sorted_order()
print("starting sorted length: %dbp"%(order[-1][1]-order[0][1]))
for c in range(cycles):
print("\n\nCycle", c)
pef=self.pick_for_extension(happy_in_order=happy_in_order,end_size=end_size,fw=True)
for x in range(reads_per_end):
if len(pef)>x:
print("read #%d fw: %d"%(x,pef[x][0]))
self.add_read(pef[x][0])
peb=self.pick_for_extension(happy_in_order=happy_in_order,end_size=end_size,fw=False)
for x in range(reads_per_end):
if len(peb)>x:
print("read #%d bw: %d"%(x,peb[x][0]))
self.add_read(peb[x][0])
order=self.sorted_order()
print("Sorted length: %dbp"%(order[-1][1]-order[0][1]))
def pick_for_thickening():
def remove_read():
```
Read nodes (rtn): [A, B, C, D, F, G, I, K]
Order nodes: [P, Q, R, S, T, A, B, D, E, F]
tnodes_in_order= [(0, A), (1, B), (3, D), (4, F)]
return (4, 4, 4)
```python=
hs=HappySorter(tgsu,trg_ntu);hs.start_from_read(491667)
hs.extension_cycles(2,6)
hs.add_read(3785826)
hs.add_read(13253989)
hs.write_connected("hsbefore.gfa")
hs.add_read(14272921)
hs.write_connected("hsafter.gfa")
```
```python=
#lrr2 will get modified and will have only threads with their happy nodes
lrr2.simple_thread_reads()
SDG.make_all_threads_happy(lrr2,trg_ntu)
htrg=lrr.dg_from_threads(False,True,10)#this is largely irrelevant for the happysorter, but useful
hs=HappySorter(htgs,trg_ntu);hs.start_from_read(11676965)
hs.extension_cycles(10,50,4,20)
```
```python=
def all_nodes_and_distances_across_rid(nv,rid,max_distance=10000000000):
nodes=[]
d=-nv.size()
while True:
try:
lf=[l for l in nv.next() if l.support().id==rid][0]
except:
break
d+=nv.size()+lf.distance()
if d>max_distance: break
nodes.append((d,lf.node().node_id()))
nv=lf.node()
return nodes
def get_all_nodes_up_to_dist(fnv,dist):
all_nodes=Counter()
all_rids=set()
for x in fnv.next()+fnv.prev():
all_rids.add(x.support().id)
for rid in all_rids:
for ndist,nid in all_nodes_and_distances_across_rid(fnv,rid,dist):
all_nodes[nid]+=1
for ndist,nid in all_nodes_and_distances_across_rid(fnv.rc(),rid,dist):
all_nodes[-nid]+=1
return all_nodes
def nodes_by_dist_and_support(fnv,dist,links):
nodes=[]
for nid,count in get_all_nodes_up_to_dist(fnv,dist).items():
if count>=links: nodes.append(nid)
return nodes
```
```python=
def rtg_to_tube(rtg):
tube={"nodes":[],"tracks":[]}
tids=set()
for xx in sorted(SDG.sort_cc( rtg,rtg.get_connected_component(rtg.get_all_nodeviews(include_disconnected=False)[0].node_id(),signed_nodes=True) ).items(),key=lambda x:x[1]):
x=rtg.get_nodeview(abs(xx[0]))
tube["nodes"].append( dict(name = "%d"%x.node_id(), width = 10, seq = "seq%d"%x.node_id() ) )
for l in x.next()+x.prev():
tids.add(l.support().id)
for tnum,tid in enumerate(tids):
tube["tracks"].append(dict(
id="%d"%tnum,
name="%d"%tid,
indexOfFirstBase=1,
sequence=[str(abs(np.node)) for np in rtg.get_thread(tid)],
freq=1
))
return tube
```
# Cleaning a raw graph
## Removing nodes from threads where they OBVIOUSLY do not belong
```python=
from collections import Counter
def score_threads(rtg,nid,support=1):
n1tn=rtg.node_thread_neighbours(nid)
for tid in rtg.node_threads(nid):
c=Counter([n1tn[abs(np.node)]>support for np in rtg.get_thread(tid)])
c[True]-=1
print("%10d %3d %3d %.2f%%"%(tid,c[True],c[False],c[True]*100.0/(c[True]+c[False])))
def clean_threads(rtg,nid,support=1):
n1tn=rtg.node_thread_neighbours(nid)
to_pop=[]
for tid in rtg.node_threads(nid):
c=Counter([n1tn[abs(np.node)]>support for np in rtg.get_thread(tid)])
c[True]-=1
if c[True]<3: to_pop.append(tid)
for tid in to_pop:
rtg.pop_node(nid,tid)
return len(to_pop)
```
# from here on
```python=
def sort_rtg(rtg):
cc=rtg.get_connected_component( rtg.get_all_nodeviews(include_disconnected=False)[0].node_id(), signed_nodes=True)
#print(cc)
order = sorted(SDG.sort_cc(rtg,cc).items(),key=lambda x:x[1])
#print(order)
if len(order)==0: print("WARNING: graph can't be sorted")
elif len(rtg.get_all_nodeviews(include_disconnected=False))!=len(order):
print("WARNING: there's more than one component the graph")
return []
return order
def thread_internal_linkage(rtg,tid):
nids=[x.node for x in rtg.get_thread(tid)]
l=[]
for n1 in nids:
l.append([])
for n2 in nids:
l[-1].append(len(set(rtg.node_threads(n1))&set(rtg.node_threads(n2))))
return l
def thread_linkage_support(rtg,tid):
t=rtg.get_thread(tid)
ot={}
for i in range(len(t)):
for nt in rtg.node_threads(t[i].node):
if nt==tid: continue
if nt not in ot: ot[nt]=[]
ot[nt].append(i)
covd=[0 for x in t]
for x in ot.values():
if len(x)>1:
covd[x[0]+1]+=1
covd[x[-1]]-=1
cc=0
tls=[]
for x in covd:
cc+=x
tls.append(cc)
return tls
def is_chimeric(rtg,tid,min_support=1):
if [x for x in thread_linkage_support(rtg,tid)[1:-1] if x<min_support]:
return True
return False
def add_threads_from_node(rtgdest,rtgsource,node_id, discard_chimeric=True):
print("adding threads from node %d"%node_id)
tids=rtgsource.node_threads(node_id)
if not tids:
print("node has no threads in source rtg")
return
for tid in tids:
if rtgdest.get_thread(tid):
print(" x skipping thread %d already present in dest rtg"%tid)
elif discard_chimeric and is_chimeric(rtgsource,tid):
print(" x skipping thread %d flagged as chimeric"%tid)
else:
print(" v adding thread %d"%tid)
rtgdest.add_thread(tid, rtgsource.get_thread(tid), False, 2)
```
```python=
rpl=rtg.clean_repeat_nodes_popping_list()
len(rpl) # --> 7318860
rtg.apply_popping_list(rpl)
pl=rtg.clean_all_nodes_popping_list()
len(pl) # --> 14176414 (tiene que ser mas!)
rtg.apply_popping_list(pl)
rtg.dump("rtg_rep_and_clean.dump")
lrtg=SDG.ReadThreadsGraph(ws.sdg)
add_threads_from_node(lrtg,rtg,1)
sort_rtg(lrtg) #-> check orders for other nodes ordered here.
```
```python=
from collections import Counter
class RTGSorter():
def __init__(self,all_rtg,safe_node_cv=5,thread_anchor_nodes=5,last_safe_nodes=10):
self.safe_node_cv=safe_node_cv
self.thread_anchor_nodes=thread_anchor_nodes
self.last_safe_nodes=last_safe_nodes
self.all_rtg=all_rtg
self.rtg=SDG.ReadThreadsGraph(self.all_rtg.sdg)
self.used_threads=set()
def find_safe_nodes(self):
safe_nodes=[nv.node_id() for nv in self.rtg.get_all_nodeviews(include_disconnected=False)
if len(self.rtg.node_threads(nv.node_id()))>=self.safe_node_cv]
return safe_nodes
def find_fw_extension_threads(self,rev=False,tan=None,last_safe_nodes=None):
if last_safe_nodes is None: last_safe_nodes=self.last_safe_nodes
if tan is None: tan=self.thread_anchor_nodes
sn=self.find_safe_nodes()
safesort=self.sort_safe_rtg()
#stupidly re-make this into a dict
ssd={}
if rev:
lsn=[x[0] for x in safesort[:last_safe_nodes]]
else:
lsn=[x[0] for x in safesort[-last_safe_nodes:]]
for x in safesort:
ssd[x[0]]=x[1]
c=Counter()
for x in sn:
c.update([x for x in self.all_rtg.node_threads(x) if x not in self.used_threads])
ct=[x[0] for x in c.items() if x[1]>=tan]
st=[]
for tid in ct:
thread=[x.node for x in self.all_rtg.get_thread(tid)]
#print("\nanalysing candidate thread",tid," length=",len(thread))
if rev: thread=[-x for x in thread[::-1]]
f=[]
b=[]
for n in enumerate(thread):
if n[1] in ssd:
f.append(n)
if -n[1] in ssd:
b.append(n)
#print("B:",b)
#print("F:",f)
if len(f)==0:
st.append((tid,b[0][0], len([1 for x in thread if len(self.rtg.node_threads(x))==tan-1]), len([1 for x in thread if x in lsn or -x in lsn])))
if len(b)==0:
st.append((tid,len(thread)-f[-1][0]+1, len([1 for x in thread if len(self.rtg.node_threads(x))==tan-1]), len([1 for x in thread if x in lsn or -x in lsn])))
return sorted(st,key=lambda x:(-x[3],-x[2],-x[1]))
def thread_linkage_support(self,rtg,tid):
t=rtg.get_thread(tid)
ot={}
for i in range(len(t)):
for nt in rtg.node_threads(t[i].node):
if nt==tid: continue
if nt not in ot: ot[nt]=[]
ot[nt].append(i)
covd=[0 for x in t]
for x in ot.values():
if len(x)>1:
covd[x[0]+1]+=1
covd[x[-1]]-=1
cc=0
tls=[]
for x in covd:
cc+=x
tls.append(cc)
return tls
def is_chimeric(self,tid,min_support=1):
if [x for x in self.thread_linkage_support(self.all_rtg,tid)[1:-1] if x<min_support]:
return True
return False
def add_thread(self,tid,discard_chimeric=True):
if tid in self.used_threads or self.rtg.get_thread(tid):
print(" x skipping thread %d already present in dest rtg"%tid)
return False
elif discard_chimeric and self.is_chimeric(tid):
print(" x skipping thread %d flagged as chimeric"%tid)
return False
else:
print(" v adding thread %d"%tid)
self.rtg.add_thread(tid, self.all_rtg.get_thread(tid), False, 2)
self.used_threads.add(tid)
return True
def remove_thread(self,tid):
self.used_threads.remove(tid)
self.rtg.remove_thread(tid)
def sort_rtg(self,rtg=None):
if rtg is None: rtg=self.rtg
all_connected=rtg.get_all_nodeviews(include_disconnected=False)
cc=rtg.get_connected_component( all_connected[0].node_id(), signed_nodes=True)
order = sorted(SDG.sort_cc(rtg,cc).items(),key=lambda x:x[1])
if len(order)==0: print("WARNING: graph can't be sorted")
elif len(all_connected)!=len(order):
print("WARNING: there's more than one component the graph")
return []
return order
def sort_safe_rtg(self,safe_cv=None):
if safe_cv is None: safe_cv=self.safe_node_cv
srtg=SDG.ReadThreadsGraph(self.rtg.sdg)
for tid in self.used_threads:
srtg.add_thread(tid, self.rtg.get_thread(tid), False, 2)
to_disconnect=[nv.node_id() for nv in srtg.get_all_nodeviews(include_disconnected=False) if len(srtg.node_threads(nv.node_id()))<safe_cv]
for nid in to_disconnect:
srtg.pop_node_from_all(nid)
return self.sort_rtg(srtg)
def auto_extend(self,steps=10000):
chimeric_reads=set()
for step in range(steps):
print("\nAuto extension step",step)
ssort=self.sort_safe_rtg()
print("Starting with a safe graph from %d to %d, total distance %d bp" % (ssort[0][0],ssort[-1][0],ssort[-1][1]-ssort[0][1]))
extended=0
for fc in [ [x[0] for x in self.find_fw_extension_threads() if x not in chimeric_reads], [x[0] for x in self.find_fw_extension_threads(rev=True) if x not in chimeric_reads]]:
for tid in fc:
if tid in chimeric_reads: continue
if not self.add_thread(tid):
chimeric_reads.add(tid)
continue
ssort=self.sort_safe_rtg()
if ssort==[]:
print("last thread ruined safe order! marking as chimeric and backtracking")
chimeric_reads.add(tid)
self.remove_thread(tid)
continue
print("Extended safe graph from %d to %d, total distance %d bp" % (ssort[0][0],ssort[-1][0],ssort[-1][1]-ssort[0][1]))
extended+=1
break
if extended==0:
print("can't find a thread to add on any end, quitting")
```
```python=
from collections import Counter
class RTGThreadsCluster:
def __init__(self,all_rtg):
self.all_rtg=all_rtg
self.rtg=SDG.ReadThreadsGraph(self.all_rtg.sdg)
self.threads=set()
self.used_threads=set()
self.chimeric_threads=set()
def thread_linkage_support(self,rtg,tid):
t=rtg.get_thread(tid)
ot={}
for i in range(len(t)):
for nt in rtg.node_threads(t[i].node):
if nt==tid: continue
if nt not in ot: ot[nt]=[]
ot[nt].append(i)
covd=[0 for x in t]
for x in ot.values():
if len(x)>1:
covd[x[0]+1]+=1
covd[x[-1]]-=1
cc=0
tls=[]
for x in covd:
cc+=x
tls.append(cc)
return tls
def is_chimeric(self,tid,min_support=1):
if [x for x in self.thread_linkage_support(self.all_rtg,tid)[1:-1] if x<min_support]:
return True
return False
def add_thread(self,tid,discard_chimeric=True):
if tid in self.used_threads or self.rtg.get_thread(tid):
print(" x skipping thread %d already present in dest rtg"%tid)
return False
elif discard_chimeric and self.is_chimeric(tid):
print(" x skipping thread %d flagged as chimeric"%tid)
self.chimeric_threads.add(tid)
return False
else:
print(" v adding thread %d"%tid)
self.rtg.add_thread(tid, self.all_rtg.get_thread(tid), False, 2)
self.used_threads.add(tid)
return True
def remove_thread(self,tid):
self.used_threads.remove(tid)
self.rtg.remove_thread(tid)
def find_threads_to_add(self, safe_cvg=3, min_nodes=5, max_span=20):
print("Collecting safe nodes and candidate threads...")
tc=Counter()
nc=Counter()
for nv in self.rtg.get_all_nodeviews(include_disconnected=False):
nt=self.rtg.node_threads( nv.node_id() )
nc[nv.node_id()]=len(nt)
if len(nt)<safe_cvg: continue
for tid in self.all_rtg.node_threads( nv.node_id() ):
if tid in nt or tid in self.chimeric_threads: continue
tc[tid]+=1
print("Checking threads to add")
tcc=0
to_add=[]
for tid,c in tc.items():
if c<min_nodes: continue
tcc+=1
tpa=[x[0] for x in enumerate(self.all_rtg.get_thread(tid)) if nc[x[1].node]>=safe_cvg]
for i in range(len(tpa)-min_nodes):
if tpa[i+min_nodes]-tpa[i]<=max_span:
to_add.append(tid)
break
print("There was",tcc,"candidates")
print("There was",len(to_add),"candidates with",min_nodes,"in",max_span,"consecutive nodes included")
return to_add
```
```python=
from statistics import median
def make_summary_connection_graph(mdg,links_to_connect=3):
"""returns a new dg with only one connection per node end pair (i.e. mdg to dg)"""
conn={}
dg=SDG.DistanceGraph(mdg.sdg)
for nv in mdg.get_all_nodeviews(include_disconnected=False,both_directions=True):
for ln in nv.next():
ck=(-nv.node_id(),ln.node().node_id())
if ck[0]<=ck[1]:
if ck not in conn: conn[ck]=[]
conn[ck].append(ln.distance())
for c,d in conn.items():
if len(d)>=links_to_connect:
dg.add_link(c[0],c[1],int(median(d)))
return dg
```
Check if a node's thread can be clustered in independent haplotypes (i.e. false anchor!)
```python=
#use this on every node, if the clusters many and complicated... pop this shit!
def cluster_threads_in_node(rtg,nid,min_shared=4):
tnodes={}
for tid in rtg.node_threads(nid):
tnodes[tid]=set([abs(x.node) for x in rtg.get_thread(tid) if abs(x.node)!=abs(nid)])
tconns={}
for tid1,nodes1 in tnodes.items():
tconns[tid1]=[]
for tid2,nodes2 in tnodes.items():
if len(nodes1&nodes2)<min_shared: continue
tconns[tid1].append(tid2)
used_tids=set()
clusters=[]
for tid,conns in tconns.items():
if tid in used_tids: continue
clusters.append(set([tid]))
used_tids.add(tid)
to_explore_tids=conns
while to_explore_tids:
ntet=[]
for tid2 in to_explore_tids:
if tid2 in used_tids: continue
clusters[-1].add(tid2)
used_tids.add(tid2)
for tid3 in tconns[tid2]:
if tid3 not in used_tids:
to_explore_tids.append(tid3)
to_explore_tids=ntet
return clusters
```
# Computing multi-clustering nodes fast
- Create all sets of nodes for all threads
- For every node:
- Cluster threads by node sets, ordered by size
- Compare size of first and second clusters:
- one significant cluster: remove node from threads not in the significant cluster.
- two or more significant clusters: remove node from all threads.
```python=
def threads_to_nodesets(rtg):
nodesets={}
for nv in rtg.get_all_nodeviews(include_disconnected=False):
nid=nv.node_id()
for tid in rtg.node_threads(nid):
if tid not in nodesets:
nodesets[tid]=set()
nodesets[tid].add(nid)
return nodesets
def cluster_threads_in_node(rtg,tnodes,nid,min_shared=4):
tconns={}
ntids=rtg.node_threads(nid)
for tid1 in ntids:
nodes1=tnodes[tid1]
tconns[tid1]=[]
for tid2 in ntids:
nodes2=tnodes[tid2]
if len(nodes1&nodes2)<min_shared: continue
tconns[tid1].append(tid2)
used_tids=set()
clusters=[]
for tid,conns in tconns.items():
if tid in used_tids: continue
clusters.append(set([tid]))
used_tids.add(tid)
to_explore_tids=conns
while to_explore_tids:
ntet=[]
for tid2 in to_explore_tids:
if tid2 in used_tids: continue
clusters[-1].add(tid2)
used_tids.add(tid2)
for tid3 in tconns[tid2]:
if tid3 not in used_tids:
to_explore_tids.append(tid3)
to_explore_tids=ntet
return clusters
```
```python=
rtg #this graph has all the threads
single_node_rtg=SDG.ReadThreadsGraph(ws.sdg)
for tid in trg.node_threads(nid):
single_node_rtg.add_thread(tid, rtg.get_thread(tid))
single_node_rtg.write_to_gfa1("single_node.gfa",selected_nodes=[x.node_id() for x in single_node_rtg.get_all_nodeviews(include_disconnected=False)])
```
```python=
def sift_candidates(
```
Sign in with Wallet
Connect another wallet