Merge pull request #196 from dvitale199/pipeline

adding option to maf prune during lambda calculation

Author: dvitale199

genotools/gwas.py

@@ -27,7 +27,7 @@
 plink2_exec = check_plink2()
 
 class Assoc:
-    def __init__(self, geno_path=None, out_path=None, pca=10, build='hg38', gwas=True, pheno_name='PHENO1', covar_path=None, covar_names=None):
+    def __init__(self, geno_path=None, out_path=None, pca=10, build='hg38', gwas=True, pheno_name='PHENO1', covar_path=None, covar_names=None, maf_lambdas=False):
         self.geno_path = geno_path
         self.out_path = out_path
         self.pca = pca
@@ -36,6 +36,7 @@ def __init__(self, geno_path=None, out_path=None, pca=10, build='hg38', gwas=Tru
         self.pheno_name = pheno_name
         self.covar_path = covar_path
         self.covar_names = covar_names
+        self.maf_lambdas = maf_lambdas
 
 
     def write_exclusion_file(self):
@@ -276,16 +277,34 @@ def run_gwas(self, covars=True):
             # add pruning step here (pre lambdas)
             gwas_df_add = gwas_df.loc[gwas_df.TEST=='ADD']
 
+            if self.maf_lambdas:
+                gwas_df_maf = gwas_df_add[gwas_df_add['A1_FREQ']>=0.01] 
+                gwas_df_maf = gwas_df_maf[gwas_df_maf['A1_FREQ']<=0.99]
+                gwas_df_maf = gwas_df_maf[~gwas_df_maf['P'].isna()]
+
+                lambda_maf_dict = Assoc.calculate_inflation(gwas_df_maf.P, normalize=False)
+                lambda1000_maf_dict = Assoc.calculate_inflation(gwas_df_maf.P, normalize=True, ncases=ncases, ncontrols=ncontrols)
+
             # calculate inflation
             lambda_dict = Assoc.calculate_inflation(gwas_df_add.P, normalize=False)
             lambda1000_dict = Assoc.calculate_inflation(gwas_df_add.P, normalize=True, ncases=ncases, ncontrols=ncontrols)
 
-            metrics_dict = {
-                'lambda': lambda_dict['metrics']['inflation'],
-                'lambda1000': lambda1000_dict['metrics']['inflation'],
-                'cases': ncases,
-                'controls': ncontrols
+            if self.maf_lambdas:
+                metrics_dict = {
+                    'lambda': lambda_dict['metrics']['inflation'],
+                    'lambda1000': lambda1000_dict['metrics']['inflation'],
+                    'lambda_maf': lambda_maf_dict['metrics']['inflation'],
+                    'lambda1000_maf': lambda1000_maf_dict['metrics']['inflation'],
+                    'cases': ncases,
+                    'controls': ncontrols
                 }
+            else:
+                metrics_dict = {
+                    'lambda': lambda_dict['metrics']['inflation'],
+                    'lambda1000': lambda1000_dict['metrics']['inflation'],
+                    'cases': ncases,
+                    'controls': ncontrols
+                    }
 
             process_complete = True
 

genotools/pipeline.py

@@ -76,6 +76,7 @@ def gt_argparse():
     parser.add_argument('--gwas', type=str, nargs='?', default='False', const='True', help='Run GWAS')
     parser.add_argument('--covars', type=str, nargs='?', default=None, const=None, help='Path to external covars')
     parser.add_argument('--covar_names', type=str, nargs='?', default=None, const=None, help='Covar names to use from external file')
+    parser.add_argument('--maf_lambdas', type=str, nargs='?', default='False', const='True', help='MAF prune before lambda calculations')
 
 
     # parse args and turn into dict
@@ -219,6 +220,7 @@ def execute_pipeline(steps, steps_dict, geno_path, out_path, samp_qc, var_qc, as
                 assoc.gwas = args['gwas']
                 assoc.covar_path = args['covars']
                 assoc.covar_names = args['covar_names']
+                assoc.maf_lambdas = args['maf_lambdas']
                 out_dict[step] = steps_dict[step]()
 
             pass_fail[step] = {'status':out_dict[step]['pass'], 'input':step_input, 'output':step_output}