Workflow improvement: support unaligned FASTQ inputs

workflow/leviosam2.py

@@ -0,0 +1,255 @@
+'''
+Nae-Chyun Chen
+Johns Hopkins University
+2021-2022
+'''
+
+import argparse
+import pathlib
+import subprocess
+import typing
+
+
+def parse_args() -> argparse.Namespace:
+    parser = argparse.ArgumentParser()
+    parser.add_argument('--leviosam2_binary',
+                        default='leviosam2',
+                        type=str,
+                        help='Path to the leviosam2 executable')
+    parser.add_argument('--samtools_binary',
+                        default='samtools',
+                        type=str,
+                        help='Path to the samtools executable')
+    parser.add_argument('--gnu_time_binary',
+                        default='gtime',
+                        type=str,
+                        help=('Path to the GNU Time executable '
+                              '(see https://www.gnu.org/software/time/)'))
+    parser.add_argument('--measure_time',
+                        action='store_true',
+                        help=('Activate to measure time with '
+                              '`--gnu_time_binary`'))
+    parser.add_argument('--keep_tmp_files',
+                        action='store_true',
+                        help=('Activate to keep temp files '
+                              'generated by the workflow'))
+    parser.add_argument('-t',
+                        '--num_threads',
+                        type=int,
+                        default=8,
+                        help='Number of threads to use')
+    parser.add_argument('--sequence_type',
+                        type=str,
+                        required=True,
+                        choices=['ilmn_pe', 'ilmn_se', 'pb_hifi', 'ont'],
+                        help='Type of sequence data')
+    parser.add_argument(
+        '--aligner',
+        type=str,
+        required=True,
+        choices=['bowtie2', 'bwamem', 'bwamem2', 'minimap2', 'winnowmap2'])
+    parser.add_argument('--source_label',
+                        type=str,
+                        default='source',
+                        help='Label of the source reference')
+    parser.add_argument('--target_label',
+                        type=str,
+                        default='target',
+                        help='Label of the target reference')
+    parser.add_argument('--read_group', type=str, help='Read group string')
+    parser.add_argument('-g',
+                        '--lift_max_gap',
+                        type=int,
+                        help='[lift] Max chain gap size allowed')
+    parser.add_argument('--lift_commit_min_mapq',
+                        type=int,
+                        help='[lift] Min MAPQ to commit')
+    parser.add_argument('--lift_commit_min_score',
+                        type=int,
+                        help='[lift] Min alignment score (AS:i tag) to commit')
+    parser.add_argument('--lift_commit_max_frac_clipped',
+                        type=float,
+                        help='[lift] Min fraction of clipped bases to commit')
+    parser.add_argument('--lift_commit_max_isize',
+                        type=int,
+                        help='[lift] Max template length (isize) to commit')
+    parser.add_argument('--lift_commit_max_hdist',
+                        type=int,
+                        help='[lift] Max edit distance (NM:i tag) to commit')
+    parser.add_argument('--lift_bed_commit_source',
+                        type=str,
+                        help=('[lift] Path to a BED (source coordinates) '
+                              'where reads in the regions are always '
+                              'committed (often for suppress annotations)'))
+    parser.add_argument('--lift_bed_defer_target',
+                        type=str,
+                        help=('[lift] Path to a BED (target cooridnates'
+                              'where reads in the regions are always '
+                              'deferred'))
+    parser.add_argument('--lift_realign_config',
+                        type=str,
+                        help=('[lift] Path to the config file for '
+                              'realignment'))
+    parser.add_argument('-i',
+                        '--input_alignment',
+                        type=str,
+                        required=True,
+                        help='Path to the input SAM/BAM/CRAM file')
+    parser.add_argument('-o',
+                        '--out_prefix',
+                        type=str,
+                        required=True,
+                        help='Output prefix')
+    parser.add_argument('-C',
+                        '--leviosam2_index',
+                        type=str,
+                        required=True,
+                        help='Path to the leviosam2 index')
+    parser.add_argument('-f',
+                        '--target_fasta',
+                        type=str,
+                        required=True,
+                        help='Path to the target reference (FASTA file)')
+    parser.add_argument('--dryrun',
+                        action='store_true',
+                        help='Activate the dryrun mode')
+    parser.add_argument('--forcerun',
+                        action='store_true',
+                        help='Activate the forcerun mode. Rerun everything')
+    # parser.add_argument()
+    # parser.add_argument()
+    # parser.add_argument()
+    # parser.add_argument()
+    # parser.add_argument()
+
+    args = parser.parse_args()
+    return args
+
+
+def run_leviosam2(time_cmd: str, leviosam2: str, clft: str, fn_input: str,
+                  out_prefix: str, num_threads: int,
+                  lift_commit_min_mapq: typing.Union[int, None],
+                  lift_commit_min_score: typing.Union[int, None],
+                  lift_commit_max_frac_clipped: typing.Union[float, None],
+                  lift_commit_max_isize: typing.Union[int, None],
+                  lift_commit_max_hdist: typing.Union[int, None],
+                  lift_max_gap: typing.Union[int, None],
+                  lift_bed_commit_source: str, lift_bed_defer_target: str,
+                  lift_realign_config: str, target_fasta: str, dryrun: bool):
+    lift_commit_min_mapq_arg = f'-S mapq:{lift_commit_min_mapq} ' \
+        if lift_commit_min_mapq else ''
+    lift_commit_min_score_arg = f'-S aln_score:{lift_commit_min_score} ' \
+        if lift_commit_min_score else ''
+    lift_commit_max_frac_clipped_arg = \
+        f'-S clipped_frac:{lift_commit_max_frac_clipped} ' \
+            if lift_commit_max_frac_clipped else ''
+    lift_commit_max_isize_arg = f'-S isize:{lift_commit_max_isize} ' \
+        if lift_commit_max_isize else ''
+    lift_commit_max_hdist_arg = f'-S hdist:{lift_commit_max_hdist} ' \
+        if lift_commit_max_hdist else ''
+    lift_max_gap_arg = f'-G {lift_max_gap} ' if lift_max_gap else ''
+    lift_bed_commit_source_arg = f'-r {lift_bed_commit_source} ' \
+        if lift_bed_commit_source else ''
+    lift_bed_defer_target_arg = f'-D {lift_bed_defer_target} ' \
+        if lift_bed_defer_target else ''
+    lift_realign_config_arg = f'-x {lift_realign_config} ' \
+        if lift_realign_config else ''
+
+    cmd = (f'{time_cmd} {leviosam2} lift -C {clft} '
+           f'-a {fn_input} -p {out_prefix}'
+           f'-t {num_threads} -m -f {target_fasta} '
+           f'{lift_commit_min_mapq_arg}'
+           f'{lift_commit_min_score_arg}'
+           f'{lift_commit_max_frac_clipped_arg}'
+           f'{lift_commit_max_hdist_arg}'
+           f'{lift_commit_max_isize_arg}'
+           f'{lift_max_gap_arg}'
+           f'{lift_bed_commit_source_arg}'
+           f'{lift_bed_defer_target_arg}'
+           f'{lift_realign_config_arg}')
+    if dryrun:
+        print(cmd)
+    else:
+        subprocess.run([cmd], shell=True)
+
+
+def run_sort_committed(time_cmd: str, samtools: str, num_threads: int,
+                       out_prefix: str, dryrun: bool, forcerun: bool):
+    fn_committed = pathlib.Path(f'{out_prefix}-committed.bam')
+    fn_committed_sorted = pathlib.Path(f'{out_prefix}-committed-sorted.bam')
+
+    cmd = (f'{time_cmd} {samtools} sort -@ {num_threads} '
+           f'-o {fn_committed_sorted} {fn_committed}')
+    if dryrun:
+        print(cmd)
+    else:
+        if not fn_committed.is_file():
+            raise FileNotFoundError(f'{fn_committed} is not a file')
+        if (not forcerun) and fn_committed_sorted.is_file():
+            print(f'[Info] Skip sort_committed -- '
+                  f'{fn_committed_sorted} exists')
+            return
+        subprocess.run([cmd], shell=True)
+
+
+def run_collate():
+    '''
+    if [ ! -s ${PFX}-paired-deferred-R1.fq.gz ]; then
+        ${MT} ${LEVIOSAM} collate \
+        -a ${PFX}-committed-sorted.bam -b ${PFX}-deferred.bam -p ${PFX}-paired
+    fi
+    '''
+    pass
+
+
+def run_workflow(args: argparse.Namespace):
+    time_cmd = ''
+    if args.measure_time:
+        time_cmd = f'{args.gnu_time_binary} -v -ao {args.out_prefix}.time_log'
+    run_leviosam2(
+        time_cmd=time_cmd,
+        leviosam2=args.leviosam2_binary,
+        clft=args.leviosam2_index,
+        fn_input=args.input_alignment,
+        out_prefix=args.out_prefix,
+        num_threads=args.num_threads,
+        lift_commit_min_mapq=args.lift_commit_min_mapq,
+        lift_commit_min_score=args.lift_commit_min_score,
+        lift_commit_max_frac_clipped=args.lift_commit_max_frac_clipped,
+        lift_commit_max_isize=args.lift_commit_max_isize,
+        lift_commit_max_hdist=args.lift_commit_max_hdist,
+        lift_max_gap=args.lift_max_gap,
+        lift_bed_commit_source=args.lift_bed_commit_source,
+        lift_bed_defer_target=args.lift_bed_defer_target,
+        lift_realign_config=args.lift_realign_config,
+        target_fasta=args.target_fasta,
+        dryrun=args.dryrun)
+
+    run_sort_committed(time_cmd=time_cmd,
+                       samtools=args.samtools_binary,
+                       num_threads=args.num_threads,
+                       out_prefix=args.out_prefix,
+                       dryrun=args.dryrun,
+                       forcerun=args.forcerun)
+
+    if args.sequence_type in ['ilmn_pe']:
+        run_collate_pe()
+        run_realign_deferred_pe()
+        run_refflow_merge_pe()
+        run_merge_pe()
+    else:
+        run_bam_to_fastq_se()
+        run_realign_deferred_se()
+        run_merge_se()
+
+    run_index()
+    run_clean()
+
+
+def validate_binary():
+    args.gnu_time_binary
+
+
+if __name__ == '__main__':
+    args = parse_args()
+    run_workflow(args)