use het sites

msk-access · May 22, 2024 · 2af3b32 · 2af3b32
1 parent a0750e1
commit 2af3b32
Show file tree

Hide file tree

Showing 3 changed files with 11 additions and 4 deletions.
diff --git a/biometrics/biometrics.py b/biometrics/biometrics.py
@@ -118,7 +118,8 @@ def run_genotyping(args, samples):
         discordance_threshold=args.discordance_threshold,
         threads=args.threads,
         zmin=args.zmin,
-        zmax=args.zmax)
+        zmax=args.zmax,
+        het=args.het)
     cluster_handler = Cluster(args.discordance_threshold)
     comparisons = genotyper.compare_samples(samples)
 

diff --git a/biometrics/cli.py b/biometrics/cli.py
@@ -216,6 +216,9 @@ def get_args():
     parser_genotype.add_argument(
         '--zmax', type=float,
         help='''Maximum z value for the colorscale on the heatmap.''')
+    parser_genotype.add_argument(
+        '--het', type=bool,
+        help='''Use Heterozygous sites too to calculate discordance rate''')
 
     # cluster parser
 

diff --git a/biometrics/genotype.py b/biometrics/genotype.py
@@ -13,14 +13,15 @@
 
 class Genotyper:
 
-    def __init__(self, no_db_compare, discordance_threshold=0.05, threads=1, zmin=None, zmax=None):
+    def __init__(self, no_db_compare, discordance_threshold=0.05, threads=1, zmin=None, zmax=None,het=False):
         self.no_db_compare = no_db_compare
         self.discordance_threshold = discordance_threshold
         self.threads = threads
         self.zmax = zmax
         self.zmin = zmin
         self.sample_type_ratio = 1
         self.comparisons = None
+        self.het=het
 
     def are_samples_same_group(self, sample1, sample2):
 
@@ -253,8 +254,10 @@ def compare_samples(self, samples):
         comparisons = pd.DataFrame(comparisons)
 
         # compute discordance rate
-
-        comparisons['DiscordanceRate'] = comparisons['HomozygousMismatch'] / (comparisons['HomozygousInRef'] + EPSILON)
+        if self.het:
+            comparisons['DiscordanceRate'] = comparisons['HomozygousMismatch'] / (comparisons['HomozygousInRef'] + EPSILON)
+        else:
+            comparisons['DiscordanceRate'] = comparisons['HomozygousMismatch'] + comparisons['HeterozygousMismatch'] / (comparisons['TotalMatch'] + EPSILON)
 
         # data['DiscordanceRate'] = data['DiscordanceRate'].map(lambda x: round(x, 6))
         comparisons.loc[comparisons['HomozygousInRef'] < 10, 'DiscordanceRate'] = np.nan