From dad23302039fda55ca82b708239335ba5a6f28d6 Mon Sep 17 00:00:00 2001 From: nicolesylvester Date: Fri, 12 Jul 2024 16:15:41 -0700 Subject: [PATCH] initial analyze amplicons pipeline --- q2_asap/analyzeAmplicons.py | 17 +++++++++--- q2_asap/analyzeAmplicons_pipeline.py | 30 +++++++++++++++++++++ q2_asap/plugin_setup.py | 38 ++++++++++++++++++++++++--- q2_asap/tests/test_methods.py | 39 +++++++++++----------------- 4 files changed, 92 insertions(+), 32 deletions(-) create mode 100644 q2_asap/analyzeAmplicons_pipeline.py diff --git a/q2_asap/analyzeAmplicons.py b/q2_asap/analyzeAmplicons.py index 65acbc7..df31bf8 100644 --- a/q2_asap/analyzeAmplicons.py +++ b/q2_asap/analyzeAmplicons.py @@ -48,11 +48,20 @@ def analyzeAmplicons(sequences: CasavaOneEightSingleLanePerSampleDirFmt, name: s output_dir_bam = BAMSortedAndIndexedDirFmt() output_dir_bwa = BWAIndexDirFmt() output_dir_xml = ASAPXMLOutputDirFmt() + + + result = subprocess.run(['ls', '-l', f'{sequences}'], stdout=subprocess.PIPE, text=True) + print(result) + print(os.listdir(f'{sequences}')) + subprocess.run(['mkdir', '-p', 'sequences'], check=True) + shutil.copytree(f'{sequences}', 'sequences', dirs_exist_ok=True) # create analyzeAmplicons command - command = f'analyzeAmplicons -r {sequences} -n {name} -o {temp_dir}/asap_output -d {depth} -b {breadth} -j /scratch/nsylvester/SARS2_variant_detection.json \ + command = f'analyzeAmplicons -r sequences -n {name} -o {temp_dir}/asap_output -d {depth} -b {breadth} -j /home/nsylvester/scratch/q2-asap/q2_asap/tests/data/SARS2_variant_detection.json \ --min_base_qual {min_base_qual} --consensus-proportion {consensus_proportion} --fill-gaps {fill_gaps} -a {aligner} --aligner-args {aligner_args}' + print(command) + # combine conda environment and command TODO: fix conda environment shell_script= f""" conda run -p /home/dlemmer/.conda/envs/asap {command} @@ -75,14 +84,14 @@ def analyzeAmplicons(sequences: CasavaOneEightSingleLanePerSampleDirFmt, name: s for file_name in os.listdir(asap_output_dir): file_path = os.path.join(asap_output_dir, file_name) if re.search(r'\.(amb|ann|bwt|pac|sa|fasta)$', file_name): - shutil.move(file_path, Path(output_dir_bwa.path) / file_name) + shutil.copy(file_path, Path(output_dir_bwa.path) / file_name) elif re.search(r'\.xml$', file_name): - shutil.move(file_path, Path(output_dir_xml.path) / file_name) + shutil.copy(file_path, Path(output_dir_xml.path) / file_name) # move bam files for file_name in os.listdir(os.path.join(asap_output_dir, "bwamem")): file_path = os.path.join(asap_output_dir, "bwamem", file_name) - shutil.move(file_path, Path(output_dir_bam.path) / file_name) + shutil.copy(file_path, Path(output_dir_bam.path) / file_name) # remove temp directory # shutil.rmtree(temp_dir) diff --git a/q2_asap/analyzeAmplicons_pipeline.py b/q2_asap/analyzeAmplicons_pipeline.py new file mode 100644 index 0000000..5b8be2b --- /dev/null +++ b/q2_asap/analyzeAmplicons_pipeline.py @@ -0,0 +1,30 @@ +#TODO: config file for bbduk/fastp/trimmer params +#TODO: bamprocessor, outputCombiner +def analyzeAmplicons_pipeline(ctx, sequences, ref_sequence, config_file_path): + + #TODO: long string of params in json file that we will use as trimmer/aligner_params + config = open(config_file_path) + data = json.load(config) + + trimmer = data['ASAP_parameters']['trimmer']['name'] + trimmer_action = ctx.get_action('trimmers', trimmer) + + aligner_index = data['ASAP_parameters']['aligner']['index'] + aligner_index_action = ctx.get_action('aligners', aligner_index) + + aligner = data['ASAP_parameters']['aligner']['name'] + aligner_action = ctx.get_action('aligners', aligner) + + # initialize results + results = [] + + trimmer_results = trimmer_action(sequences) + # results.append(trimmer_results) + + aligner_index_result, = aligner_index_action(ref_sequence) + # results.append(aligner_index_result) + + aligner_result, = aligner_action(sequences = sequences, reference_index = aligner_index_result) + results.append(aligner_result) + + return tuple(results) \ No newline at end of file diff --git a/q2_asap/plugin_setup.py b/q2_asap/plugin_setup.py index d975300..f8135d0 100644 --- a/q2_asap/plugin_setup.py +++ b/q2_asap/plugin_setup.py @@ -18,8 +18,11 @@ from ._formats import ASAPXMLFormat, ASAPHTMLFormat, ASAPXMLOutputDirFmt, ASAPHTMLOutputDirFmt from ._types import ASAPXML, ASAPHTML from q2_nasp2_types.index import BWAIndex -from q2_nasp2_types.alignment import BAM -from q2_types.feature_data import FeatureData +from q2_nasp2_types.alignment import BAMIndexed, SAM +from q2_types.feature_data import FeatureData, Sequence +from q2_asap.analyzeAmplicons_pipeline import analyzeAmplicons_pipeline +from q2_types.bowtie2 import Bowtie2Index + plugin = Plugin( @@ -49,7 +52,7 @@ 'aligner_args': Str }, outputs=[ - ('output_bams', FeatureData[BAM]), + ('output_bams', FeatureData[BAMIndexed]), ('bwa_index', BWAIndex), ('asap_xmls', ASAPXML) ], @@ -59,7 +62,7 @@ 'depth': 'minimum read depth required to consider a position covered. [default: 100]', 'breadth': 'minimum breadth of coverage required to consider an amplicon as present. [default: 0.8]', 'min_base_qual': 'what is the minimum base quality score (BQS) to use a position (Phred scale, i.e. 10=90, 20=99, 30=99.9 accuracy', - 'consensus_proportion': 'minimum proportion required to call at base at that position, else 'N'. [default: 0.8]', + 'consensus_proportion': 'minimum proportion required to call at base at that position, else "N". [default: 0.8]', 'fill_gaps': 'fill no coverage gaps in the consensus sequence [default: False], optional parameter is the character to use for filling [defaut: n]', 'aligner': 'aligner to use for read mapping, supports bowtie2, novoalign, and bwa. [default: bowtie2]', 'aligner_args': "additional arguments to pass to the aligner, enclosed in ''."}, @@ -87,3 +90,30 @@ ) +# register analyzeAmplicon pipeline +plugin.pipelines.register_function( + function=analyzeAmplicons_pipeline, + inputs={ + 'sequences': (SampleData[SequencesWithQuality] | SampleData[PairedEndSequencesWithQuality]), + 'ref_sequence': FeatureData[Sequence], + }, + parameters={'config_file_path': Str}, + outputs=[ + # ('trimmer_results', (SampleData[SequencesWithQuality] | SampleData[PairedEndSequencesWithQuality])), + # ('aligner_index_result', (BWAIndex | Bowtie2Index)), + ('aligner_result', FeatureData[SAM]) + ], + input_descriptions={ + 'sequences': 'the sequences to be analyzed', + 'ref_sequence': 'The reference sequence used to build bwa index.', + }, + parameter_descriptions={'config_file_path': 'path to config file that holds all params'}, + output_descriptions={ + # 'trimmer_results': 'The results after completing trimming', + # 'aligner_index_result': 'The resulting aligner index', + 'aligner_result': 'The result of aligning reads with specified aligner' + }, + name='Analyze Amplicons Pipeline', + description=("A pipeline to run analyze amplicons") +) + diff --git a/q2_asap/tests/test_methods.py b/q2_asap/tests/test_methods.py index b84a72f..7d26bbf 100644 --- a/q2_asap/tests/test_methods.py +++ b/q2_asap/tests/test_methods.py @@ -7,38 +7,29 @@ # ---------------------------------------------------------------------------- import pandas as pd -import pandas.testing as pdt from qiime2.plugin.testing import TestPluginBase from qiime2.plugin.util import transform -from q2_types.feature_table import BIOMV100Format +from qiime2 import Artifact -from q2_asap._methods import duplicate_table +# from q2_asap.analyzeAmplicons_pipeline import analyzeAmplicons_pipeline - -class DuplicateTableTests(TestPluginBase): +class TestAnalyzeAmpliconPipeline(TestPluginBase): package = 'q2_asap.tests' - def test_simple1(self): - in_table = pd.DataFrame( - [[1, 2, 3, 4, 5], [9, 10, 11, 12, 13]], - columns=['abc', 'def', 'jkl', 'mno', 'pqr'], - index=['sample-1', 'sample-2']) - observed = duplicate_table(in_table) - - expected = in_table + def test_analyzeAmplicon_pipeline(self): + # access the pipeline as QIIME 2 sees it, + # for correct assignment of `ctx` variable + analyzeAmplicons_pipeline = self.plugin.pipelines['analyzeAmplicons_pipeline'] - pdt.assert_frame_equal(observed, expected) + #import artifact for reference sequence + ref_sequence_artifact = Artifact.import_data('FeatureData[Sequence]', '/home/nsylvester/scratch/q2-asap/q2_asap/tests/data/wuhan_sequence.fasta') - def test_simple2(self): - # test table duplication with table loaded from file this time - # (for demonstration purposes) - in_table = transform( - self.get_data_path('table-1.biom'), - from_type=BIOMV100Format, - to_type=pd.DataFrame) - observed = duplicate_table(in_table) + # get sequences (paired-end-demux.qza) + sequences_artifact = Artifact.load('q2_asap/tests/data/paired-end-demux.qza') - expected = in_table + # run the pipeline + results = analyzeAmplicons_pipeline(sequences = sequences_artifact, ref_sequence = ref_sequence_artifact) - pdt.assert_frame_equal(observed, expected) + # assert output + self.assertTrue(len(results) > 0)