Rapid works now

Author: aleksandr-medvedev-msc

find

@@ -0,0 +1,4 @@
+#!/bin/bash
+
+pwd
+docker run --rm -it -v `pwd`/workdir:/media -v /etc/localtime:/etc/localtime:ro genx_relatives:latest launcher.py find "$@"

launcher.py

@@ -201,7 +201,57 @@ def get_parser_args():
                 Smaller values of it tend to give more distant matches than default 500 and more false-positives.
             """
     )
-
+    
+    parser.add_argument(
+        '--rapid-error-rate',
+        default=0.0025,
+        type=float,
+        help="""
+                Minimum number of SNPs in IBD segment.
+                Smaller values of it tend to give more distant matches than default 500 and more false-positives.
+            """
+    )
+        
+    parser.add_argument(
+        '--rapid-min-snp',
+        default=500,
+        type=int,
+        help="""
+                Minimum number of SNPs in IBD segment.
+                Smaller values of it tend to give more distant matches than default 500 and more false-positives.
+            """
+    )
+            
+    parser.add_argument(
+        '--rapid-num-runs',
+        default=10,
+        type=int,
+        help="""
+                Minimum number of SNPs in IBD segment.
+                Smaller values of it tend to give more distant matches than default 500 and more false-positives.
+            """
+    )
+                
+    parser.add_argument(
+        '--rapid-num-success',
+        default=2,
+        type=int,
+        help="""
+                Minimum number of SNPs in IBD segment.
+                Smaller values of it tend to give more distant matches than default 500 and more false-positives.
+            """
+    )
+    
+    parser.add_argument(
+        '--rapid-seg-len',
+        default=5.0,
+        type=float,
+        help="""
+                Minimum number of SNPs in IBD segment.
+                Smaller values of it tend to give more distant matches than default 500 and more false-positives.
+            """
+    )
+    
     parser.add_argument(
         '--stat-file',
         default='stat_file.txt',
@@ -407,8 +457,8 @@ def copy_input(input_dir, working_dir, samples_file):
         config_dict['ref_dir'] = args.ref_directory
     if args.chip:
         config_dict['chip'] = args.chip
-    if args.flow not in ['ibis', 'ibis-king', 'germline-king']:
-        raise ValueError(f'--flow can be one of the ["ibis", "ibis-king", "germline-king"] and not {args.flow}')
+    if args.flow not in ['ibis', 'ibis-king', 'germline-king', 'rapid']:
+        raise ValueError(f'--flow can be one of the ["ibis", "ibis-king", "germline-king", "rapid"] and not {args.flow}')
     config_dict['flow'] = args.flow
     if args.command in ['preprocess', 'simulate', 'reference', 'bundle', 'simbig', 'remove_relatives']:
         config_dict['remove_imputation'] = args.remove_imputation
@@ -421,6 +471,11 @@ def copy_input(input_dir, working_dir, samples_file):
     config_dict['ibis_seg_len'] = args.ibis_seg_len
     config_dict['ibis_min_snp'] = args.ibis_min_snp
 
+    config_dict['rapid_error_rate'] = args.rapid_error_rate
+    config_dict['rapid_min_snp'] = args.rapid_min_snp
+    config_dict['rapid_num_runs'] = args.rapid_num_runs
+    config_dict['rapid_num_success'] = args.rapid_num_success
+    config_dict['rapid_seg_len'] = args.rapid_seg_len
     config_dict['missing_samples'] = args.missing_samples
     config_dict['alt_hom_samples'] = args.alt_hom_samples
     config_dict['het_samples'] = args.het_samples

rules/relatives_rapid.smk

@@ -6,9 +6,9 @@ rule index_and_split:
     input: 
         vcf = "preprocessed/data.vcf.gz"
     output: 
-        "vcf/for_rapid_{chrom}.vcf.gz"
+        vcf = "vcf/for_rapid_{chrom}.vcf.gz"
     conda:
-        "../envs/bcftools.yaml"
+        "bcftools"
     log:
         "logs/vcf/index_and_split-{chrom}.log"
     benchmark:
@@ -25,7 +25,7 @@ rule estimate_rapid_params:
     conda:
         "rapid_params"
     log:
-        "logs/rapid/estimate_rapid_params.log"
+        "logs/rapid/estimate_rapid_params-{chrom}.log"
     output:
         rapid_params = "rapid/params_{chrom}"
     params:
@@ -40,10 +40,10 @@ rule estimate_rapid_params:
 rule interpolate:
     input:
         vcf = rules.index_and_split.output['vcf'],
-        cmmap=cmmap
+        cmmap=CMMAP
     output: "cm/chr{chrom}.cm.g"
     conda:
-        "../envs/interpolation.yaml"
+        "interpolation"
     log:
         "logs/cm/interpolate-{chrom}.log"
     script:
@@ -58,7 +58,7 @@ rule erase_dosages:
     params:
         vcf='vcf/erased_chr{chrom}'
     conda:
-        '../envs/bcftools.yaml'
+        'bcftools'
     shell:
         "bcftools annotate -x 'fmt' {input.vcf} -O z -o {output.vcf}"
 '''
@@ -74,50 +74,32 @@ rule rapid:
         window_size=3,
         output_folder='rapid'
     output:
-        "rapid/results_{chrom}.max.gz"
+        seg = "rapid/results_{chrom}.max.gz"
+    conda:
+        "rapid"
     shell:
         """
-            rapid -i {input.vcf} -g {input.g} -d {params.min_cm_len} -o {params.output_folder} \
-            -w {params.window_size} -r {params.num_runs} -s {params.num_success}"
+            rapidibd -i {input.vcf} -g {input.g} -d {params.min_cm_len} -o {params.output_folder}/chr{wildcards.chrom} \
+            -w {params.window_size} -r {params.num_runs} -s {params.num_success}
         """
 
 
 rule merge_rapid_segments:
     input:
-        seg = expand('rapid/results_{chrom}.max.gz', chrom=CHROMOSOMES)
+        seg = expand('rapid/chr{chrom}/results.max.gz', chrom=CHROMOSOMES)
     output:
         seg = "rapid/results.max"
     shell:
         """
             for file in {input}
             do
-                cat ${{file}} >> {output}
+                zcat ${{file}} >> {output}
             done
         """
 
-checkpoint transform_ibis_segments:
-    input:
-        ibd = rules.merge_rapid_segments.output.seg,
-        fam = "preprocessed/data.fam"
-    output:
-        bucket_dir = directory("ibd")
-    log:
-        "logs/ibis/transform_ibis_segments.log"
-    conda:
-        "evaluation"
-    script:
-        "../scripts/transform_ibis_segments.py"
-
-
-def aggregate_input(wildcards):
-    checkpoints.transform_ibis_segments.get()
-    ids = glob_wildcards(f"ibd/{{id}}.tsv").id
-    return expand(f"ibd/{{id}}.tsv", id=ids)
-
-
 rule ersa:
     input:
-        ibd = aggregate_input
+        ibd = rules.merge_rapid_segments.output['seg']
     output:
         "ersa/relatives.tsv"
     conda:
@@ -155,8 +137,8 @@ rule ersa:
 
 rule postprocess_ersa:
     input:
-        ibd=rules.ibis.output['ibd'],
-        ersa=rules.ersa.output[0]
+        ibd = rules.merge_rapid_segments.output['seg'],
+        ersa = rules.ersa.output[0]
     output: "results/relatives.tsv"
     conda: "evaluation"
     log: "logs/merge/postprocess-ersa.log"

workflows/find/Snakefile

@@ -43,3 +43,5 @@ elif flow == 'ibis-king':
     include: "../../rules/relatives_ibis_king.smk"
 elif flow == 'germline-king':
     include: "../../rules/relatives.smk"
+elif flow == 'rapid':
+    include: "../../rules/relatives_rapid.smk"