Skip to content

Commit

Permalink
Bug fix, scheme was not getting forced
Browse files Browse the repository at this point in the history
  • Loading branch information
martinghunt committed Oct 25, 2024
1 parent 1acbfd7 commit 16026a3
Show file tree
Hide file tree
Showing 2 changed files with 13 additions and 2 deletions.
1 change: 1 addition & 0 deletions viridian/one_sample_pipeline.py
Original file line number Diff line number Diff line change
Expand Up @@ -343,6 +343,7 @@ def detect_amplicon_scheme(self):
min_depth_cutoff=self.coverage_min_x,
min_percent_genome_cutoff=self.coverage_min_pc,
max_primer_dist=constants.SCHEME_ID_PRIMER_WITHIN_END,
force_scheme_choice=self.force_amp_scheme,
)
self.log_dict["reads"] = id_results["read_counts"]
self.log_dict["amplicons"] = id_results["amplicons"]
Expand Down
14 changes: 12 additions & 2 deletions viridian/scheme_id.py
Original file line number Diff line number Diff line change
Expand Up @@ -512,6 +512,7 @@ def analyse_bam(
min_depth_cutoff=20,
min_percent_genome_cutoff=50.0,
min_primer_hits=constants.READ_SAMPLE_MIN_PRIMER_HITS,
force_scheme_choice=None,
):
json_dict = {
"amplicons": None,
Expand Down Expand Up @@ -570,7 +571,11 @@ def analyse_bam(
cumulative_score_plot(schemes, score_plot, title=sample_name)
scores = get_scores_from_schemes(schemes)
json_dict["scheme_choice"] = scores
best_scheme_name = scores["best_scheme"]
if force_scheme_choice is None:
best_scheme_name = scores["best_scheme"]
else:
best_scheme_name = force_scheme_choice

best_scheme = schemes.get(best_scheme_name, None)

per_pos_depth_plot = os.path.join(outdir, "depth_across_genome.pdf")
Expand All @@ -585,7 +590,12 @@ def analyse_bam(

if best_scheme is None:
raise Exception("Error getting best amplicon scheme. Cannot continue")
logging.info(f"{LOG_PREFIX} Best scheme: {best_scheme_name}")
if force_scheme_choice is None:
logging.info(f"{LOG_PREFIX} Best scheme: {best_scheme_name}")
else:
logging.info(
f"{LOG_PREFIX} Using forced scheme: {best_scheme_name} (may not have the best score)"
)
logging.info(f"{LOG_PREFIX} Counting primer matches for scheme {best_scheme_name}")
best_scheme.count_primer_hits(
left_primer_hits, right_primer_hits, max_dist=max_primer_dist
Expand Down

0 comments on commit 16026a3

Please sign in to comment.