Skip to content

Commit

Permalink
ycoverage zeroDivision handling
Browse files Browse the repository at this point in the history
  • Loading branch information
mr-eyes committed Oct 28, 2024
1 parent fe02979 commit dce3886
Show file tree
Hide file tree
Showing 2 changed files with 21 additions and 17 deletions.
31 changes: 19 additions & 12 deletions src/snipe/api/multisig_reference_QC.py
Original file line number Diff line number Diff line change
Expand Up @@ -401,7 +401,7 @@ def process_sample(self, sample_sig: SnipeSig, predict_extra_folds: Optional[Lis
})

# ============= GENOME STATS =============

self.logger.debug("Calculating genome statistics.")
# Compute intersection of sample and reference genome
self.logger.debug("Type of sample_sig: %s | Type of reference_sig: %s", sample_sig.sigtype, self.reference_sig.sigtype)
Expand Down Expand Up @@ -480,21 +480,21 @@ def process_sample(self, sample_sig: SnipeSig, predict_extra_folds: Optional[Lis
distance_to_wgs = abs(relative_total_abundance - self.relative_total_abundance_grey_zone[0])
distance_to_wxs = abs(relative_total_abundance - self.relative_total_abundance_grey_zone[1])
predicted_assay_type = "WGS" if distance_to_wgs < distance_to_wxs else "WXS"

self.logger.debug("Predicted assay type: %s", predicted_assay_type)

else:
self.logger.debug("No amplicon signature provided.")

# ============= Contamination/Error STATS =============

self.logger.debug("Calculuating error and contamination indices.")

sample_nonref = sample_sig - self.reference_sig
sample_nonref_singletons = sample_nonref.count_singletons()
sample_nonref_non_singletons = sample_nonref.total_abundance - sample_nonref_singletons
sample_total_abundance = sample_sig.total_abundance

predicted_error_index = (
sample_nonref_singletons / sample_total_abundance
if sample_total_abundance is not None and sample_total_abundance > 0 else 0
Expand Down Expand Up @@ -577,7 +577,7 @@ def process_sample(self, sample_sig: SnipeSig, predict_extra_folds: Optional[Lis
})

advanced_stats.update(amplicon_stats)

# ============= CHR STATS =============

def sort_chromosomes(chrom_dict):
Expand Down Expand Up @@ -681,7 +681,10 @@ def sort_chromosomes(chrom_dict):
self.logger.warning("Autosomal mean abundance is zero. Setting chrX Ploidy score to zero to avoid division by zero.")
xploidy_score = 0.0
else:
xploidy_score = (xchr_mean_abundance / autosomal_mean_abundance) if len(specific_xchr_sig) > 0 else 0.0
xploidy_score = (
(xchr_mean_abundance / autosomal_mean_abundance)
if len(specific_xchr_sig) > 0 and autosomal_mean_abundance > 0 else 0.0
)

self.logger.debug("Calculated chrX Ploidy score: %.4f", xploidy_score)
sex_stats.update({"chrX Ploidy score": xploidy_score})
Expand Down Expand Up @@ -709,9 +712,13 @@ def sort_chromosomes(chrom_dict):
self.logger.warning("Insufficient k-mers for chrY Coverage score calculation. Setting chrY Coverage score to zero.")
ycoverage = 0.0
else:
ycoverage = ((len(ychr_in_sample) / len(ychr_specific_kmers)) / (len(sample_autosomal_sig) / len(autosomals_specific_kmers))
if len(ychr_specific_kmers) > 0 and len(autosomals_specific_kmers) > 0 else 0
try:
ycoverage = (
(len(ychr_in_sample) / len(ychr_specific_kmers)) /
(len(sample_autosomal_sig) / len(autosomals_specific_kmers))
)
except (ZeroDivisionError, TypeError):
ycoverage = 0.0


self.logger.debug("Calculated chrY Coverage score: %.4f", ycoverage)
Expand Down Expand Up @@ -771,7 +778,7 @@ def sort_chromosomes(chrom_dict):
vars_nonref_stats["unexplained variance mean abundance"] = sample_nonref.mean_abundance
vars_nonref_stats["unexplained variance median abundance"] = sample_nonref.median_abundance
vars_nonref_stats["unexplained variance fraction of total abundance"] = (
sample_nonref.total_abundance / sample_nonref_total_abundance
sample_nonref.total_abundance / sample_nonref_total_abundance
if sample_nonref_total_abundance > 0 and sample_nonref.total_abundance is not None else 0.0
)

Expand All @@ -785,7 +792,7 @@ def sort_chromosomes(chrom_dict):

# ============= Coverage Prediction (ROI) =============

if predict_extra_folds and genome_stats["Genome coverage index"]:
if predict_extra_folds and genome_stats["Genome coverage index"] > 0.01:
predicted_fold_coverage = {}
predicted_fold_delta_coverage = {}
nparts = 30
Expand All @@ -798,7 +805,7 @@ def sort_chromosomes(chrom_dict):

# Get sample signature intersected with the reference
_sample_sig_genome = sample_sig & self.reference_sig

hashes = _sample_sig_genome.hashes
abundances = _sample_sig_genome.abundances
N = len(hashes)
Expand Down
7 changes: 2 additions & 5 deletions src/snipe/api/snipe_sig.py
Original file line number Diff line number Diff line change
Expand Up @@ -1449,11 +1449,8 @@ def total_abundance(self) -> int:
int: Total abundance.
"""
self._validate_abundance_operation(None, "calculate total abundance")

total = self.total_abundance
if total == 0:
self.logger.error("Total abundance is zero, which may lead to division by zero.")
raise ZeroDivisionError("Total abundance is zero.")
total = int(np.sum(self._abundances))
self.logger.debug("Total abundance: %d", total)
return total

@property
Expand Down

0 comments on commit dce3886

Please sign in to comment.