Skip to content

Commit 8f801d5

Browse files
JannesSPJannesSP
committed
update
1 parent b53159b commit 8f801d5

File tree

2 files changed

+37
-32
lines changed

2 files changed

+37
-32
lines changed

‎src/python/misc/csv_to_ms_heatmap.py‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -185,9 +185,9 @@ def plot_heatmap(ams: pd.DataFrame, output_file: str):
185185
+ [""]
186186
)
187187
dataset_labels = [
188-
r"$H.\ sapiens$", r"$E.\ coli$", "SARS-CoV-2", "IVT", r"IVT $H.\ sapiens$", "m1Y",
189-
r"$H.\ sapiens$", r"$S.\ cerevisiae$", "CEVd", "IVT", "psU",
190-
r"$H.\ sapiens$", "Zymo HMW", r"$S.\ Aureus$", r"$P.\ Anserina$", "5mC", "tool average"
188+
r"$H.\ sapiens$", r"$E.\ coli$", "SARS-CoV-2", "IVT", r"IVT $H.\ sapiens$", r"m$^1$Ψ",
189+
r"$H.\ sapiens$", r"$S.\ cerevisiae$", "CEVd", "IVT", "Ψ",
190+
r"$H.\ sapiens$", "Zymo HMW", r"$S.\ Aureus$", r"$P.\ Anserina$", r"$^5$mC", "tool average"
191191
]
192192

193193
# Set the dataset labels

‎src/python/plot/plotToolSegments.py‎

Lines changed: 34 additions & 29 deletions
Original file line numberDiff line numberDiff line change
@@ -162,7 +162,7 @@ def readBasecalls(file: str, readid: str) -> str:
162162
with AlignmentFile(file, "rb", check_sq=False) as bamfile:
163163
for read in bamfile.fetch(until_eof=True):
164164
if readid == read.query_name:
165-
return read.query_sequence
165+
return read.query_sequence, read.get_tag("pi") if read.has_tag("pi") else ""
166166
return None
167167

168168
def main() -> None:
@@ -177,19 +177,22 @@ def main() -> None:
177177

178178
args = parse()
179179

180-
read = readBasecalls(args.basecalls, args.readid)
180+
readid, signalid = readBasecalls(args.basecalls, args.readid)
181181

182182
tools = {
183183
"Dynamont" : np.array(readDynamont(args.dynamont, args.readid)),
184184
"Uncalled4" : np.array(readUncalled4(args.uncalled4, args.readid)),
185-
"f5c Resquiggle" : np.array(readF5CResquiggle(args.f5cresquiggle, args.readid, read, args.k)),
185+
"f5c Resquiggle" : np.array(readF5CResquiggle(args.f5cresquiggle, args.readid, readid, args.k)),
186186
"f5c Eventalign" : np.array(readF5CEventalign(args.f5ceventalign, splitext(args.f5ceventalign)[0] + '.sum', args.readid)),
187187
# "Tombo" : readTombo(args.tombo, args.readid),
188188
"Dorado" : np.array(readDorado(args.dorado, args.readid)),
189189
}
190190

191191
outpath = args.out
192-
signal = r5.read(args.pod5).getpASignal(args.readid)
192+
try:
193+
signal = r5.read(args.pod5).getpASignal(args.readid)
194+
except:
195+
signal = r5.read(args.pod5).getSignal(signalid)
193196

194197
fig, ax = plt.subplots(nrows = 5, figsize=(110,15), dpi=300)
195198
fig.suptitle("Segmentation of the same read by different tools")
@@ -210,18 +213,18 @@ def main() -> None:
210213
ax[i].vlines([start, end], ymin=min(signal), ymax=max(signal), colors=basecolors[base], linestyles='--', linewidth=0.7)
211214
ax[i].add_patch(Rectangle((start, min(signal)), end - start, max(signal)-min(signal), alpha=0.4, edgecolor=basecolors[base], facecolor=basecolors[base]))
212215

213-
if tool == "Dynamont":
214-
prob_ax = ax[i].twinx()
215-
borderProbs, acgt_probs = getDynamontProbs()
216-
# prob_ax.plot(np.arange(len(borderProbs)), borderProbs, c='red', linewidth=0.5, linestyle='-')
217-
prob_ax.plot(acgt_probs[:, 0], c=basecolors['A'], linewidth=1, linestyle='-', label="A probability")
218-
prob_ax.plot(acgt_probs[:, 1], c=basecolors['C'], linewidth=1, linestyle='-', label="C probability")
219-
prob_ax.plot(acgt_probs[:, 2], c=basecolors['G'], linewidth=1, linestyle='-', label="G probability")
220-
prob_ax.plot(acgt_probs[:, 3], c=basecolors['T'], linewidth=1, linestyle='-', label="T probability")
221-
prob_ax.set_ylabel("Nucleotide Probabilities", fontsize=10)
222-
prob_ax.set_ylim((-0.1, 7))
223-
prob_ax.set_yticks([0, 1])
224-
prob_ax.legend(loc='upper right', fontsize=10)
216+
# if tool == "Dynamont":
217+
# prob_ax = ax[i].twinx()
218+
# borderProbs, acgt_probs = getDynamontProbs()
219+
# # prob_ax.plot(np.arange(len(borderProbs)), borderProbs, c='red', linewidth=0.5, linestyle='-')
220+
# prob_ax.plot(acgt_probs[:, 0], c=basecolors['A'], linewidth=1, linestyle='-', label="A probability")
221+
# prob_ax.plot(acgt_probs[:, 1], c=basecolors['C'], linewidth=1, linestyle='-', label="C probability")
222+
# prob_ax.plot(acgt_probs[:, 2], c=basecolors['G'], linewidth=1, linestyle='-', label="G probability")
223+
# prob_ax.plot(acgt_probs[:, 3], c=basecolors['T'], linewidth=1, linestyle='-', label="T probability")
224+
# prob_ax.set_ylabel("Nucleotide Probabilities", fontsize=10)
225+
# prob_ax.set_ylim((-0.1, 7))
226+
# prob_ax.set_yticks([0, 1])
227+
# prob_ax.legend(loc='upper right', fontsize=10)
225228

226229
plt.tight_layout()
227230
plt.savefig(join(outpath, args.readid + "_tool_segmentation.svg"), dpi=300)
@@ -248,23 +251,25 @@ def main() -> None:
248251
ax[i].vlines([start], ymin=min(signal), ymax=max(signal), colors=basecolors[base], linestyles='--', linewidth=0.7)
249252
ax[i].add_patch(Rectangle((start, min(signal)), end - start, max(signal)-min(signal), alpha=0.4, edgecolor=basecolors[base], facecolor=basecolors[base]))
250253

251-
if tool == "Dynamont":
252-
prob_ax = ax[i].twinx()
253-
# prob_ax.plot(np.arange(len(borderProbs)), borderProbs, c='red', linewidth=0.5, linestyle='-')
254-
prob_ax.plot(acgt_probs[:, 0], c=basecolors['A'], linewidth=1, linestyle='-', label="A probability")
255-
prob_ax.plot(acgt_probs[:, 1], c=basecolors['C'], linewidth=1, linestyle='-', label="C probability")
256-
prob_ax.plot(acgt_probs[:, 2], c=basecolors['G'], linewidth=1, linestyle='-', label="G probability")
257-
prob_ax.plot(acgt_probs[:, 3], c=basecolors['T'], linewidth=1, linestyle='-', label="T probability")
258-
prob_ax.set_ylabel("Nucleotide Probabilities", fontsize=10)
259-
prob_ax.set_ylim((-0.1, 7))
260-
prob_ax.set_yticks([0, 1])
261-
prob_ax.legend(loc='upper right', fontsize=10)
254+
# if tool == "Dynamont":
255+
# prob_ax = ax[i].twinx()
256+
# # prob_ax.plot(np.arange(len(borderProbs)), borderProbs, c='red', linewidth=0.5, linestyle='-')
257+
# prob_ax.plot(acgt_probs[:, 0], c=basecolors['A'], linewidth=1, linestyle='-', label="A probability")
258+
# prob_ax.plot(acgt_probs[:, 1], c=basecolors['C'], linewidth=1, linestyle='-', label="C probability")
259+
# prob_ax.plot(acgt_probs[:, 2], c=basecolors['G'], linewidth=1, linestyle='-', label="G probability")
260+
# prob_ax.plot(acgt_probs[:, 3], c=basecolors['T'], linewidth=1, linestyle='-', label="T probability")
261+
# prob_ax.set_ylabel("Nucleotide Probabilities", fontsize=10)
262+
# prob_ax.set_ylim((-0.1, 7))
263+
# prob_ax.set_yticks([0, 1])
264+
# prob_ax.legend(loc='upper right', fontsize=10)
262265

263266
ax[i].set_xticks(np.arange(0, len(signal), 500))
264267
# ax[i].set_xlim((19000, 22000)) #! for c28
265-
ax[i].set_xlim((7400, 8000))
268+
ax[i].set_xlim((7400, 8000)) #! for paper h_sapiens 131ee77f-...
269+
# ax[i].set_xlim((4500, 6000))
270+
# ax[i].set_ylim((40, 130))
266271

267-
plt.tight_layout()
272+
# plt.tight_layout()
268273
plt.savefig(join(outpath, args.readid + "_tool_segmentation_region.svg"), dpi=300)
269274
plt.savefig(join(outpath, args.readid + "_tool_segmentation_region.pdf"), dpi=300)
270275
plt.savefig(join(outpath, args.readid + "_tool_segmentation_region.png"), dpi=300)

0 commit comments

Comments
 (0)