multi-way contacts file

[erbon7]

Hi, I was wondering if the multi-way contacts derived from the alignment file are stored in a dedicated file and if yes which one ?
The "alignments.parquet" file contains only filtered pairwise alignments if I understood correctly ?
Thx

[eharr]

@erbon7 sorry for the delayed response. In version 0.3.0 there's *.contacts.parquet file with all pairwise contacts. You can use the read_id column to group these to see all the contacts for a given contcatemer/read. Is this what you meant?

[erbon7]

Thanks Eoghan for the answer. I just have another question. In the preprint there is a graph analysis procedure to select the best alignments. Is the code to do that part of the pore-c distribution ?

[eharr]

Yes - there the 2 functions linked here carry it out:

pore-c/pore_c/analyses/alignments.py

Lines 181 to 230 in a865b7e

	def create_align_graph(aligns, gap_fn):
	# we'll visit the alignments in order of increasing endpoint on the read, need to keep
	# the ids as the index in the original list of aligns for filtering later
	aligns = aligns[["read_start", "read_end", "read_length", "align_score"]].copy().sort_values(["read_end"])
	node_ids = list(aligns.index)
	graph = nx.DiGraph()
	# initialise graph with root and sink node, and one for each alignment
	# edges in the graph will represent transitions from one alignment segment
	# to the next
	graph.add_nodes_from(["ROOT", "SINK"] + node_ids)
	for align in aligns.itertuples():
	align_idx = align.Index
	align_score = align.align_score
	gap_penalty_start = gap_fn(align.read_start)
	graph.add_edge("ROOT", align_idx, weight=gap_penalty_start - align_score)
	gap_penalty_end = gap_fn(int(align.read_length - align.read_end))
	graph.add_edge(align_idx, "SINK", weight=gap_penalty_end)
	# for each pair of aligned segments add an edge
	for idx_a, align_idx_a in enumerate(node_ids[:-1]):
	align_a_end = aligns.at[align_idx_a, "read_end"]
	for align_idx_b in node_ids[idx_a + 1 :]: # noqa: E203 black does this
	align_b_score = aligns.at[align_idx_b, "align_score"]
	align_b_read_start = aligns.at[align_idx_b, "read_start"]
	gap_penalty = gap_fn(abs(int(align_b_read_start) - int(align_a_end)))
	graph.add_edge(align_idx_a, align_idx_b, weight=gap_penalty - align_b_score)
	return graph
	def filter_shortest_path(read_df, aligner="minimap2"):
	aligns = read_df[read_df.pass_filter]
	num_aligns = len(aligns)
	if num_aligns < 2:
	# can't build a graph, so nothing is filtered by this method
	return read_df
	if aligner == "minimap2":
	gap_fn = minimap_gapscore
	elif aligner == "bwa":
	gap_fn = bwa_gapscore
	else:
	raise ValueError(f"Unrecognised aligner: {aligner}")
	graph = create_align_graph(aligns, gap_fn)
	distance, shortest_path = nx.single_source_bellman_ford(graph, "ROOT", "SINK")
	for idx in aligns.index:
	if idx not in shortest_path:
	read_df.at[idx, "pass_filter"] = False
	read_df.at[idx, "filter_reason"] = "not_on_shortest_path"
	return read_df