diff --git a/.gitignore b/.gitignore index 9285a09..354c86a 100644 --- a/.gitignore +++ b/.gitignore @@ -185,4 +185,5 @@ oryza-terrier.final.TEs.csv drosophila-terrier.final.TEs.csv lightning_logs/ comparison-test-data/get-rules.ipynb -compare.sh \ No newline at end of file +compare.sh +drosophila.final.TEs.fa.gz \ No newline at end of file diff --git a/README.rst b/README.rst index 95456ce..5168b90 100644 --- a/README.rst +++ b/README.rst @@ -78,7 +78,7 @@ To run inference on a FASTA file, run this command: .. code-block:: bash - terrier --file INPUT.fa --output-fasta OUTPUT.fa + terrier --input INPUT.fa --output-fasta OUTPUT.fa That will add the classification to after the sequence ID in the `OUTPUT.fa` FASTA file. @@ -86,15 +86,21 @@ If you want to save the probabilities for all classes run this: .. code-block:: bash - terrier --file INPUT.fa --output-csv OUTPUT.csv + terrier --input INPUT.fa --output-csv OUTPUT.csv The columns will be the probability of each classification and the rows correspond to each sequence in ``INPUT.fa``. +You can also use a URL as the input: + +.. code-block:: bash + + terrier --input https://example.com/INPUT.fasta.gz --output-fasta OUTPUT.fa + If you want to output a visualization of the prediction probabilities: .. code-block:: bash - terrier --file INPUT.fa --image-dir OUTPUT-IMAGES/ + terrier --input INPUT.fa --image-dir OUTPUT-IMAGES/ The outputs for the above can be combined together. For more options run diff --git a/docs/preprocessing.rst b/docs/preprocessing.rst index 44870d5..8f71575 100644 --- a/docs/preprocessing.rst +++ b/docs/preprocessing.rst @@ -28,7 +28,7 @@ These two files can be generated from the Repbase database using the ``terrier-t .. code-block:: bash - terrier-tools preprocess --repbase $REPBASE_DIR --seqbank $REPBASE_DIR/Repbase-seqbank.sb --seqtree $REPBASE_DIR/Repbase-seqtree.st + terrier-tools preprocess --input $REPBASE_DIR --seqbank $REPBASE_DIR/Repbase-seqbank.sb --seqtree $REPBASE_DIR/Repbase-seqtree.st This will create a SeqBank file called ``Repbase-seqbank.sb`` and a SeqTree files called ``Repbase-seqtree.st`` and place them the ``$REPBASE_DIR``. @@ -46,13 +46,13 @@ Now you are ready to train Terrier using the SeqBank and SeqTree files you have Optional: Display the SeqTree ------------------------------ -You can list the number of accessions for each node in the SeqTree file with this command: +You can list the number of sequences for each node in the SeqTree file with this command: .. code-block:: bash seqtree render $REPBASE_DIR/Repbase-seqtree.st --print --count -That will output a tree with the number of accessions like this: +That will output a tree with the number of sequences like this: .. code-block:: text @@ -124,4 +124,33 @@ This will create an HTML file with the Sunburst chart of the SeqTree like this: You can open the HTML file in a browser to view the chart. -You can also output the SeqTree with a .png, .svg, or .pdf extension by changing the extension of the output file. \ No newline at end of file +You can also output the SeqTree with a .png, .svg, or .pdf extension by changing the extension of the output file. + +Custom Datasets +---------------- + +You can create a custom repeat library in FASTA format, with the classification of each sequence like this: + +.. code-block:: text + + >SeqID#DNA/Academ + ACTGACTGACTG... + +Or with the classification separated with a tab character like this: + +.. code-block:: text + + >SeqID LTR/Caulimovirus + ACTGACTGACTG... + +Then preprocess like this: + +.. code-block:: bash + + terrier-tools preprocess --input custom.fasta --seqbank custom-seqbank.sb --seqtree custom-seqtree.st + +You can include Repbase with your custom dataset like this: + +.. code-block:: bash + + terrier-tools preprocess --input $REPBASE_DIR --input custom.fasta --seqbank combined-seqbank.sb --seqtree combined-seqtree.st \ No newline at end of file diff --git a/docs/reproduction.rst b/docs/reproduction.rst index 2c03567..ab1669c 100644 --- a/docs/reproduction.rst +++ b/docs/reproduction.rst @@ -15,15 +15,12 @@ Fruit Fly Genome Bickmann et al. (2023) provide Transposable Elements (TE) models of a fruit-fly genome. Download it with the following command: -.. code-block:: bash - - wget https://raw.githubusercontent.com/IOB-Muenster/TEclass2/refs/heads/main/tests/Drosophila_melanogaster.fasta Run inference using Terrier like this: .. code-block:: bash - terrier --file Drosophila_melanogaster.fasta \ + terrier --input https://raw.githubusercontent.com/IOB-Muenster/TEclass2/refs/heads/main/tests/drosophila.final.TEs.fa \ --output-csv drosophila-terrier.final.TEs.csv \ --min-length 0 \ --threshold 0 @@ -108,15 +105,12 @@ Rice Genome Bickmann et al. (2023) also provide Transposable Elements (TE) models of a rice genome. Download it with the following command: -.. code-block:: bash - - wget https://raw.githubusercontent.com/IOB-Muenster/TEclass2/refs/heads/main/tests/Oryza_sativa.fasta - Run inference using Terrier like this: .. code-block:: bash - terrier --file Oryza_sativa.fasta --output-csv oryza-terrier.final.TEs.csv --threshold 0 + terrier --intput https://raw.githubusercontent.com/IOB-Muenster/TEclass2/refs/heads/main/tests/oryza.final.TEs.fa \ + --output-csv oryza-terrier.final.TEs.csv --threshold 0 .. note:: @@ -203,7 +197,7 @@ Run inference using Terrier like this: .. code-block:: bash - terrier --file Homo_sapiens.fasta --output-csv Terrier-human.csv --threshold 0 + terrier --input Homo_sapiens.fasta --output-csv Terrier-human.csv --threshold 0 Now evaluate the results with the following command: @@ -243,7 +237,7 @@ Run inference using Terrier like this: .. code-block:: bash - terrier --file Mus_musculus.fasta --output-csv Terrier-mouse.csv --threshold 0 + terrier --input Mus_musculus.fasta --output-csv Terrier-mouse.csv --threshold 0 Now evaluate the results with the following command: diff --git a/docs/training.rst b/docs/training.rst index 3a7bc82..759d2cc 100644 --- a/docs/training.rst +++ b/docs/training.rst @@ -7,7 +7,7 @@ After performing the instructions on the :ref:`preprocessing:Preprocessing` page To train Terrier, you will need to use the `terrier-tools` CLI utility. -To use the same hyperparameters as in the main release of Terrier, you can run the following command: +To train with the default settings of Terrier, you can run the following command: .. code-block:: bash @@ -17,6 +17,19 @@ To use the same hyperparameters as in the main release of Terrier, you can run t --seqtree $SEQTREE \ --seqbank $SEQBANK +If you want to train using the pretrained Terrier model weights as a starting point, you can add the ``--pretrained`` flag: + +.. code-block:: bash + + SEQBANK=$REPBASE_DIR/Repbase-seqbank.sb + SEQTREE=$REPBASE_DIR/Repbase-seqtree.st + terrier-tools train \ + --seqtree $SEQTREE \ + --seqbank $SEQBANK \ + --pretrained default + +You can replace the word ``default`` with a path to a checkpoint file if you have one or to a URL to a checkpoint file. + You can see other command-line options by running: .. code-block:: bash diff --git a/poetry.lock b/poetry.lock index 591925a..14e7905 100644 --- a/poetry.lock +++ b/poetry.lock @@ -375,14 +375,14 @@ chardet = ">=3.0.2" [[package]] name = "bio-corgi" -version = "0.5.0a2" +version = "0.5.1" description = "Classifier for ORganelle Genomes Inter alia" optional = false python-versions = "<3.13,>=3.10" groups = ["main"] files = [ - {file = "bio_corgi-0.5.0a2-py3-none-any.whl", hash = "sha256:ead9af52004ba7f581080819456c6aa157b4068683fb7c8551fcad6d81a6e1f1"}, - {file = "bio_corgi-0.5.0a2.tar.gz", hash = "sha256:4c1a6839fef5efff05f475e36b44c34b0939a4c504e5a96f03e1ca2ebdbd6e4a"}, + {file = "bio_corgi-0.5.1-py3-none-any.whl", hash = "sha256:2f3c513e995ffada505c6a63b1db81bad718dc1e4ff39fcb2b9c3a3c95b5edb2"}, + {file = "bio_corgi-0.5.1.tar.gz", hash = "sha256:2d78ab8b6e807c47eb0e63bb8f2a4e065ae9d9a2b3a6a514169391b946e0723b"}, ] [package.dependencies] @@ -5804,4 +5804,4 @@ propcache = ">=0.2.1" [metadata] lock-version = "2.1" python-versions = ">=3.10,<3.13" -content-hash = "30d6d7342864455c1ffe49ced6bbc7b9ce977ff624ad83910d6e771084e45b16" +content-hash = "833ed9506890909ce15aea6a429061fff472694cf28d0fd4f4e0d68755a8ae2c" diff --git a/pyproject.toml b/pyproject.toml index bfe9f74..b533720 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [tool.poetry] name = "bio-terrier" -version = "0.3.3" +version = "0.3.4" description = "Transposable Element Repeat Result classifIER" authors = ["Robert Turnbull "] license = "Apache-2.0" @@ -26,7 +26,7 @@ python = ">=3.10,<3.13" numpy = "<2.0.0" pandas = "<=2.2.0" h5py = ">=3.8.0" -bio-corgi = ">=0.5.0a2" +bio-corgi = ">=0.5.1" pyfastx = ">=1.1.0" toml = ">=0.10.2" diff --git a/terrier/apps.py b/terrier/apps.py index 55c605d..cb25c5e 100644 --- a/terrier/apps.py +++ b/terrier/apps.py @@ -219,23 +219,12 @@ def get_prediction_probability(row): def checkpoint(self, checkpoint:Path=None) -> str: return checkpoint or "https://github.com/rbturnbull/terrier/releases/download/v0.2.0/terrier-0.2.0.ckpt" - @ta.tool - def create_repeatmasker_seqtree(self, output:Path, repbase:Path, label_smoothing:float=0.0, gamma:float=0.0, partitions:int=5): - from .repeatmasker import create_repeatmasker_seqtree - return create_repeatmasker_seqtree( - output=output, - repbase=repbase, - label_smoothing=label_smoothing, - gamma=gamma, - partitions=partitions, - ) - @ta.tool def preprocess( self, - repbase:Path=ta.Param(..., help="The path to the RepBase fasta directory."), - seqbank:Path=ta.Param(..., help="The path to save the new SeqBank file."), - seqtree:Path=ta.Param(..., help="The path to save the new SeqTree file."), + input:list[str]=ta.Param(..., help="The path to a FASTA file, URL to a FASTA file, multiple FASTA files or a directory of FASTA files (e.g. the RepBase FASTA directory)"), + seqbank:Path=ta.Param(None, help="The path to save the new SeqBank file."), + seqtree:Path=ta.Param(None, help="The path to save the new SeqTree file."), label_smoothing:float=0.0, gamma:float=0.0, partitions:int=5, @@ -243,23 +232,26 @@ def preprocess( from seqbank import SeqBank from .repeatmasker import create_repeatmasker_seqtree - seqbank = SeqBank(path=seqbank, write=True) - assert repbase is not None - repbase = Path(repbase) - assert repbase.exists() + assert seqbank or seqtree, "You must provide either a --seqbank or --seqtree output path (usually both)." - # Create the seqbank from the FASTA files with .ref extension - files = list(repbase.glob('*.ref')) - seqbank.add_files(files, format="fasta") + fasta_paths = self.find_fasta_paths(input) + + # Create the seqbank from the FASTA files + if seqbank: + seqbank = SeqBank(path=seqbank, write=True) + seqbank.add_files(fasta_paths, format="fasta") # Create the seqtree - return create_repeatmasker_seqtree( - output=seqtree, - repbase=repbase, - label_smoothing=label_smoothing, - gamma=gamma, - partitions=partitions, - ) + if seqtree: + seqtree_path = Path(seqtree) + seqtree = create_repeatmasker_seqtree( + fasta_paths=fasta_paths, + label_smoothing=label_smoothing, + gamma=gamma, + partitions=partitions, + ) + seqtree.save(seqtree_path) + seqtree.render(print=1, count=True) @ta.tool def evaluate( diff --git a/terrier/repeatmasker.py b/terrier/repeatmasker.py index c7e6108..5363b20 100644 --- a/terrier/repeatmasker.py +++ b/terrier/repeatmasker.py @@ -4,9 +4,42 @@ import toml from corgi.seqtree import SeqTree from collections import Counter +import gzip -def create_repeatmasker_seqtree(output:Path, repbase:Path, label_smoothing:float=0.0, gamma:float=0.0, partitions:int=5): +def open_maybe_gz(file:Path): + if file.name.endswith('.gz'): + return gzip.open(file, "rt") + else: + return open(file, "r") + + +def get_verbatim_classification(path:Path, record) -> str: + accession = record.id + + description = record.description + + if "#" in description: + description = description[ description.find("#")+1 : ] + + components = description.split("\t") + if len(components) >= 2: + return components[1] + + if path.name == "simple.ref": + return "Simple Repeat" + elif accession.startswith("SINE_"): + return "SINE" + else: + return description + + +def create_repeatmasker_seqtree( + fasta_paths:list[Path], + label_smoothing:float=0.0, + gamma:float=0.0, + partitions:int=5, +) -> SeqTree: with open(Path(__file__).parent/"data/repbase-to-repeatmasker.toml", "r") as f: mapping = toml.load(f) @@ -21,45 +54,49 @@ def create_repeatmasker_seqtree(output:Path, repbase:Path, label_smoothing:float # Read files count = 0 - for file in repbase.glob('*.ref'): - with open(file) as f: + for file in fasta_paths: + with open_maybe_gz(file) as f: for record in SeqIO.parse(f, "fasta"): partition = count % partitions accession = record.id - components = record.description.split("\t") - if len(components) != 3: - if file.name == "simple.ref": - classification = "Simple Repeat" - elif accession.startswith("SINE_"): - classification = "SINE" - else: - continue - else: - classification = components[1] - - if classification not in mapping: + classification = get_verbatim_classification(file, record) + + if classification in mapping: + mapped_counter.update([classification]) + classification = mapping[classification] + + if classification not in mapping.values(): not_mapped_counter.update([classification]) continue - mapped_counter.update([classification]) - - repeat_name = mapping[classification] - if repeat_name == "Unknown": + if classification == "Unknown": continue - if repeat_name not in classification_nodes: - components = repeat_name.split("/") + if classification not in classification_nodes: + components = classification.split("/") repeat_type = components[0] repeat_subtype = components[1] if len(components) > 1 else "" if repeat_type not in classification_nodes: - classification_nodes[repeat_type] = SoftmaxNode(repeat_type, parent=classification_tree, label_smoothing=label_smoothing, gamma=gamma, repeat_masker_name=repeat_type) + classification_nodes[repeat_type] = SoftmaxNode( + repeat_type, + parent=classification_tree, + label_smoothing=label_smoothing, + gamma=gamma, + repeat_masker_name=repeat_type, + ) repeat_type_node = classification_nodes[repeat_type] if repeat_subtype: - classification_nodes[repeat_name] = SoftmaxNode(repeat_subtype, parent=repeat_type_node, label_smoothing=label_smoothing, gamma=gamma, repeat_masker_name=repeat_name) - - node = classification_nodes[repeat_name] + classification_nodes[classification] = SoftmaxNode( + repeat_subtype, + parent=repeat_type_node, + label_smoothing=label_smoothing, + gamma=gamma, + repeat_masker_name=classification, + ) + + node = classification_nodes[classification] try: seqtree.add(accession, node, partition) @@ -68,13 +105,11 @@ def create_repeatmasker_seqtree(output:Path, repbase:Path, label_smoothing:float count += 1 - print("repbase,count,mapped,repeat_masker") + print("provided,count,mapped,repeat_masker") for classification,count in mapped_counter.most_common(): print(classification,count,1, mapping[classification], sep=",") for classification,count in not_mapped_counter.most_common(): print(classification,count,0, "", sep=",") - seqtree.save(output) - seqtree.classification_tree.render(print=1) - + return seqtree diff --git a/tests/test-data/custom.fna b/tests/test-data/custom.fna new file mode 100644 index 0000000..fe6abc5 --- /dev/null +++ b/tests/test-data/custom.fna @@ -0,0 +1,4 @@ +>I123#DNA/Academ +ACGT +>XXadsfadsf#LTR/Caulimovirus +AGCT \ No newline at end of file diff --git a/tests/test-data/custom.fna.gz b/tests/test-data/custom.fna.gz new file mode 100644 index 0000000..c3c7f38 Binary files /dev/null and b/tests/test-data/custom.fna.gz differ diff --git a/tests/test-data/repbase.ref b/tests/test-data/repbase.ref new file mode 100644 index 0000000..6cb5718 --- /dev/null +++ b/tests/test-data/repbase.ref @@ -0,0 +1,6 @@ +>IS905 DNA transposon Lactococcus lactis +ACGT +>BAGGINS1 Loa Drosophila melanogaster +ACGT +>SINE_DFSFDFS +ACGT \ No newline at end of file diff --git a/tests/test-data/simple.ref b/tests/test-data/simple.ref new file mode 100644 index 0000000..c243294 --- /dev/null +++ b/tests/test-data/simple.ref @@ -0,0 +1,4 @@ +>IJFDKSFF +ACGT +>GOOGLE +ACGT \ No newline at end of file diff --git a/tests/test_repeatmasker.py b/tests/test_repeatmasker.py new file mode 100644 index 0000000..82449fe --- /dev/null +++ b/tests/test_repeatmasker.py @@ -0,0 +1,109 @@ +from pathlib import Path +from Bio import SeqIO +from terrier.repeatmasker import get_verbatim_classification, create_repeatmasker_seqtree + + +TEST_DATA_DIR = Path(__file__).parent / "test-data" + +def test_get_verbatim_classification_repbase(): + fasta = TEST_DATA_DIR / "repbase.ref" + expected = [ + "DNA transposon", + "Loa", + "SINE", + ] + count = sum(1 for line in open(fasta) if line.startswith(">")) + assert count == len(expected) + with open(fasta) as f: + for record, exp in zip(SeqIO.parse(f, "fasta"), expected): + classification = get_verbatim_classification(fasta, record) + assert classification == exp + + +def test_get_verbatim_classification_simple(): + fasta = TEST_DATA_DIR / "simple.ref" + expected = [ + "Simple Repeat", + "Simple Repeat", + ] + count = sum(1 for line in open(fasta) if line.startswith(">")) + assert count == len(expected) + with open(fasta) as f: + for record, exp in zip(SeqIO.parse(f, "fasta"), expected): + classification = get_verbatim_classification(fasta, record) + assert classification == exp + + +def test_get_verbatim_classification_custom(): + expected = [ + "DNA/Academ", + "LTR/Caulimovirus", + ] + fasta = TEST_DATA_DIR / "custom.fna" + + count = sum(1 for line in open(fasta) if line.startswith(">")) + assert count == len(expected) + with open(fasta) as f: + for record, exp in zip(SeqIO.parse(f, "fasta"), expected): + classification = get_verbatim_classification(fasta, record) + assert classification == exp + + +def test_create_repeatmasker_seqtree_repbase(): + fasta = TEST_DATA_DIR / "repbase.ref" + seqtree = create_repeatmasker_seqtree([fasta]) + assert seqtree.classification_tree is not None + assert seqtree is not None + assert len(seqtree) == 3 + assert {"IS905", "BAGGINS1", "SINE_DFSFDFS"} == set(seqtree.keys()) + assert seqtree.classification_tree.render_equal( + """ + root + ├── DNA + ├── LINE + │ └── R1 + └── SINE + """ + ) + +def test_create_repeatmasker_seqtree_repbase_simple(): + seqtree = create_repeatmasker_seqtree([TEST_DATA_DIR / "repbase.ref", TEST_DATA_DIR / "simple.ref"]) + assert seqtree.classification_tree is not None + assert seqtree is not None + assert len(seqtree) == 5 + assert {"IS905", "BAGGINS1", "SINE_DFSFDFS", "IJFDKSFF", "GOOGLE"} == set(seqtree.keys()) + assert seqtree.classification_tree.render_equal( + """ + root + ├── DNA + ├── LINE + │ └── R1 + ├── SINE + └── Satellite + """ + ) + +def test_create_repeatmasker_seqtree_custom(): + seqtree = create_repeatmasker_seqtree([ + TEST_DATA_DIR / "repbase.ref", + TEST_DATA_DIR / "simple.ref", + TEST_DATA_DIR / "custom.fna", + ]) + assert seqtree.classification_tree is not None + assert seqtree is not None + assert len(seqtree) == 7 + assert {"IS905", "BAGGINS1", "SINE_DFSFDFS", "IJFDKSFF", "GOOGLE", "I123#DNA/Academ", "XXadsfadsf#LTR/Caulimovirus"} == set(seqtree.keys()) + assert seqtree.classification_tree.render_equal( + """ + root + ├── DNA + │ └── Academ + ├── LINE + │ └── R1 + ├── SINE + ├── Satellite + └── LTR + └── Caulimovirus + """ + ) + assert seqtree["I123#DNA/Academ"].node.name == "Academ"