diff --git a/environment.yaml b/environment.yaml index 4823d8d..e8cc78a 100644 --- a/environment.yaml +++ b/environment.yaml @@ -8,8 +8,10 @@ dependencies: - bcftools=1.17 - blast=2.14.1 - bwa=0.7.17 + - minimap2=2.28 - samtools=1.17 - spades=3.15.3 + - flye=2.9.4 - clustalw=2.1 - freebayes=1.3.6 - python=3 @@ -18,4 +20,4 @@ dependencies: - pandas=2.0.3 - seaborn=0.12.2 - pip: - - git+https://github.com/BCCDC-PHL/FluViewer.git@v0.1.11-3 + - git+https://github.com/BCCDC-PHL/FluViewer.git diff --git a/fluviewer/analysis.py b/fluviewer/analysis.py index a209177..375a08a 100644 --- a/fluviewer/analysis.py +++ b/fluviewer/analysis.py @@ -97,8 +97,6 @@ def normalize_depth( inputs: dict, outdir: Path, out_name: str, - fwd_reads_raw: Path, - rev_reads_raw: Path, depth: int, max_memory: int, ): @@ -113,10 +111,6 @@ def normalize_depth( :type outdir: str :param out_name: Name of the output directory. :type out_name: str - :param fwd_reads_raw: Path to the raw forward reads. - :type fwd_reads_raw: str - :param rev_reads_raw: Path to the raw reverse reads. - :type rev_reads_raw: str :param depth: Target depth of coverage. :type depth: int :param max_memory: Maximum memory to allocate to BBNorm. @@ -131,6 +125,8 @@ def normalize_depth( input_reads_fwd = inputs.get('input_reads_fwd', None) input_reads_rev = inputs.get('input_reads_rev', None) + input_reads_long = inputs.get('input_reads_long', None) + analysis_summary = { 'timestamp_analysis_start': timestamp_analysis_start, 'inputs': inputs, @@ -138,15 +134,23 @@ def normalize_depth( normalized_reads_fwd = os.path.join(outdir, f'{out_name}-normalized_R1.fastq') normalized_reads_rev = os.path.join(outdir, f'{out_name}-normalized_R2.fastq') - - terminal_command = (f'bbnorm.sh in={input_reads_fwd} in2={input_reads_rev} ' - f'out={normalized_reads_fwd} out2={normalized_reads_rev} target={depth}') + normalized_reads_long = os.path.join(outdir, f'{out_name}-normalized_RL.fastq') + + if input_reads_fwd and input_reads_rev: + terminal_command = (f'bbnorm.sh in={input_reads_fwd} in2={input_reads_rev} ' + f'out={normalized_reads_fwd} out2={normalized_reads_rev} target={depth}') + elif input_reads_long: + terminal_command = (f'bbnorm.sh in={input_reads_long} out={normalized_reads_long} target={depth}') + terminal_command = (terminal_command + f' -Xmx{max_memory}g' if max_memory is not None else terminal_command) # add gzip compression to output files - terminal_command += f'\n\ngzip -f {normalized_reads_fwd}\n' - terminal_command += f'\ngzip -f {normalized_reads_rev}\n' + if input_reads_fwd and input_reads_rev: + terminal_command += f'\n\ngzip -f {normalized_reads_fwd}\n' + terminal_command += f'\ngzip -f {normalized_reads_rev}\n' + elif input_reads_long: + terminal_command += f'\n\ngzip -f {normalized_reads_long}\n' process_name = 'bbnorm' error_code = 2 @@ -179,10 +183,15 @@ def normalize_depth( timestamp_analysis_complete = datetime.datetime.now().isoformat() - outputs = { - 'normalized_reads_fwd': os.path.abspath(normalized_reads_fwd) + '.gz', - 'normalized_reads_rev': os.path.abspath(normalized_reads_rev) + '.gz', - } + if input_reads_fwd and input_reads_rev: + outputs = { + 'normalized_reads_fwd': os.path.abspath(normalized_reads_fwd) + '.gz', + 'normalized_reads_rev': os.path.abspath(normalized_reads_rev) + '.gz', + } + elif input_reads_long: + outputs = { + 'normalized_reads_long': os.path.abspath(normalized_reads_long) + '.gz', + } analysis_summary = { 'process_name': process_name, @@ -204,6 +213,7 @@ def assemble_contigs( inputs: dict, outdir: Path, out_name: str, + threads: int, ): """ Normalized, downsampled reads are assembled de novo into contigs @@ -228,34 +238,67 @@ def assemble_contigs( logs_dir = os.path.join(outdir, 'logs') os.makedirs(logs_dir, exist_ok=True) - - spades_output = os.path.join(outdir, 'spades_output') - fwd_reads = inputs.get('reads_fwd', None) - rev_reads = inputs.get('reads_rev', None) - os.makedirs(spades_output, exist_ok=True) + if 'reads_fwd' in inputs and 'reads_rev' in inputs: + spades_output = os.path.join(outdir, 'spades_output') + fwd_reads = inputs.get('reads_fwd', None) + rev_reads = inputs.get('reads_rev', None) - terminal_command = (f'spades.py --rnaviral --isolate -1 {fwd_reads} ' - f'-2 {rev_reads} -o {spades_output}') + os.makedirs(spades_output, exist_ok=True) - process_name = 'spades' - error_code = 3 + terminal_command = (f'spades.py --threads {threads} --rnaviral -1 {fwd_reads} ' + f'-2 {rev_reads} -o {spades_output}') - script_file = os.path.join(outdir, f'{process_name}_script.sh') - return_code = run(terminal_command, outdir, out_name, process_name, error_code) - analysis_summary['return_code'] = return_code - if not os.path.isfile(os.path.join(spades_output, 'contigs.fasta')): - log.error('No contigs assembled! Aborting analysis.') - error_code = 4 - analysis_summary['return_code'] = error_code - analysis_summary['error_message'] = error_messages_by_code[error_code] - analysis_summary['inputs'] = inputs - return analysis_summary + process_name = 'spades' + error_code = 3 + + script_file = os.path.join(outdir, f'{process_name}_script.sh') + return_code = run(terminal_command, outdir, out_name, process_name, error_code) + analysis_summary['return_code'] = return_code + if not os.path.isfile(os.path.join(spades_output, 'contigs.fasta')): + log.error('No contigs assembled! Aborting analysis.') + error_code = 4 + analysis_summary['return_code'] = error_code + analysis_summary['error_message'] = error_messages_by_code[error_code] + analysis_summary['inputs'] = inputs + return analysis_summary - num_contigs = 0 - src_contigs_path = os.path.join(spades_output, 'contigs.fasta') - dest_contigs_path = os.path.join(outdir, f'{out_name}_contigs.fasta') - shutil.copy(src_contigs_path, dest_contigs_path) + num_contigs = 0 + src_contigs_path = os.path.join(spades_output, 'contigs.fasta') + dest_contigs_path = os.path.join(outdir, f'{out_name}_contigs.fasta') + shutil.copy(src_contigs_path, dest_contigs_path) + + elif 'reads_long' in inputs: + # Note: Trying out flye for long-read assembly. + # May need to adjust approach later, just testing + # this approach for now. + + flye_output = os.path.join(outdir, 'flye_output') + long_reads = inputs.get('reads_long', None) + + os.makedirs(flye_output, exist_ok=True) + + terminal_command = (f'flye --threads {threads} --nano-raw {long_reads} --genome-size 13.5k --out-dir {flye_output}') + + process_name = 'flye' + error_code = 3 + script_file = os.path.join(outdir, f'{process_name}_script.sh') + return_code = run(terminal_command, outdir, out_name, process_name, error_code) + analysis_summary['return_code'] = return_code + + if not os.path.isfile(os.path.join(flye_output, 'assembly.fasta')): + log.error('No contigs assembled! Aborting analysis.') + error_code = 4 + analysis_summary['return_code'] = error_code + analysis_summary['error_message'] = error_messages_by_code[error_code] + analysis_summary['inputs'] = inputs + return analysis_summary + + num_contigs = 0 + src_contigs_path = os.path.join(flye_output, 'assembly.fasta') + dest_contigs_path = os.path.join(outdir, f'{out_name}_contigs.fasta') + shutil.copy(src_contigs_path, dest_contigs_path) + with open(dest_contigs_path, 'r') as f: for line in f: if line.startswith('>'): @@ -1126,7 +1169,7 @@ def make_map_ref(data_frame): return mapping_refs_path -def map_reads(inputs, outdir, out_name, min_qual): +def map_reads(inputs, outdir, out_name, min_qual, threads): """ Reads are mapped to the mapping references (produced by make_mapping_refs func) using BWA mem. The alignment is filtered to retain only paired reads, then sorted and indexed. @@ -1152,64 +1195,103 @@ def map_reads(inputs, outdir, out_name, min_qual): mapping_refs_path = make_mapping_refs(inputs, outdir, out_name) - terminal_command = (f'bwa index {mapping_refs_path}') - process_name = 'bwa_index' - error_code = 14 - return_code = run(terminal_command, outdir, out_name, process_name, error_code) - if return_code != 0: - log.error(f'Error running BWA index (Exit status: {return_code})') - analysis_summary['return_code'] = error_code - analysis_summary['error_message'] = error_messages_by_code[error_code] - return analysis_summary + if 'reads_fwd' in inputs and 'reads_rev' in inputs: + terminal_command = (f'bwa index {mapping_refs_path}') + process_name = 'bwa_index' + error_code = 14 + return_code = run(terminal_command, outdir, out_name, process_name, error_code) + if return_code != 0: + log.error(f'Error running BWA index (Exit status: {return_code})') + analysis_summary['return_code'] = error_code + analysis_summary['error_message'] = error_messages_by_code[error_code] + return analysis_summary - fwd_reads = inputs.get('reads_fwd', None) - rev_reads = inputs.get('reads_rev', None) - alignment_path = os.path.join(outdir, f'{out_name}_alignment.sam') - terminal_command = (f'bwa mem {mapping_refs_path} {fwd_reads} {rev_reads} ' - f'> {alignment_path}') - process_name = 'bwa_mem' - error_code = 15 - return_code = run(terminal_command, outdir, out_name, process_name, error_code) - if return_code != 0: - log.error(f'Error running BWA mem (Exit status: {return_code})') - analysis_summary['return_code'] = error_code - analysis_summary['error_message'] = error_messages_by_code[error_code] - return analysis_summary + fwd_reads = inputs.get('reads_fwd', None) + rev_reads = inputs.get('reads_rev', None) + alignment_path = os.path.join(outdir, f'{out_name}_alignment.sam') + terminal_command = (f'bwa mem -t {threads} {mapping_refs_path} {fwd_reads} {rev_reads} ' + f'> {alignment_path}') + process_name = 'bwa_mem' + error_code = 15 + return_code = run(terminal_command, outdir, out_name, process_name, error_code) + if return_code != 0: + log.error(f'Error running BWA mem (Exit status: {return_code})') + analysis_summary['return_code'] = error_code + analysis_summary['error_message'] = error_messages_by_code[error_code] + return analysis_summary - filtered_alignment_path = os.path.join(outdir, f'{out_name}_alignment.bam') + filtered_alignment_path = os.path.join(outdir, f'{out_name}_alignment.bam') - samtools_require_flags = [ - 'PAIRED', # 0x1 - ] - samtools_require_flags_str = ','.join(samtools_require_flags) - samtools_exclude_flags = [ - 'UNMAP', # 0x4 - 'MUNMAP', # 0x8 - 'SECONDARY', # 0x100 - 'QCFAIL', # 0x200 - 'SUPPLEMENTARY', # 0x800 - ] - samtools_exclude_flags_str = ','.join(samtools_exclude_flags) - - log.info(f'Filtering alignment with sam require flags: {samtools_require_flags_str}.') - log.info(f'Filtering alignment with sam exclude flags: {samtools_exclude_flags_str}.') - log.info(f'See: http://www.htslib.org/doc/samtools-flags.html for info on sam flags.') - log.info('Removing unmapped reads, secondary alignments, and supplementary alignments.') - log.info(f'Applying minimum mapping quality: {min_qual}') - terminal_command = (f'samtools view ' - f'--require-flags {samtools_require_flags_str} ' - f'--exclude-flags {samtools_exclude_flags_str} ' - f'--min-MQ {min_qual} ' - f'--with-header {alignment_path} | samtools sort -o {filtered_alignment_path}') - process_name = 'samtools_view' - error_code = 16 - return_code = run(terminal_command, outdir, out_name, process_name, error_code) - if return_code != 0: - log.error(f'Error running samtools view (Exit status: {return_code})') - analysis_summary['return_code'] = error_code - analysis_summary['error_message'] = error_messages_by_code[error_code] - return analysis_summary + samtools_require_flags = [ + 'PAIRED', # 0x1 + ] + samtools_require_flags_str = ','.join(samtools_require_flags) + samtools_exclude_flags = [ + 'UNMAP', # 0x4 + 'MUNMAP', # 0x8 + 'SECONDARY', # 0x100 + 'QCFAIL', # 0x200 + 'SUPPLEMENTARY', # 0x800 + ] + samtools_exclude_flags_str = ','.join(samtools_exclude_flags) + + log.info(f'Filtering alignment with sam require flags: {samtools_require_flags_str}.') + log.info(f'Filtering alignment with sam exclude flags: {samtools_exclude_flags_str}.') + log.info(f'See: http://www.htslib.org/doc/samtools-flags.html for info on sam flags.') + log.info('Removing unmapped reads, secondary alignments, and supplementary alignments.') + log.info(f'Applying minimum mapping quality: {min_qual}') + terminal_command = (f'samtools view ' + f'--require-flags {samtools_require_flags_str} ' + f'--exclude-flags {samtools_exclude_flags_str} ' + f'--min-MQ {min_qual} ' + f'--with-header {alignment_path} | samtools sort -o {filtered_alignment_path}') + process_name = 'samtools_view' + error_code = 16 + return_code = run(terminal_command, outdir, out_name, process_name, error_code) + if return_code != 0: + log.error(f'Error running samtools view (Exit status: {return_code})') + analysis_summary['return_code'] = error_code + analysis_summary['error_message'] = error_messages_by_code[error_code] + return analysis_summary + elif 'reads_long' in inputs: + reads_long = inputs.get('reads_long', None) + alignment_path = os.path.join(outdir, f'{out_name}_alignment.sam') + terminal_command = (f'minimap2 -t {threads} -ax map-ont {mapping_refs_path} {reads_long} > {alignment_path}') + process_name = 'minimap2' + error_code = 15 + return_code = run(terminal_command, outdir, out_name, process_name, error_code) + if return_code != 0: + log.error(f'Error running minimap2 (Exit status: {return_code})') + analysis_summary['return_code'] = error_code + analysis_summary['error_message'] = error_messages_by_code[error_code] + return analysis_summary + + filtered_alignment_path = os.path.join(outdir, f'{out_name}_alignment.bam') + samtools_exclude_flags = [ + 'UNMAP', # 0x4 + 'SECONDARY', # 0x100 + 'QCFAIL', # 0x200 + 'SUPPLEMENTARY', # 0x800 + ] + samtools_exclude_flags_str = ','.join(samtools_exclude_flags) + log.info(f'Filtering alignment with sam exclude flags: {samtools_exclude_flags_str}.') + log.info(f'See: http://www.htslib.org/doc/samtools-flags.html for info on sam flags.') + log.info('Removing unmapped reads, secondary alignments, and supplementary alignments.') + log.info(f'Applying minimum mapping quality: {min_qual}') + terminal_command = (f'samtools view ' + f'--exclude-flags {samtools_exclude_flags_str} ' + f'--min-MQ {min_qual} ' + f'--with-header {alignment_path} | samtools sort -o {filtered_alignment_path}') + process_name = 'samtools_view' + error_code = 16 + return_code = run(terminal_command, outdir, out_name, process_name, error_code) + if return_code != 0: + log.error(f'Error running samtools view (Exit status: {return_code})') + analysis_summary['return_code'] = error_code + analysis_summary['error_message'] = error_messages_by_code[error_code] + return analysis_summary + log.info(f'Indexing alignment...') terminal_command = (f'samtools index {filtered_alignment_path}') process_name = 'samtools_index' diff --git a/fluviewer/cli_args.py b/fluviewer/cli_args.py index 449c803..28cbdd8 100644 --- a/fluviewer/cli_args.py +++ b/fluviewer/cli_args.py @@ -12,8 +12,9 @@ def parse_args(): :rtype: dict """ parser = argparse.ArgumentParser(description='BCCDC-PHL/FluViewer: Influenza A virus consensus sequence generation and variant calling') - parser.add_argument('-f', '--forward-reads', type=Path, required=True, help='Path to FASTQ file containing forward reads.') - parser.add_argument('-r', '--reverse-reads', type=Path, required=True, help='Path to FASTQ file containing reverse reads.') + parser.add_argument('-f', '--forward-reads', type=Path, help='Path to FASTQ file containing forward reads.') + parser.add_argument('-r', '--reverse-reads', type=Path, help='Path to FASTQ file containing reverse reads.') + parser.add_argument('--long-reads', type=Path, help='Path to FASTQ file containing long reads.') parser.add_argument('-d', '--db', type=Path, required=True, help='Path to FASTA file containing FluViewer database.') parser.add_argument('-o', '--outdir', type=Path, help='Output directory. (default=FluViewer_)') parser.add_argument('-n', '--output-name', type=str, required=True, help='Output name. Includes this name in output files, and in consensus sequence headers.') @@ -47,10 +48,12 @@ def validate_args(args): :type args: argparse.Namespace """ independent_validation_rules = { - 'forward_reads': { 'validation_fn': lambda x: os.path.isfile(x), + 'forward_reads': { 'validation_fn': lambda x: os.path.isfile(x) if x is not None else True, 'error_msg': 'Input file does not exist: {0}' }, - 'reverse_reads': { 'validation_fn': lambda x: os.path.isfile(x), + 'reverse_reads': { 'validation_fn': lambda x: os.path.isfile(x) if x is not None else True, 'error_msg': 'Input file does not exist: {0}' }, + 'long_reads': { 'validation_fn': lambda x: os.path.isfile(x) if x is not None else True, + 'error_msg': 'Input file does not exist: {0}' }, 'db': { 'validation_fn': lambda x: os.path.isfile(x), 'error_msg': 'Input file does not exist: {0}' }, 'outdir': { 'validation_fn': lambda x: True, diff --git a/fluviewer/fluviewer.py b/fluviewer/fluviewer.py index 3be4111..99bf1bc 100644 --- a/fluviewer/fluviewer.py +++ b/fluviewer/fluviewer.py @@ -41,7 +41,10 @@ def main(): exit(1) else: os.makedirs(args.outdir, exist_ok=True) - + + if (args.forward_reads or args.reverse_reads) and args.long_reads: + print('Error: Combining short and long reads is not supported. Please provide either short reads or long reads, but not both.', sys.stderr) + exit(1) log_level = getattr(logging, args.log_level.upper()) print(log_level) @@ -57,8 +60,11 @@ def main(): version_split = version.split('-') log.info(f'Derived from: KevinKuchinski/FluViewer v{version_split[0]}') log.info(f'Inputs:') - log.info(f"Fwd reads: {args.forward_reads}") - log.info(f"Rev reads: {args.reverse_reads}") + if args.forward_reads and args.reverse_reads: + log.info(f"Fwd reads: {args.forward_reads}") + log.info(f"Rev reads: {args.reverse_reads}") + elif args.long_reads: + log.info(f"Long reads: {args.long_reads}") log.info(f"Reference sequences: {args.db}") log.info(f"Outputs:") @@ -75,7 +81,7 @@ def main(): log.info(f"Target depth for pre-normalization of reads: {args.target_depth}") log.info(f"Coverage depth limit for variant calling: {args.coverage_limit}") - + database.check_database( args.db, args.outdir, @@ -110,17 +116,23 @@ def main(): log.info(f'Beginning analysis stage: {current_analysis_stage}') log.info(f'Output directory: {current_analysis_stage_outdir}') - current_analysis_stage_inputs = { - 'input_reads_fwd': os.path.abspath(args.forward_reads), - 'input_reads_rev': os.path.abspath(args.reverse_reads), - } + if args.forward_reads and args.reverse_reads: + current_analysis_stage_inputs = { + 'input_reads_fwd': os.path.abspath(args.forward_reads), + 'input_reads_rev': os.path.abspath(args.reverse_reads), + } + elif args.long_reads: + current_analysis_stage_inputs = { + 'input_reads_long': os.path.abspath(args.long_reads), + } + else: + log.error('No input reads provided.') + exit(1) normalize_depth_analysis_summary = analysis.normalize_depth( current_analysis_stage_inputs, current_analysis_stage_outdir, args.output_name, - os.path.abspath(args.forward_reads), - os.path.abspath(args.reverse_reads), args.target_depth, args.max_memory, ) @@ -158,15 +170,25 @@ def main(): current_analysis_stage_index = analysis_stages.index(current_analysis_stage) if not args.skip_depth_normalization: - current_analysis_stage_inputs = { - 'reads_fwd': normalize_depth_analysis_summary['outputs']['normalized_reads_fwd'], - 'reads_rev': normalize_depth_analysis_summary['outputs']['normalized_reads_rev'], - } + if args.forward_reads and args.reverse_reads: + current_analysis_stage_inputs = { + 'reads_fwd': normalize_depth_analysis_summary['outputs']['normalized_reads_fwd'], + 'reads_rev': normalize_depth_analysis_summary['outputs']['normalized_reads_rev'], + } + elif args.long_reads: + current_analysis_stage_inputs = { + 'reads_long': normalize_depth_analysis_summary['outputs']['normalized_reads_long'], + } else: - current_analysis_stage_inputs = { - 'reads_fwd': os.path.abspath(args.forward_reads), - 'reads_rev': os.path.abspath(args.reverse_reads), - } + if args.forward_reads and args.reverse_reads: + current_analysis_stage_inputs = { + 'reads_fwd': os.path.abspath(args.forward_reads), + 'reads_rev': os.path.abspath(args.reverse_reads), + } + elif args.long_reads: + current_analysis_stage_inputs = { + 'reads_long': os.path.abspath(args.long_reads), + } current_analysis_stage_outdir = os.path.join(args.outdir, 'analysis_by_stage', f'{current_analysis_stage_index:02}_{current_analysis_stage}') current_analysis_stage_outdir = os.path.abspath(current_analysis_stage_outdir) @@ -177,6 +199,7 @@ def main(): current_analysis_stage_inputs, current_analysis_stage_outdir, args.output_name, + args.threads, ) if assemble_contigs_analysis_summary['return_code'] != 0: log.error(f'Error in analysis stage: {current_analysis_stage}') @@ -336,15 +359,25 @@ def main(): } if args.skip_depth_normalization: - current_analysis_stage_inputs.update({ - 'reads_fwd': os.path.abspath(args.forward_reads), - 'reads_rev': os.path.abspath(args.reverse_reads), - }) + if args.forward_reads and args.reverse_reads: + current_analysis_stage_inputs.update({ + 'reads_fwd': os.path.abspath(args.forward_reads), + 'reads_rev': os.path.abspath(args.reverse_reads), + }) + elif args.long_reads: + current_analysis_stage_inputs.update({ + 'reads_long': os.path.abspath(args.long_reads), + }) else: - current_analysis_stage_inputs.update({ - 'reads_fwd': normalize_depth_analysis_summary['outputs']['normalized_reads_fwd'], - 'reads_rev': normalize_depth_analysis_summary['outputs']['normalized_reads_rev'], - }) + if args.forward_reads and args.reverse_reads: + current_analysis_stage_inputs.update({ + 'reads_fwd': normalize_depth_analysis_summary['outputs']['normalized_reads_fwd'], + 'reads_rev': normalize_depth_analysis_summary['outputs']['normalized_reads_rev'], + }) + elif args.long_reads: + current_analysis_stage_inputs.update({ + 'reads_long': normalize_depth_analysis_summary['outputs']['normalized_reads_long'], + }) log.info(f'Beginning analysis stage: {current_analysis_stage}') @@ -353,6 +386,7 @@ def main(): current_analysis_stage_outdir, args.output_name, args.min_mapping_quality, + args.threads, ) if map_reads_analysis_summary['return_code'] != 0: log.error(f'Error in analysis stage: {current_analysis_stage}')