finalized the script for strand bias

TORCH-Consortium · Jun 18, 2024 · 508d9ff · 508d9ff
1 parent 9271829
commit 508d9ff
Show file tree

Hide file tree

Showing 4 changed files with 20 additions and 9 deletions.
diff --git a/bin/reduce_strand_bias.py b/bin/reduce_strand_bias.py
@@ -1,12 +1,14 @@
 #!/usr/bin/env python3
 
+
 import argparse
 import vcf
-from scipy.stats import binom_test
-
+from scipy.stats import binom_test, binomtest
 
-def filter_vcf_file(vcf_input, vcf_output, pval=0.5):
+import warnings
+warnings.filterwarnings("ignore", category=DeprecationWarning)
 
+def filter_vcf_file(vcf_input, vcf_output, pval):
     # Open the VCF file for reading
 
     vcf_reader = vcf.Reader(open(vcf_input, 'r'))
@@ -33,7 +35,6 @@ def filter_vcf_file(vcf_input, vcf_output, pval=0.5):
             if p_value >= pval:
                 filtered_records.append(record)
 
-
     # Write the filtered records to a new VCF file
     vcf_writer = vcf.Writer(open(vcf_output, 'w'), vcf_reader)
 
@@ -46,11 +47,13 @@ def filter_vcf_file(vcf_input, vcf_output, pval=0.5):
 
 if __name__ == '__main__':
     parser = argparse.ArgumentParser(description='Reduce strand bias in VCF file')
+    parser.add_argument('pval', metavar='pval', type=str, help='The pval to be used for filtering')
     parser.add_argument('input', metavar='input_vcf', type=str, help='The initial VCF file')
     parser.add_argument('output', metavar='output_vcf', type=str, help='The output VCF filed')
     args = vars(parser.parse_args())
 
+    input_pval = float(args['pval'])
     input_vcf_file = args['input']
     output_vcf_file = args['output']
 
-    filter_vcf_file(input_vcf_file, output_vcf_file)
+    filter_vcf_file(input_vcf_file, output_vcf_file, input_pval)
diff --git a/default_params.config b/default_params.config
@@ -39,6 +39,8 @@ ntm_fraction_cutoff = 0.20
 site_representation_cutoff = 0.95
 
 
+strand_bias_cutoff = 0.05
+
 // ##### Partial workflows ##### 
 
 // Set this to true if you'd like to only validate input fastqs and check their FASTQC reports

diff --git a/modules/utils/reformat_lofreq.nf b/modules/utils/reformat_lofreq.nf
@@ -6,20 +6,26 @@ process UTILS_REFORMAT_LOFREQ {
         tuple val(sampleName), path(lofreqVcf)
 
     output:
-        tuple val(sampleName), path("*LoFreq.Reformat.vcf")
+        tuple val(sampleName), path("*lofreq.reformat.corrected.vcf")
 
     script:
 
         """
         reformat_lofreq.py ${lofreqVcf} \\
             ${sampleName} \\
-            ${sampleName}.LoFreq.Reformat.vcf
+            ${sampleName}.lofreq.reformat.vcf
+
+        reduce_strand_bias.py \\
+            ${params.strand_bias_cutoff} \\
+            ${sampleName}.lofreq.reformat.vcf  \\
+            ${sampleName}.lofreq.reformat.corrected.vcf 
         """
 
     stub: 
 
         """
-        touch ${sampleName}.LoFreq.Reformat.vcf
+        touch ${sampleName}.lofreq.reformat.vcf
+        touch ${sampleName}.lofreq.reformat.corrected.vcf 
         """ 
 
 }
diff --git a/workflows/call_wf.nf b/workflows/call_wf.nf
@@ -115,7 +115,7 @@ workflow CALL_WF {
 
         BGZIP__LOFREQ(UTILS_REFORMAT_LOFREQ.out)
 
-        GATK_INDEX_FEATURE_FILE__LOFREQ(BGZIP__LOFREQ.out, 'LoFreq.Reformat')
+        GATK_INDEX_FEATURE_FILE__LOFREQ(BGZIP__LOFREQ.out, 'lofreq.reformat.corrected')
 
         //----------------------------------------------------------------------------------
         // STATS