Skip to content

Commit 3326156

Browse files
aafshinfardaafshinfard
authored
sqcosbin and k3bin (#146)
* added sqcosbin and k3bin * changes to include distributed+sqcosbin (ext) * happy pylint * bug fix * added cost
1 parent 4458a4c commit 3326156

2 files changed

Lines changed: 66 additions & 6 deletions

File tree

‎bin/physlr-make‎

Lines changed: 20 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -184,10 +184,11 @@ help:
184184
@echo " min_component_size minimum number of barcodes in a backbone [50]."
185185
@echo " minimum_barcode_multiplicity minimum number of minimizers per barcode [10]."
186186
@echo " maximum_barcode_multiplicity maximum number of minimizers per barcode [5000]."
187-
@echo " mol_strategy molecule separation strategy [distributed]. Available options are bc, bc+k3, distributed."
187+
@echo " mol_strategy molecule separation strategy [distributed]. Available options are bc, bc+k3, distributed, ext."
188188
@echo " bc (biconnected componenets) is the least conservative and is only suitable for datasets with low barcode multiplicity."
189-
@echo " bc+k3 (biconnected componenets + k-3 cliques) is most conservative and requires more time."
190-
@echo " distributed is a modified version of bc+k3 that is faster but less conservative."
189+
@echo " bc+k3 (biconnected componenets + k-3 cliques) is more conservative than bc and requires more time."
190+
@echo " distributed is a modified version of bc+k3 that is faster than bc+k3 but may be more (or even less) conservative."
191+
@echo " ext (extensive) mixes distributed with a modified version of sqcos (cosine similarity of squared adjacency matrix) which makes it more conservative."
191192
@echo " bloom_filter_size size of bloom filter [10000000000] (10G)."
192193
@echo " arcs Use ARCS to augment scaffolds (only compatible with ARCS v1.1.1) [false]."
193194
@echo ""
@@ -989,6 +990,10 @@ min_path_size=200
989990
%.overlap.m$m.k3.mol.backbone.map.$(draft).n10.bed: %.overlap.m$m.k3.mol.backbone.path %.tsv $(draft).k$k-w$w.physlr.tsv
990991
$(time) $(python) $(bin)/physlr map -V$V -n10 $^ >$@
991992

993+
# Map the draft assembly to the backbone graph and output BED (K3-cliques and sqcosbin).
994+
%.overlap.m$m.ext.mol.backbone.map.$(draft).n10.bed: %.overlap.m$m.ext.mol.backbone.path %.tsv $(draft).k$k-w$w.physlr.tsv
995+
$(time) $(python) $(bin)/physlr map -V$V -n10 $^ >$@
996+
992997
# Map the draft assembly to the Louvain backbone graph and output BED.
993998
%.overlap.m$m.louvain.mol.backbone.map.$(draft).n10.bed: %.overlap.m$m.louvain.mol.backbone.path %.tsv $(draft).k$k-w$w.physlr.tsv
994999
$(time) $(python) $(bin)/physlr map -V$V -n10 $^ >$@
@@ -1001,6 +1006,10 @@ min_path_size=200
10011006
%.overlap.m$m.k3.mol.backbone.map.$(ref).n10.bed: %.overlap.m$m.k3.mol.backbone.path %.tsv $(name)/$(ref).k$k-w$w.physlr.tsv
10021007
$(time) $(python) $(bin)/physlr map -V$V -n10 $^ >$@
10031008

1009+
# Map the reference assembly to the backbone graph and output BED (K3-cliques and sqcosbin).
1010+
%.overlap.m$m.ext.mol.backbone.map.$(ref).n10.bed: %.overlap.m$m.ext.mol.backbone.path %.tsv $(name)/$(ref).k$k-w$w.physlr.tsv
1011+
$(time) $(python) $(bin)/physlr map -V$V -n10 $^ >$@
1012+
10041013
# Map the reference assembly to the Louvain backbone graph and output BED.
10051014
%.overlap.m$m.louvain.mol.backbone.map.$(ref).n10.bed: %.overlap.m$m.louvain.mol.backbone.path %.tsv $(name)/$(ref).k$k-w$w.physlr.tsv
10061015
$(time) $(python) $(bin)/physlr map -V$V -n10 $^ >$@
@@ -1013,6 +1022,10 @@ min_path_size=200
10131022
%.overlap.m$m.k3.mol.backbone.map-mkt.$(draft).n10.bed: %.overlap.m$m.k3.mol.backbone.path %.tsv $(draft).k$k-w$w.physlr.tsv
10141023
$(time) $(python) $(bin)/physlr map-mkt -V$V -n10 $^ >$@
10151024

1025+
# Map the draft assembly to the `K3-cliques and sqcosbin` backbone graph using physlr map-mkt and output BED.
1026+
%.overlap.m$m.ext.mol.backbone.map-mkt.$(draft).n10.bed: %.overlap.m$m.ext.mol.backbone.path %.tsv $(draft).k$k-w$w.physlr.tsv
1027+
$(time) $(python) $(bin)/physlr map-mkt -V$V -n10 $^ >$@
1028+
10161029
# Map the draft assembly to the Louvain backbone graph using physlr map-mkt and output BED.
10171030
%.overlap.m$m.louvain.mol.backbone.map-mkt.$(draft).n10.bed: %.overlap.m$m.louvain.mol.backbone.path %.tsv $(draft).k$k-w$w.physlr.tsv
10181031
$(time) $(python) $(bin)/physlr map-mkt -V$V -n10 $^ >$@
@@ -1104,6 +1117,10 @@ g=$(shell awk '{x += $$2} END{print x}' $(name)/$(ref).fa.fai)
11041117
%.k3.mol.tsv: %.tsv
11051118
$(time) $(python) $(bin)/physlr molecules -V$V -t$t --separation-strategy=bc+k3 $< >$@
11061119

1120+
# Separate barcodes into molecules using K3-clique and sqcosbin community detection.
1121+
%.ext.mol.tsv: %.tsv
1122+
$(time) $(python) $(bin)/physlr molecules -V$V -t$t --separation-strategy=distributed+sqcosbin $< >$@
1123+
11071124
# Separate barcodes into molecules using Louvain community detection.
11081125
%.louvain.mol.tsv: %.tsv
11091126
$(time) $(python) $(bin)/physlr molecules -V$V -t$t --separation-strategy=bc+louvain $< >$@

‎physlr/physlr.py‎

Lines changed: 46 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1468,7 +1468,13 @@ def detect_communities_cosine_of_squared(g, node_set, squaring=True, threshold=0
14681468
import numpy as np
14691469
from sklearn.metrics.pairwise import cosine_similarity
14701470

1471+
if len(node_set) < Physlr.args.skip_small:
1472+
return [set(node_set)]
1473+
if len(node_set) == 1:
1474+
return [set(node_set)]
1475+
14711476
communities = []
1477+
14721478
if len(node_set) > 1:
14731479
adj_array = nx.adjacency_matrix(g.subgraph(node_set)).toarray()
14741480
if squaring:
@@ -1586,6 +1592,19 @@ def determine_molecules(g, u, junctions, strategy):
15861592
for component in communities:
15871593
communities_temp.extend(
15881594
Physlr.detect_communities_k_clique(g, component, k=3))
1595+
elif algorithm == "k3bin":
1596+
for component in communities:
1597+
communities_temp.extend(
1598+
Physlr.merge_communities(
1599+
g, [cluster
1600+
for bin_set in
1601+
Physlr.partition_subgraph_into_bins_randomly(
1602+
component)
1603+
for cluster in
1604+
Physlr.detect_communities_k_clique(g, bin_set, k=3)
1605+
]
1606+
)
1607+
)
15891608
elif algorithm == "k4":
15901609
for component in communities:
15911610
communities_temp.extend(
@@ -1599,6 +1618,20 @@ def determine_molecules(g, u, junctions, strategy):
15991618
for component in communities:
16001619
communities_temp.extend(
16011620
Physlr.detect_communities_cosine_of_squared(g, component))
1621+
elif algorithm == "sqcosbin":
1622+
for component in communities:
1623+
communities_temp.extend(
1624+
Physlr.merge_communities(
1625+
g, [cluster
1626+
for bin_set in
1627+
Physlr.partition_subgraph_into_bins_randomly(
1628+
component)
1629+
for cluster in
1630+
Physlr.detect_communities_cosine_of_squared(
1631+
g, bin_set, squaring=True, threshold=Physlr.args.sqcost)
1632+
]
1633+
)
1634+
)
16021635
elif algorithm == "louvain":
16031636
for component in communities:
16041637
communities_temp.extend(
@@ -1622,7 +1655,8 @@ def determine_molecules_process(u):
16221655

16231656
def physlr_molecules(self):
16241657
"Separate barcodes into molecules."
1625-
alg_white_list = {"bc", "cn2", "cn3", "k3", "k4", "cos", "sqcos", "louvain", "distributed"}
1658+
alg_white_list = {"bc", "cn2", "cn3", "k3", "k3bin", "k4",
1659+
"cos", "sqcos", "sqcosbin", "louvain", "distributed"}
16261660
alg_list = self.args.strategy.split("+")
16271661
if not alg_list:
16281662
sys.exit("Error: physlr molecule: missing parameter --separation-strategy")
@@ -2605,8 +2639,11 @@ def parse_arguments():
26052639
argparser.add_argument(
26062640
"--separation-strategy", action="store", dest="strategy", default="bc+k3",
26072641
help="strategy for barcode to molecule separation [bc+k3]. Use a combination"
2608-
" of bc, k3, cos, sqcos, louvain, and distributed"
2609-
" concatenated plus sign (example:bc+k3+bc)")
2642+
" of bc, k3, k3bin, cos, sqcos, sqcosbin, louvain, and distributed"
2643+
" concatenated via plus sign (example:bc+k3+bc)")
2644+
argparser.add_argument(
2645+
"--separation-skip-small", action="store", dest="skip_small", type=int, default=10,
2646+
help="Skip splitting the barcode if neighborhood subgraph is smaller than this [0]")
26102647
argparser.add_argument(
26112648
"--coef", action="store", dest="coef", type=float, default=1.5,
26122649
help="ignore minimizers that occur in Q3+c*(Q3-Q1) or more barcodes [0]")
@@ -2701,6 +2738,12 @@ def parse_arguments():
27012738
argparser.add_argument(
27022739
"--gap-size", action="store", dest="gap_size", type=int, default=100,
27032740
help="gap size used in scaffolding [100].")
2741+
argparser.add_argument(
2742+
"--cost", action="store", dest="cost", type=float, default=0.5,
2743+
help="threshold for `cos` molecule separation [0.5].")
2744+
argparser.add_argument(
2745+
"--sqcost", action="store", dest="sqcost", type=float, default=0.75,
2746+
help="threshold for `sqcos` molecule separation [0.75].")
27042747
argparser.add_argument(
27052748
"--minimizer-overlap", action="store", dest="minimizer_overlap", type=float, default=0,
27062749
help="Percent of edges to remove [0].")

0 commit comments

Comments
 (0)