diff --git a/.github/workflows/dev_docs.yaml b/.github/workflows/dev_docs.yaml
new file mode 100644
index 0000000..e09425f
--- /dev/null
+++ b/.github/workflows/dev_docs.yaml
@@ -0,0 +1,31 @@
+name: Build/Publish Develop Docs
+on:
+ push:
+ tags-ignore: ['**']
+ branches:
+ - dev
+permissions:
+ contents: write
+jobs:
+ deploy:
+ runs-on: ubuntu-latest
+ steps:
+ - uses: actions/checkout@v4
+ with:
+ fetch-depth: 0
+ - name: Configure Git Credentials
+ run: |
+ git config user.name github-actions[bot]
+ git config user.email 41898282+github-actions[bot]@users.noreply.github.com
+ - uses: actions/setup-python@v5
+ with:
+ python-version: 3.x
+ - run: echo "cache_id=$(date --utc '+%V')" >> $GITHUB_ENV
+ - uses: actions/cache@v4
+ with:
+ key: mkdocs-material-${{ env.cache_id }}
+ path: .cache
+ restore-keys: |
+ mkdocs-material-
+ - run: pip install mkdocs-material mike
+ - run: mike deploy --push Dev
diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml
index 5dd02c2..85a3784 100644
--- a/.github/workflows/main.yml
+++ b/.github/workflows/main.yml
@@ -1,9 +1,9 @@
-name: ci
+name: Build/Publish Latest Release Docs
+
on:
- push:
- branches:
- - main
- - docs_v2
+ release:
+ types: [published]
+
permissions:
contents: write
jobs:
@@ -11,6 +11,8 @@ jobs:
runs-on: ubuntu-latest
steps:
- uses: actions/checkout@v4
+ with:
+ fetch-depth: 0
- name: Configure Git Credentials
run: |
git config user.name github-actions[bot]
@@ -25,5 +27,5 @@ jobs:
path: .cache
restore-keys: |
mkdocs-material-
- - run: pip install mkdocs-material
- - run: mkdocs gh-deploy --force
+ - run: pip install mkdocs-material mike
+ - run: mike deploy --push --update-aliases ${{ github.event.release.tag_name }} latest
diff --git a/Dockerfile b/Dockerfile
new file mode 100644
index 0000000..7f6396f
--- /dev/null
+++ b/Dockerfile
@@ -0,0 +1,15 @@
+FROM mambaorg/micromamba:latest
+
+USER $MAMBA_USER
+
+COPY --chown=$MAMBA_USER:$MAMBA_USER docker_env.yml /tmp/env.yml
+RUN micromamba install -y -n base -f /tmp/env.yml && \
+ micromamba clean --all --yes
+
+COPY --chown=$MAMBA_USER:$MAMBA_USER . /tmp
+ARG MAMBA_DOCKERFILE_ACTIVATE=1
+RUN python -m pip install /tmp
+
+RUN rm -rf /tmp/*
+
+ENV PATH=${PATH}:/opt/conda/bin
\ No newline at end of file
diff --git a/README.md b/README.md
index c6b9c9b..ff4ad3c 100644
--- a/README.md
+++ b/README.md
@@ -1,95 +1,31 @@
# Small RNA Pipeline for the Hunt Lab
-This repository contains the code for the pipeline used to carry out sRNA-seq ananlysis in Vicky Hunt's lab. For usage and the tasks the pipeline can run, see the documentation, avalible in PDF, DOCX and HTML format in the `docs/` directory in the repository.
+This repository contains the code for the pipeline used to carry out sRNA-seq ananlysis in Vicky Hunt's lab. This is version 2 for version 1 see the `v1.x` branch.
-Note this software is very new, please report any bugs found to the GitHub bug tracker. Also note that it currently does not support GZipped input files, though support will be considered in the future if it is wanted enough.
+Documentation and information on how to install and run this pipeline can be found [here](https://vicky-hunt-lab.github.io/HuntLab-smallRNA/latest). Please report any bugs found or ask any other questions relating to the pipeline on the GitHub bug tracker.
## Installation
-Stable versions of the pipeline can be downloaded from releases on GitHub, unstable can be cloned using `git`.
+A containerised version of the pipeline is provided to allow for easy installation. Apptainer is recommended, as it allows for running with relative paths.
-It is recommended to use conda to install this pipeline. To perform the installation, ensure conda is installed on the machine, navigate to the directory containing the code and `environment.yml` file, then run:
+### Installing with Apptainer
+Create a SIF file with the code by running the following:
```
-$ conda env create -f environment.yml
+$ apptainer pull hlsmallrna.sif docker://kieranr51/huntlab-smallrna
```
-This will create a conda environment called huntlab-smallrna. This can then be activated each time you want to use the software with the command:
+To test it works try displaying the help message with the following command:
```
-$ conda activate huntlab-smallrna
+$ apptainer run hlsmallrna.sif hlsmallrna --help
```
-By default this provides the hlsmallrna executable. This provides all of the tasks listed in the documentation, each one can be run by typing hlsmallrna followed by the name and options for the task.
+If this displays the help message, you are good to go. All subsequent commands can be run with:
+```
+$ apptainer run -B $HOME:$HOME hlsmallrna.sif [COMMAND]
+```
-## Configuration File
-
-This pipeline can be configured with a "config.toml" file. This file is written in ["Tom's Obvious, Minimal Language"](https://toml.io/en/) (TOML) format, which is a common format for configeration. This file is split into multiple sections, with headers in square brackets. Valid keys are shown in the example file below. (Note lines starting with a # are comments like in python). An example is below, see the documentation for more comprihensive instructions.
-
-```toml
-# Configeration that doesn't fit anywhere else
-[general]
-# path to the directory to put the output in
-output_directory = "/path/to/output"
-
-# Number of threads to pass to any program that allows for multiprocessing
-threads = 8
-
-# This section contains subsitues for optional command line arguments
-[command]
-# for process command
-adapter = "ACGTTTTTT"
-front = "TTTTTACG"
-anywhere = ""
-cutoff = 25
-
-# for sort command
-min_length = 7
-max_length = 35
-
-# for unitas command
-refseq = [
- "/path/to/refseq1",
- "/path/to/refseq2"
-]
-species = "x"
-
-# for targetid command
-target_files = [
- "/path/to/targetseq1",
- "/path/to/targetseq2"
-]
-min_seq_length = 9
-num_mismatches = 0
-
-# Section containing configeration for each of the exteranl programs called while in use
-[cli-tools]
-[cli-tools.trim]
-trim_pass_threads = true
-path_to_trim = "cutadapt"
-trim_type = ""
-trim_params = []
-
-[cli-tools.fastqc]
-fastqc_pass_threads = true
-path_to_fastqc = "fastqc"
-fastqc_params = []
-
-[cli-tools.unitas]
-unitas_pass_threads = true
-path_to_unitas = "unitas"
-unitas_params = []
-
-[cli-tools.bowtie2]
-bowtie2_pass_threads = true
-path_to_bowtie2 = "bowtie2"
-bowtie2_params = []
-path_to_bowtie2_build = "bowtie2-build"
-bowtie2_build_params = []
-
-[cli-tools.samtools]
-path_to_samtools = "samtools"
-samtools_view_params = ["-h", "-F", "256", "-F", "4"]
-samtools_fastq_params = []
-```
\ No newline at end of file
+If this is not possible, other methods of installation are covered in the documentation at:
+[https://vicky-hunt-lab.github.io/HuntLab-smallRNA/latest/installation/](https://vicky-hunt-lab.github.io/HuntLab-smallRNA/latest/installation/)
\ No newline at end of file
diff --git a/bin/overlap_ss.sh b/bin/overlap_ss.sh
deleted file mode 100644
index 72ac9b6..0000000
--- a/bin/overlap_ss.sh
+++ /dev/null
@@ -1,180 +0,0 @@
-#!/bin/sh
-# Copyright 2022 Vicky Hunt Lab Members
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-# http://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software
-# distributed under the License is distributed on an "AS IS" BASIS,
-# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-# See the License for the specific language governing permissions and
-# limitations under the License.
-echo "Starting..."
-
-# echo "Enter Genome file name: "
-# read genome
-# #"celgenome.fa"
-
-# echo "First sequence file (e.g. 22nt sequences): "
-# read ms
-# #"Cel_M_col.fa"
-
-# echo "Second sequence file (e.g. 27nt sequences): "
-# read ps
-# #"Allseq_Cel_male_LF.fa"
-
-genome=$1
-ms=$2
-ps=$3
-longestseq=$4
-
-# Index the genome of interest
-bowtie2-build $genome genome
-
-# Map target sequences and passenger sequences to the genome using Bowtie2
-
-bowtie2 -x genome -U $ms -f -a -N 0 --no-1mm-upfront --norc > Main.sam
-bowtie2 -x genome -U $ps -f -a -N 0 --no-1mm-upfront --norc > Pass.sam
-
-echo "bowtie2 complete"
-
-# revcomp_rna -i $ms -o Rev_main.fa
-# revcomp_rna -i $ps -o Rev_pass.fa
-
-# grep -v ">" Rev_main.fa > ms_nohead.txt
-# grep -v ">" Rev_pass.fa > ps_nohead.txt
-
-# echo "Reverse complement files created"
-
-# Convert sam output to bam output with samtools
-
-for f1 in `ls -1 *.sam`
-do
- f2="$(basename $f1 .sam)"
- samtools view -S -b $f1 > $f2.bam
-done
-echo "sam to bam files complete"
-
-# samtools view -h Main.bam | grep -vf ms_nohead.txt | samtools view -bS -o Main_filter.bam -
-# samtools view -h Pass.bam | grep -vf ps_nohead.txt | samtools view -bS -o Pass_filter.bam -
-
-cp Main.bam Main_filter.bam
-cp Pass.bam Pass_filter.bam
-
-mkdir Originals
-mv Main.bam Originals
-mv Pass.bam Originals
-
-# Sort the bam output file with samtools
-for f1 in `ls -1 *.bam`
-do
- f2="$(basename $f1 .bam)"
- samtools sort $f1 -o $f2.sort.bam
-done
-echo "bam files sorted"
-
-# Index the sorted bam file using samtools
-for f1 in `ls -1 *.sort.bam`
-do
- samtools index $f1
-done
-echo "bam file index complete"
-
-# Convert indexed bam files to bed files using bedtools Bamtobed
-for f1 in `ls -1 *.sort.bam`
-do
- f2="$(basename $f1 .sort.bam)"
- bedtools bamtobed -i $f1 > $f2.bed
-done
-echo "bed files created"
-
-bedtools intersect -s -a Main_filter.bed -b Pass_filter.bed -wo > intersectOUTSS.txt
-
-echo "same strand intersect file created"
-
-awk '{print $1,$2,$3,$4,$7,$8,$9,$10,$13}' intersectOUTSS.txt > cleanupTEST.txt
-
-echo "File clean up complete"
-
-# Caculate length columns.
-# Column 9 will be column 3 minus column 2. Column 10 will be column 7 minus column 6
-
-awk '{$10 = $3-$2; print}' cleanupTEST.txt > length1.txt
-awk '{$11 = $7-$6; print}' length1.txt > length2.txt
-
-echo "Lengths calculated"
-
-# Calculate begining and end overhangs.
-# Column 11 will be column 6 minus column 2 - Overhang beginning
-awk '{$12 = $2-$6; print}' length2.txt > overhang1.txt
-
-# Column 12 will be column 7 minus column 3 - Overhang end
-awk '{$13 = $3-$7; print}' overhang1.txt > overhang2.txt
-
-echo "Overhangs calculated"
-
-
-awk ' $12 < 0 { print $0 ; }' overhang2.txt > beg_less0.txt
-awk ' $12 > 0 { print $0 ; }' overhang2.txt > beg_more0.txt
-
-cat beg_less0.txt beg_more0.txt > beg_file.txt
-
-awk ' $13 < 0 { print $0 ; }' overhang2.txt > end_less0.txt
-awk ' $13 > 0 { print $0 ; }' overhang2.txt > end_more0.txt
-
-cat end_less0.txt end_more0.txt > end_file.txt
-
-cat beg_file.txt end_file.txt > No_blunts.txt
-echo "Blunts removed"
-
-sort -u No_blunts.txt > No_Dups.txt
-
-echo "Duplicates removed"
-
-awk '{print $1,$2,$3,$4}' No_Dups.txt > Cut_no_dup.txt
-sort -u Cut_no_dup.txt > No_Dups_sorted.txt
-awk '{print $4}' No_Dups_sorted.txt > Read_num.txt
-sort -u Read_num.txt > Nodup_read.txt
-
-echo "Read_Number, Occurences" > Read_number_out.txt
-for read in `cat Nodup_read.txt`
-do
-awk -v num=$read '$1 == num { print $0 ; }' Read_num.txt > ReadTemp.txt
-count=`wc -l ReadTemp.txt | cut -d' ' -f1`
-echo "$read,$count" >> Read_number_out.txt
-rm ReadTemp.txt
-done
-
-awk '{print $1,$2,$3,$5,$6,$7,$9}' No_Dups.txt > OverlapOUT.txt
-
-echo "Overlap, Occurrences" > OverlapLengths.txt
-for num in $(seq 1 $longestseq);
-do
-awk -v num=$num '$7 == num { print $0 ; }' OverlapOUT.txt > OverTemp.txt
-count=`wc -l OverTemp.txt | cut -d' ' -f1`
-echo "$num,$count" >> OverlapLengths.txt
-done
-rm OverTemp.txt
-
-echo "Common Overlap Lengths counted"
-
-awk '{print $3}' OverlapOUT.txt > end.txt
-sort -u end.txt > unique.txt
-
-echo "Sequence1End, matches" > Distribution.txt
-for coord in `cat unique.txt`
-do
-awk -v num=$coord '$3 == num { print $0 ; }' OverlapOUT.txt > MatchTemp.txt
-count=`wc -l MatchTemp.txt | cut -d' ' -f1`
-echo "$coord,$count" >> Distribution.txt
-rm MatchTemp.txt
-done
-
-tar -cf Ouput_zip.tar unique.txt genome* Main* Pass* beg* end* length* overhang* cleanupTEST.txt No* ps* Rev* ms*
-rm unique.txt genome* Main* Pass* beg* end* length* overhang* cleanup.txt No* ps* Rev* ms*
-
-echo "Distribution of sequences counted"
-echo "Complete"
\ No newline at end of file
diff --git a/config.toml b/config.toml
deleted file mode 100644
index b5e5a11..0000000
--- a/config.toml
+++ /dev/null
@@ -1,28 +0,0 @@
-# Configeration file I have been using for testing
-[general]
-threads = 8
-output_directory = "./output"
-
-[command]
-target_files = [
- "/Users/kieran/Documents/BathWorkDocs/RNAPipelineProject/Reference files/miRNA_sratti_hairpin_incNovel.fasta",
- "/Users/kieran/Documents/BathWorkDocs/RNAPipelineProject/Reference files/rRNA_srattiWB.fasta",
- "/Users/kieran/Documents/BathWorkDocs/RNAPipelineProject/Reference files/strongyloides_ratti.PRJEB125.WBPS12.CDS_transcripts.fasta",
- "/Users/kieran/Documents/BathWorkDocs/RNAPipelineProject/Reference files/TE_sequences.fasta",
- "/Users/kieran/Documents/BathWorkDocs/RNAPipelineProject/Reference files/tRNA_sratti_WB.fasta"
-]
-
-[cli-tools]
-[cli-tools.trim]
-trim_pass_threads = true
-trim_type = ""
-trim_parms = []
-
-[cli-tools.fastqc]
-fastqc_pass_threads = true
-path_to_fastqc = "/Applications/FastQC.app/Contents/MacOS/fastqc"
-fastqc_params = []
-
-[cli-tools.unitas]
-path_to_unitas = ""
-unitas_params = []
\ No newline at end of file
diff --git a/docker_env.yml b/docker_env.yml
new file mode 100644
index 0000000..0824393
--- /dev/null
+++ b/docker_env.yml
@@ -0,0 +1,23 @@
+name: base
+channels:
+ - bioconda
+ - conda-forge
+dependencies:
+ - python
+ - pyyaml
+ - pysam
+ - biopython
+ - gffutils
+ - matplotlib
+ - cutadapt
+ - unitas
+ - perl-archive-extract
+ - bioconda/label/main::perl-lwp-simple
+ - bowtie2
+ - samtools
+ - bedtools
+ - eggnog-mapper
+ - scipy
+ - r
+ - bioconductor-topgo
+ - pip
\ No newline at end of file
diff --git a/docs/Documentation.docx b/docs/Documentation.docx
deleted file mode 100644
index 35491db..0000000
Binary files a/docs/Documentation.docx and /dev/null differ
diff --git a/docs/Documentation.pdf b/docs/Documentation.pdf
deleted file mode 100644
index 7eb9d7a..0000000
Binary files a/docs/Documentation.pdf and /dev/null differ
diff --git a/docs/cli-reference/build-coord-files.md b/docs/cli-reference/build-coord-files.md
new file mode 100644
index 0000000..0b74ff2
--- /dev/null
+++ b/docs/cli-reference/build-coord-files.md
@@ -0,0 +1,30 @@
+# Build Coord Files
+
+```
+usage: build_coord_files [-h] [-a] [-q] [-s] [-b] [-g] [-r] [-c]
+ [--feature FEATURE] [--genome-coords GENOME_COORDS]
+ [--band-width BAND_WIDTH] [-o OUTPUT]
+ input_file
+
+Convert common bioinformatics file formats into coordinate file to plot
+
+positional arguments:
+ input_file File to convert, attempts to autodetect type
+
+optional arguments:
+ -h, --help show this help message and exit
+ -a, --fasta Treat input as a FASTA file
+ -q, --fastq Treat input as a FASTQ file
+ -s, --sam Treat input as a SAM file
+ -b, --bam Treat input as a BAM file
+ -g, --gff Treat input as a GFF file
+ -r, --rm-fa-out Treat input as a RepeatMasker .fa.out file
+ -c, --scaffold-aware Merge scaffolds into one chromosome
+ --feature FEATURE Select a gff feature to use
+ --genome-coords GENOME_COORDS
+ File containing the coordinates of the genome
+ --band-width BAND_WIDTH
+ Size of the bands showing change in scaffolds
+ -o OUTPUT, --output OUTPUT
+ File to output to
+```
\ No newline at end of file
diff --git a/docs/cli-reference/extractnc.md b/docs/cli-reference/extractnc.md
new file mode 100644
index 0000000..ed90563
--- /dev/null
+++ b/docs/cli-reference/extractnc.md
@@ -0,0 +1,28 @@
+# ExtractNC
+
+```
+usage: extract_nc [-h] [-o OUTPUT] genome gff_file
+
+positional arguments:
+ genome FASTA containing the genome to extract from
+ gff_file GFF file containing annotations of CDS and mRNA
+ regions
+
+options:
+ -h, --help show this help message and exit
+ -o OUTPUT, --output OUTPUT
+ FASTA file to write output to
+```
+
+Input files:
+
+- `genome` - file containing the whole genome of the species of interest
+
+- `gff_file` - GFF3 file containing annotations at least for the mRNA and coding region (labelled CDS) for the genome
+
+Output files:
+
+- `noncoding.fasta` - FASTA file containing only the transcribed noncoding regions of the DNA
+
+!!! warning
+ The method used here will not work if there a coding region is labelled outside a single mRNA region. The program does a check before running to make sure this is true and throws an exception before trying to run if so, to prevent nonsense results from being produced.
diff --git a/docs/cli-reference/index.md b/docs/cli-reference/index.md
new file mode 100644
index 0000000..f8f91fd
--- /dev/null
+++ b/docs/cli-reference/index.md
@@ -0,0 +1,27 @@
+# CLI Reference
+References of how to use the command line interface of the program and notes on some ways particular commands work that could trip you up.
+
+General help message is as follows:
+
+```
+usage: hlsmallrna [-h] [-o OUTPUT] [-t THREADS] [-k] [-v] config_file
+
+Pipeline to process small RNAs - version 2
+
+positional arguments:
+ config_file Path to YAML config file
+
+options:
+ -h, --help show this help message and exit
+ -o OUTPUT, --output OUTPUT
+ Directory to write output to (will be created if it
+ doesn't exist, default: hlsmallrna_output)
+ -t THREADS, --threads THREADS
+ Number of threads to run the pipeline on (default: 4)
+ -k, --keep-files If set, keep the intermediate files, to help debug the
+ application
+ -v, --verbose If set, print the output for the intermediate commands
+```
+
+!!! note "The Config File"
+ If you are looking for instructions on writing the config file, use [the relevent page of the walkthrough](../walkthrough/config.md) to help.
\ No newline at end of file
diff --git a/docs/cli-reference/label-for-unitas.md b/docs/cli-reference/label-for-unitas.md
new file mode 100644
index 0000000..ca9411d
--- /dev/null
+++ b/docs/cli-reference/label-for-unitas.md
@@ -0,0 +1,17 @@
+# Label for Unitas
+
+```
+usage: label_for_unitas [-h] [-o OUTPUT] label file_path
+
+Program to automatically prepend unitas labels to all the sequences in a FASTA
+file
+
+positional arguments:
+ label Label to prepend to the sequences in a fasta file
+ file_path Path to the fasta file to label
+
+optional arguments:
+ -h, --help show this help message and exit
+ -o OUTPUT, --output OUTPUT
+ File to output to, defaults to labelled.fasta
+```
\ No newline at end of file
diff --git a/docs/cli-reference/overlapss.md b/docs/cli-reference/overlapss.md
new file mode 100644
index 0000000..f1039e1
--- /dev/null
+++ b/docs/cli-reference/overlapss.md
@@ -0,0 +1,16 @@
+# Overlap SS
+
+```
+usage: overlap_ss [-h] [-q] genome rna_file_1 rna_file_2
+
+Looks for overlaps of two classes of RNA on the same strand of DNA
+
+positional arguments:
+ genome FASTA file containing the genome of the organism
+ rna_file_1 File containing the RNA to use as a base
+ rna_file_2 File containing the RNA to look for overlaps with
+
+optional arguments:
+ -h, --help show this help message and exit
+ -q, --quiet Suppress output from intermediate commands
+```
\ No newline at end of file
diff --git a/docs/html/Documentation.html b/docs/html/Documentation.html
deleted file mode 100644
index d181522..0000000
--- a/docs/html/Documentation.html
+++ /dev/null
@@ -1 +0,0 @@
-Hunt Lab Small RNA Pipeline Docs Version 1.1
Hunt Lab Small RNA Pipeline Docs
Kieran Reynolds
This tool is designed to automate a lot of the common bioinformatics performed on small RNA by the Hunt Lab. It is deliberately designed to be highly configurable so it can be easily applied to a wide range of scenarios.
Installation 2
Tasks 2
Example Usage 3
Process 3
Sort 4
Unitas 5
Extract Non-Coding 6
Target Identification 6
All 7
Plotting Genome Location 7
Basic Plot 7
Dealing with Multiple Scaffolds 9
Further Processing Options 11
Signature Identification 11
AT/GC Richness 11
Upstream motif (e.g. piRNA) 11
Ping Pong signature 11
Dicer Signature (StepRNA) 12
Same Strand Overlap 13
miRNA 13
Notes on Settings 13
Configuration file 13
Manually setting parameters of internal commands 14
Suppressing excess output 14
Troubleshooting 15
CLI Reference 16
Process 16
Sort 17
ExtractNC 18
Unitas 18
TargetID 19
Label for Unitas 20
Build Coord Files 20
Overlap SS 21
Installation
You can obtain the code from https://github.com/kieranr51/HuntLab-smallRNA. To download, click the code button and copy the URL. You can then download it to a computer with git installed on it by opening the command line and running:
$ git clone <copied_url>
It is recommended to use conda to install this pipeline. To perform the installation, ensure conda is installed on the machine, navigate to the directory containing the code and environment.yml file, then run:
$ conda env create -f environment.yml
This will create a conda environment called huntlab-smallrna. This can then be activated each time you want to use the software with the command:
$ conda activate huntlab-smallrna
By default this provides the hlsmallrna executable. This provides all of the tasks listed below, each one can be run by typing hlsmallrna followed by the name and options for the task.
Tasks
Currently the following tasks have been implemented in it:
- Process - Trim adapters off raw RNA-Seq data, run a quality control report and trim of any bases that fall below a specified threshold from all the RNA
- Sort - Align the RNA against a given reference, removing any sequences that fail to align, then split them into files based on sequence length to allow for easy further processing. Can also remove sequences below a minimum length and above a maximum length if desired
- ExtractNC - Using a genome and GFF file containing annotations, extract the region that is transcribed, but not a coding region
- Unitas - Run a directory of RNA files (e.g. output of sort step) though unitas to classify the types of RNA, then combine the summaries into one spreadsheet for easy viewing and manipulation
- TargetID - Identify which of a set of potential targets a set of small RNA are targeting by reverse complimenting them and aligning them to the potential targets
More in depth documentation on each of these steps, starting with examples and suggestions for further processing. Further down there is a CLI Reference, that contains general help for each command, including what files are required and what outputs you get from each step.
Example Usage
To use this pipeline, you will need at least a FASTQ file containing raw small RNA-Seq reads. Then depending on what you want to achieve, you can run different steps of the pipeline. If you ever get stuck and need to see the options of a particular stage, all parts have their own help message, that can be seen by passing --help after the name of the stage.
Before starting though, you will need to create a configuration file that allows you to specify some parameters in advance, so you don’t have to keep retyping them. By default, it looks for these in a file in your current directory called config.toml, though alternatives can be specified by passing the name to the -C argument before the name of the stage to run. More details on this can be found in the ‘Configuration file’ section below, but for now setup a simple file to specify an output directory of output/ by putting the following in config.toml:
[general]
output_directory = "./output"
In addition, check that all of the commands can be run by just typing their name in the command line, as this is how the program tries to use them by default. If they can’t, you can set a custom path by setting the relevant path_to_ variable in the configuration file. For example if unitas is run with unitas.pl instead of just unitas add the following to the bottom of the config file:
[cli-tools]
[cli-tools.unitas]
path_to_unitas = "unitas.pl"
Process
The first stage is Process. This is for if you have raw data that hasn’t had it’s adapters trimmed and been quality filtered, so if your data has had this done to it already, this step can be skipped. To run you need a FASTQ file containing raw small RNA and the sequence and end of any adapters attached during sequencing.
Once you have this data, to include the trimming in the run, you need to set one of three flags with the adapter sequence, depending on the position of the particular adapter.
If your adapter is on the 3’ end only, use -a; if it is from the 5’ end only, use -g and if it could be on either end use -b. If you don’t need adapters trimmed, just don’t specify any of them and the step will be skipped. For example, if your small RNAs are in smallRNA.fastq and you want to trim the adapter sequence AGCATA from the 3’ end only, you would run:
$ hlsmallrna process -a AGCATA smallRNA.fastq
In addition, this step runs FastQC to produce a report on the quality of the small RNA and automatically cuts anything with a lower quartile of quality below a cutoff. By default, this is set to 20, but can be changed by the user using the -c flag. If you do not want this, you can set -c 0 and the FastQC step will be skipped for efficiency. For example if you want to change the cutoff for the last command to 5, you could run:
$ hlsmallrna process -c 5 -a AGCATA smallRNA.fastq
This produces a number of files in the output directory, importantly cut_sequences.fastq has the sequences that have adapters trimmed and low quality parts removed. (See CLI Reference section for information on other output for all of the commands)
Sort
The second stage is Sort. This stage aligns the small RNAs to the genome of the organism and filters out any that don’t align, to remove contamination. Then it splits the remaining reads into FASTQ files by length, to allow for easier classification.
To run this, you need a FASTQ file of processed small RNA and a FASTA file containing the genome of the organism of interest. For example, if your small RNA are in smallRNA.fastq and your genome is in genome.fasta you can run:
$ hlsmallrna sort smallRNA.fastq genome.fasta
If you are only interested in a subset of lengths of small RNA, this step can be set to restrict the lengths it outputs with -l for minimum and -x for maximum. For example, if you only want small RNA with lengths between 12 and 30, you could run:
$ hlsmallrna sort -l 12 -x 30 smallRNA.fastq genome.fasta
By default this aligns the small RNA to the genome with bowtie2’s default settings, which allow different amounts of mismatches depending on how long the RNA is. Details of this can be found in the score configuration options in the manual, that can be found at: https://bowtie-bio.sourceforge.net/bowtie2/manual.shtml. To deal with cases where we may want to do a more stringent run, the -m or --ref-mismatches argument can be supplied with the number of mismatches for bowtie2 to allow. For example, if you want to perform an operation using the same files as above, but only include small RNA that map exactly to the genome, you could run:
$ hlsmallrna sort -m 0 smallRNA.fastq genome.fasta
This command produces a few useful files in the output directory:
- binned_rna/ - directory containing one fastq file for each length of small RNA
- rna_length_report.csv - table showing a summary of the RNAs by length and first base
- baseplot.png - plot of the length and first base of the RNAs (this is sometimes wrong, replot from the csv values for now, will fix matplotlib code soon)
Unitas
The third stage is Unitas. This runs unitas on a directory containing a range of lengths of small RNA to classify the origin of the RNA. To do this, you will need files containing reference sets of the genome targets you wish to find, with labels in their IDs and a directory containing small RNA of different lengths, similar to the output of sort.
You can set the reference file using the -r command line argument, but I think it is easier to set them in the configuration file, as it makes the command shorter to type. If your reference files are in the same directory you are running in and called ref1.fasta, ref2.fasta and ref3.fasta, you can add the following to the config file:
[command]
refseq = [
"ref1.fasta",
"ref2.fasta",
"ref3.fasta"
]
Note that all of these reference sequences should be labelled with the name of the category they belong to and the | (bar or pipe character) before the start of their ID. If this hasn’t been done and all of the sequences in a FASTA file should be labelled the same, you can use the label_for_unitas script that was installed with the small RNA pipeline. For example if you want to label the sequences in unlabelledTEs.fasta with the label TE, with the result produced in teRef.fasta, you could run:
$ label_for_unitas "TE" unlabelledTEs.fasta -o teRef.fasta
Once the correct reference file are set, if your small RNA, split by length, are in output/binned_rna, you can run:
$ hlsmallrna unitas output/binned_rna
This produces a useful summary and graph in the output directory:
- unitasGraph.png - graph showing the categories allocated by unitas against the length of the small RNA (this is sometimes wrong, replot from the csv values for now, will fix matplotlib code soon)
- unitas_graph_data.csv - summary of the raw numbers produced from the unitas runs, used to produce the graph
If you want to get small RNA that are derived from genes, you will need to obtain both of the mRNA and CDS regions and pass them to speciled arguments in the pipeline so they can be combined. This can either be done using the config file as followers or with the corresponding command line arguments (--cds and --unspliced-transcriptome):
[command]
refseq = [
"ref1.fasta",
"ref2.fasta",
"ref3.fasta"
]
cds = “cds_sequences.fasta”
unspliced_transcriptome = “unspliced_transcriptome.fasta”
Extract Non-Coding
Sometimes you want to see if a small RNA comes from a non-coding region of the transcriptome, so the pipeline contains a utility command to automatically extract this using a genome and some annotations. At minimum, these annotations need to contain a transcribed region labelled with ‘mRNA’ and a coding region within the transcribed region labelled ‘CDS’. Most of the ones on Wormbase Parasite are fine for this purpose.
As an example, to extract the regions from genome.fasta using the annotations in annotations.gff you could run:
$ fungi
This will produce a FASTA file containing the non-coding regions, complete with unitas labels applied, called noncoding.fasta in the output directory.
Target Identification
The fourth and final main step is Target Identification. This takes a set of small RNA and a set of potential targets and produces a list of genome features that are targeted by each small RNA.
As with unitas, you can set target files using the command line option -t but it is easier to add them to the configuration file. For example if your target files are target1.fasta, target2.fasta and target3.fasta, they can be set by adding the following in the section below [command]:
target_files = [
"target1.fasta",
"target2.fasta",
"target3.fasta"
]
Then, if your small RNA are in smallRNA.fastq you can run the following:
$ hlsmallrna targetid smallRNA.fastq
Make sure the sequences in both smallRNA.fastq and each of the reference files have unique sequence ids, otherwise the command will fail.
If you are working with very short sequences (<5 bases), you will need to reduce the minimum length with the -m option. Conversely, if you know your minimum length is longer than this, this step can be sped up by increasing it. For example, if you know the shortest sequence you are dealing with is 10 bases long, you could run:
$ hlsmallrna targetid -m 10 smallRNA.fastq
Though this is not necessary and the program will produce a correct result for sequences longer than 5 without it, so only add if you need the extra speed.
This produces rna_target_list.csv which contains a list of small RNA target pairs and some general information about each pair’s alignment.
If you are dealing with a small RNA that doesn’t need perfect complementarity, you can add the --num-mismatches option to set the number of mismatches the alignment should allow. For example, if you want to allow up to four mismatches, you could run:
$ hlsmallrna targetid --num-mismatches 4 smallRNA.fastq
All
Since the first three steps are often run one after the other, a shorthand command all is provided for convenience. This takes all of the arguments from the Process, Sort and Unitas steps, then runs them, piping the result from one into the next. All files produced by those steps are then produced in the output directory.
For example if you have small RNAs in smallRNA.fastq, a genome in genome.fasta and refseq set in the config file and you want to trim the adapter sequence ACATA from both ends, but not run quality filtering you could run:
$ hlsmallrna all -c 0 -b ACATA smallRNA.fastq genome.fasta
Plotting Genome Location
Basic Plot
A common task you might want to do is plotting the location of things within the genome, for example, genome features you are interested in or the positions of small RNAs you have aligned to the genome. To make this easy, a script is included in this environment called build_coord_files. This automatically generates tsv files for use with the R package chromPlot. This section explains how you can generate a plot with it.
Firstly, you need at least two files: a full genome in a FASTA or FASTQ file and a set of locations you are interested in, either as an alignment output (in SAM or BAM format) or as feature annotations in GFF format. To start produce the genome coordinate file with the command: (assuming genome.fasta is your genome file)
$ build_coord_files -o genome_coord.tsv genome.fasta
By default, the script will figure out what type of output to produce by the file extension. But if it can’t, the user will need to provide a flag to tell it what type of input you are using (in this case --fasta). Next we produce the coordinates file for the item of interest, assuming it is in alignment.sam you can run:
$ build_coord_files -o alignment_coord.tsv alignment.sam
This will produce two files, one for coordinates on the sense strand and one for coordinates on the antisense strand. Finally, we can use these files to plot the result. If it isn’t already, install chromPlot with bioconducter using the following R command:
> BiocManager::install("chromPlot")
Then you should be able to create a plot using the following R code (of two example chromosomes called ‘SRAE_chr1’ and ‘SRAE_chr2’):
chrom_coords <- read.csv("/path/to/genome_coord.tsv", sep = "\t")
sense_pos <- read.csv("/path/to/sense_alignment_coord.tsv", sep = "\t")
antisense_pos <- read.csv("/path/to/antisense_alignment_coord.tsv", sep = "\t")
library(chromPlot)
chromPlot(
gaps = chrom_coords, annot1 = sense_pos, annot2 = antisense_pos, plotRndchr = T,
chr = c("SRAE_chr1", "SRAE_chr2"), bin = 50000,
chrSide = c(-1, 1, -1, -1, 1, -1, -1, 1)
)
The script breaks down as follows:
- First three lines load in the data we just produced with build_coord_files
- Following line loads the chromPlot library into R
- The final line plots the data with chromPlot:
- gaps defines what each chromosome should look like
- annot1 and annot2 define the annotation data to show
- plotRndchr disables the assumption in chromPlot that the chromosomes are labelled with numbers
- chr selects which chromosomes to plot by specifying an array of names
- bin sets the bin size to group the features into when plotting
- chrSide defines which side of the chromosome to plot each feature, the only change we have made from the default is setting the second item to 1 as opposed to -1. This ensures annot2 is plotted on the opposite side to annot1
The resulting plot will look something like:
Dealing with Multiple Scaffolds
In more complex cases, one or more of your chromosomes may be split into multiple scaffolds. Build_coord_files has a scaffold aware mode that can automatically try to merge these, but it will only work if your scaffold’s IDs are in the form of <chromosome_name>_scaffold<scaffold_number> e.g. SRAE_chrX_scaffold1. Once you have this, you can generate the data files with the -c flag set, as follows:
$ build_coord_files -c -o genome_coord.tsv genome.fasta
$ build_coord_files -c --genome-coords genome_coord.tsv -o alignment_coord.tsv alignment.sam
Note you need to provide the genome coords file when generating files for what you are plotting, this allows for a correction of coordinates to be made from where scaffolds are merged. In addition to the regular files, this produces scaffold_info_genome_coord.tsv, that can be used to add scaffold dividers to your plots. To do this, the following R script can be run:
chrom_coords <- read.csv("/path/to/genome_coord.tsv", sep = "\t")
sense_pos <- read.csv("/path/to/sense_alignment_coord.tsv", sep = "\t")
antisense_pos <- read.csv("/path/to/antisense_alignment_coord.tsv", sep = "\t")
scaffolds <- read.csv("/path/to/scaffold_info_genome_coord.tsv", sep = "\t")
library(chromPlot)
chromPlot(
gaps = chrom_coords, annot1 = sense_pos, annot2 = antisense_pos, plotRndchr = T,
chr = c("SRAE_chr1", "SRAE_chr2"), bin = 50000, bands = scaffolds,
chrSide = c(-1, 1, -1, -1, 1, -1, -1, 1)
)
Compared to the last script, this adds:
- The fourth line, that loads in the scaffold data
- The bands argument, that draws red lines at the end of each scaffold
Note that, by default the bands are 20,000 bp wide, as that looks good on the plot. At that size they may not appear if you make the bins too much bigger. Therefore, if you plan to make the bins bigger, pass a larger value to the --band-width flag when creating the genome file. That makes the bands bigger, which results in a more visible line on the plot.
A plot produced by this method will look something like:
Further Processing Options
Signature Identification
AT/GC Richness
Nucleotide richness and the presence of any conserved motifs is identified using WebLogo (weblogo.threeplusone.com). Sequences need to be the same length and can be either uploaded using a fasta file or can be pasted in.
Upstream motif (e.g. piRNA)
map to genome with bowtie2, convert sam to bam with samtools, convert bamToBed with bedtools
$ bowtie2 -x genome.fasta -f -U sequences.fasta -S sequences_mapped.sam
$ samtools view sequences_mapped.sam -o sequences_mapped.bam
$ bamToBed -i sequences.bam > sequences.bed
extract flanking sequences with flankBed in bedtools (need to use bedtools format genome file - generate with faidx and cut)
$ samtools faidx genome.fasta
$ cut -f 1,2 genome.fa.fai > chrom.sizes
$ flankBed -i sequences.bed -g chrom.sizes -b 100 > sequences_100nt_flank.bed
or 60 upstream only
$ flankBed -i sequences.bed -g chrom.sizes -l 60 -r 0 > sequnences_60nt.bed
extract fasta sequences from bed fastFromBed with bedtools
$ fastaFromBed -fi genome.fasta -bed sequences_100nt_Flank.bed -fo sequences_100nt_Flank_seqs.fasta
Ping Pong signature
To find ping-pong signature, Unitas with the option -pp can be used to find 5’ overlaps of mapped sequence reads to the genome and calculate a Z-score for the enrichment of 10 bp overlaps.
$ perl unitas.pl -pp -input map.file (SAM or ELAND3 format) –species x –refseq reference.fasta
Dicer Signature (StepRNA)
To install stepRNA, see the documentation on the GitHub page (https://github.com/bmm514/stepRNA).
stepRNA is able to identify a Dicer processing signature from a small RNA sequence dataset. It requires two FASTA files as input that have had their adapters already trimmed from the ends.
$ stepRNA -r/--reference REFERENCE.fa -q/--reads READS.fa
- READS.fa should contain all of the small RNA sequencing reads.
- REFERENCE.fa is usually filtered to investigate small RNAs of interest e.g. 26G sRNAs.
Identical reads in the READ and REFERENCE files can also be removed before searching for a Dicer signature, if desired, by using -e/--remove_exact.
Importantly the FASTA headers must be unique - this can be done by stepRNA with -u/--make_unique.
stepRNA will then generate:
- BAM alignment files (with different combinations of overhang lengths)
- Overhang length CSV
- Passenger Read Length CSV
These can then be used to plot the distribution and/or further explore the reads that have aligned.
Note: If the sRNA lengths in the read file are all the same i.e. 21 nt long, stepRNA will not be able to find overhangs due to the scoring system. In this case the -m/--min_score should be set to an appropriate value; we recommend sRNA length - 7nt (e.g. 21 - 6 = 15).
$ stepRNA --help to bring up the command line help, a description of the methods can be found in the stepRNA publication (add_when_available)
stepRNA [-h] -r REFERENCE -q READS [-n NAME] [-d DIRECTORY]
[-m MIN_SCORE] [-e] [-u] [-j] [-V]
Align a reference RNA file to read sequences. Output will be a set of CSV files containing information about the length of the reads, number of reads aligned to a reference sequence and the length of overhangs of the alignment.
Reference RNA file will be automatically indexed
Optional Arguments:
-h, --help show this help message and exit
-n NAME, --name NAME Prefix for the output files
-d DIRECTORY, --directory DIRECTORY
Directory to store the output files
-m MIN_SCORE, --min_score MIN_SCORE
Minimum score to accept, default is the shortest read
length
Required Arguments:
-r REFERENCE, --reference REFERENCE
Path to the reference sequences
-q READS, --reads READS
Path to the read sequences
Flags:
-e, --remove_exact Remove exact read matches to the reference sequence
-u, --make_unique Make FASTA headers unique in reference and reads i.e.
>Read_1 >Read_2
-j, --write_json Write count dictionaries to a JSON file
-V, --version Print version number then exit.
Same Strand Overlap
If you wish to detect same strand overlaps, you can use the overlap_ss script that is installed in the same environment as this pipeline. This uses a number of bash programs to create files that show where two sets of small RNA that appear on the same strand overlap.
To use this script you need two sets of small RNA as FASTA or FASTQ files, one for the base and one for the overlap, plus the genome of the species of interest. You can then run the script in the following way:
$ overlap_ss genome.fasta smallRNA1.fastq smallRNA2.fastq
This script will then create output in the output directory specified in config.toml under the directory samestrand_overlap/.
miRNA
Notes on Settings
Configuration file
For ease of use this program requires a configuration file in TOML format to run successfully. At minimum it can be a blank file, but I strongly suggest setting at least the output_directory setting to be named something more useful than output, the default setting.
By default the program tries to use config.toml in your current directory, but it can be set to a different file using the -C option before specifying which task to run, e.g.:
$ hlsmallrna -C alternate.toml sort ...
The layout of a TOML file has a number of sections each with a number of configuration keys. Each key has an expected type (e.g. string, integer, floating point number etc.) and the program will check they are as expected before running. If one is wrong, the program will crash and alert the user of the problem. A list of these keys and example values can be found in the README.md included with the program.
Manually setting parameters of internal commands
For advanced uses, you may want to add parameters to the internal commands run. This can be done through variables in the cli-tools section of the configuration file. A full list of these parameters can be found in the README.md provided with the code.
There is a main one that you can set, called <name>_params. It is a list of custom parameters to be added on to the command when called. For example if you wanted to set min-length to 10 for FastQC you could add the following to the configuration file:
[cli-tools]
[cli-tools.fastqc]
fastqc_params = ["--min_length", "10"]
Suppressing excess output
If you do not want to view the output from each of the commands, you can pass the -q option to the main script. (Before the task selection, like with -C). It is recommended to test the command first however, as this option sometimes suppresses useful error messages.
Note this doesn’t stop all output, only that of internal commands. If you don’t want any progress messages whatsoever, only errors and warnings from python, set -q twice. E.g. -q -q or -qq.
Troubleshooting
If it is taking too long to run or producing odd data, I strongly suggest deleting or moving the output directory somewhere else, then doing a clean run. This is because some steps use the data in the output directory as input, so may be combining old and new data.
When using the targetid command, samtools is used to filter out duplicates and create the FASTQ file outputs. Samtools is known to be picky about how the input file is layed out, so will fail. The program should keep running, but if you want these files, check the following in
your input files:
- For FASTA files, there is no @ at the start of any of the header strings
- There are no duplicate headers in any of the input files
CLI Reference
References of how to use the command line interface of the program and notes on some ways particular commands work that could trip you up.
General help message is as follows:
usage: hlsmallrna [-h] [-q] [-C PATH_TO_CONFIG]
{process,sort,extractnc,unitas,targetid,all} ...
Pipeline to process small RNAs
positional arguments:
{process,sort,extractnc,unitas,targetid,all}
process Preprocessing for the RNA
sort Find RNAs that align to a genome and sort them by
length
extractnc Extract the noncoding region from a fasta with a GFF
file
unitas Run unitas on split files and merge results
targetid Align small RNA to a number of genome features to find
out what is targeted
all Run process, sort and unitas one after the other
optional arguments:
-h, --help show this help message and exit
-q, --quiet Suppress output from intermediate commands
-C PATH_TO_CONFIG, --path-to-config PATH_TO_CONFIG
Path to the TOML format config file to use
Process
Docs:
process [-h] [-a ADAPTER] [-g FRONT] [-b ANYWHERE]
[-c CUTOFF]
small_rna
positional arguments:
small_rna Path to FASTQ containing the small RNA
optional arguments:
-h, --help show this help message and exit
-a ADAPTER, --adapter ADAPTER
Sequence of the adapter to remove from the 3' end
-g FRONT, --front FRONT
Sequence of the adapter to remove from the 5' end
-b ANYWHERE, --anywhere ANYWHERE
Sequence of the adapters to remove from both ends
-c CUTOFF, --cutoff CUTOFF
Quality cutoff to trim RNA sequences at
Input files:
Small_rna - fastq file containing raw RNA-seq data
Output files:
Fastqc/ - directory containing the raw output of FastQC
Trimmed_rna.fastq - file containing the raw RNA with adapters removed (only produced if an adapter sequence is provided)
Cut_sequences.fastq - sequences with low quality parts removed
Note that the parts of this step can easily be skipped. If you want to skip adapter trimming, don’t specify any adapters with -a, -g or -c. If you don’t want to run fastQC set -c 0. The program knows that in these cases, the steps do not have to be run.
Sort
sort [-h] [-d CDS] [-l MIN_LENGTH] [-x MAX_LENGTH]
[-m REF_MISMATCHES] [--disable-alignment]
small_rna [genome]
positional arguments:
small_rna Path to FASTQ containing the small RNA
genome Genome to align against
optional arguments:
-h, --help show this help message and exit
-d CDS, --cds CDS Optional CDS region, also align this to the CDS region
as well as the genome
-l MIN_LENGTH, --min-length MIN_LENGTH
Minimum length to bin
-x MAX_LENGTH, --max-length MAX_LENGTH
Maximum length to bin
-m REF_MISMATCHES, --ref-mismatches REF_MISMATCHES
Number of mismatches to use in bowtie2, None for
default behaviour
--disable-alignment Skip the alignment to the reference genome step
Input files:
Small_rna - fastq file containing RNA with adapters removed
Genome - Reference genome of the species you are using to align against in fasta format
Output files:
Bbmap_index/ - contains the index of the reference genome created by bbmap
Mapped_sequences.fastq - sequences successfully mapped to the reference by bbmap
binned_rna/ - directory containing one fastq file for each length of sequence contained in the bbmap output
Rna_length_report.csv - table showing a summary of the RNAs by length and first base
Baseplot.png - plot of the length and first base of the RNAs, using the data in rna_length_report.csv
Baseplot_data.csv - raw data used to make baseplot.png to allow for easy regraphing
ExtractNC
Docs:
extarctnc [-h] genome gff_file
positional arguments:
genome FASTA containing the genome to extract from
gff_file GFF file containing annotations of CDS and mRNA regions
optional arguments:
-h, --help show this help message and exit
Input files:
Genome - file containing the whole genome of the species of interest
Gff_file - GFF3 file containing annotations at least for the mRNA and coding region (labelled CDS) for the genome
Output files:
Noncoding.fasta - FASTA file containing only the transcribed noncoding regions of the DNA
Note the method used here will not work if there a coding region is labelled outside a single mRNA region. The program does a check before running to make sure this is true and throws an exception before trying to run if so, to prevent nonsense results from being produced.
Unitas
Docs:
unitas [-h] [-d CDS] [-u UNSPLICED_TRANSCRIPTOME]
[-r [REFSEQ [REFSEQ ...]]] [-s SPECIES]
path_to_rnas
positional arguments:
path_to_rnas Path to the folder with varying length RNAs in
optional arguments:
-h, --help show this help message and exit
-d CDS, --cds CDS Optional CDS region, passed to unitas
-u UNSPLICED_TRANSCRIPTOME, --unspliced-transcriptome UNSPLICED_TRANSCRIPTOME
Optional, unspliced transcriptome, passed to unitas
-r [REFSEQ [REFSEQ ...]], --refseq [REFSEQ [REFSEQ ...]]
References for use with unitas
-s SPECIES, --species SPECIES
Species to set in unitas arguments
Input files:
Path_to_rnas - path to the directory containing fastq files of each set of small RNAs you want to procesus in unitas
Note: to get unitas to work, you need to set either the species to a compatible species or provide one or more files of reference small RNAs. These can either be done with the -s and -r optional arguments or by adding them to your configuration file (see below).
CDS - Coding sequence, pipeline automatically labels with Gene and combines with the unspliced transcriptome, optional
Unspliced_transcriptome - Unspliced transcriptome, pipeline automatically labels with Gene and combines with the CDS, optional
Output files:
unitas/ - directory containing the results for all the unitas runs performed by the program
Unitas_summery.csv - CSV file containing the combined summaries of all of the unitas runs (plus some extra numbers). It is recommended to open it in a spreadsheet package (e.g. microsoft excel, libreoffice calc) as there is a lot of data and it may look a mess in a text editor.
UnitasGraph.png - graph showing the categories allocated by unitas against the length of the small RNA
Unitas_graph_data.csv - raw data used for creating unitasGraph.png to allow for easy regraphing
TargetID
Docs:
targetid [-h] [-m MIN_SEQ_LENGTH] [-t TARGET_FILES [TARGET_FILES ...]]
[--num-mismatches NUM_MISMATCHES]
small_rna
positional arguments:
small_rna Path to the FASTQ containing the small RNA to find targets of
optional arguments:
-h, --help show this help message and exit
-m MIN_SEQ_LENGTH, --min-seq-length MIN_SEQ_LENGTH
Minimum sequence length to properly align
-t TARGET_FILES [TARGET_FILES ...], --target-files TARGET_FILES [TARGET_FILES ...]
Files containing genome features that could be targeted
--num-mismatches NUM_MISMATCHES
Number of mismatches to allow in the alignment, defaults to 0
Input:
small_rna - fastq file containing small RNA to look for targets of
target_files - one or more fastq files containing a list of potential target sequences
Output:
Bowtie_indexes/ - indexes produced of the target files with bowtie2
target_alignments/ - SAM files containing the results of alignment attempts and FASTQ files containing all the small RNA successfully aligned against the targets. One set of files are produced for each target file provided
rna_target_list.csv - list of RNA target pairs
Note that the samtools view command used in this step is the only one that has defaults pre-set in the samtools_view_params config key, namly -h -F 256 -F 4. They will be removed if you set the value differently in the config file. To only add parameters, you will need to set these ones again before the extra ones you want to set. E.g. samtools_view_params = ["-h", "-F", "256", "-F", "4", ...]
Make sure all of your target files are labelled with what they are, by prepending the name and | (the bar or pipe character) to the ID in the fasta file. If you don’t do this, it will not be nicely sorted in the results correctly.
Label for Unitas
usage: label_for_unitas [-h] [-o OUTPUT] label file_path
Program to automatically prepend unitas labels to all the sequences in a FASTA
file
positional arguments:
label Label to prepend to the sequences in a fasta file
file_path Path to the fasta file to label
optional arguments:
-h, --help show this help message and exit
-o OUTPUT, --output OUTPUT
File to output to, defaults to labelled.fasta
Build Coord Files
usage: build_coord_files [-h] [-a] [-q] [-s] [-b] [-g] [-r] [-c]
[--feature FEATURE] [--genome-coords GENOME_COORDS]
[--band-width BAND_WIDTH] [-o OUTPUT]
input_file
Convert common bioinformatics file formats into coordinate file to plot
positional arguments:
input_file File to convert, attempts to autodetect type
optional arguments:
-h, --help show this help message and exit
-a, --fasta Treat input as a FASTA file
-q, --fastq Treat input as a FASTQ file
-s, --sam Treat input as a SAM file
-b, --bam Treat input as a BAM file
-g, --gff Treat input as a GFF file
-r, --rm-fa-out Treat input as a RepeatMasker .fa.out file
-c, --scaffold-aware Merge scaffolds into one chromosome
--feature FEATURE Select a gff feature to use
--genome-coords GENOME_COORDS
File containing the coordinates of the genome
--band-width BAND_WIDTH
Size of the bands showing change in scaffolds
-o OUTPUT, --output OUTPUT
File to output to
Overlap SS
usage: overlap_ss [-h] [-q] genome rna_file_1 rna_file_2
Looks for overlaps of two classes of RNA on the same strand of DNA
positional arguments:
genome FASTA file containing the genome of the organism
rna_file_1 File containing the RNA to use as a base
rna_file_2 File containing the RNA to look for overlaps with
optional arguments:
-h, --help show this help message and exit
-q, --quiet Suppress output from intermediate commands
\ No newline at end of file
diff --git a/docs/html/images/image1.png b/docs/html/images/image1.png
deleted file mode 100644
index c90ed2e..0000000
Binary files a/docs/html/images/image1.png and /dev/null differ
diff --git a/docs/html/images/image2.png b/docs/html/images/image2.png
deleted file mode 100644
index c9bc24d..0000000
Binary files a/docs/html/images/image2.png and /dev/null differ
diff --git a/docs/images/chromplot1.jpg b/docs/images/chromplot1.jpg
new file mode 100644
index 0000000..cc8c28c
Binary files /dev/null and b/docs/images/chromplot1.jpg differ
diff --git a/docs/images/chromplot2.jpg b/docs/images/chromplot2.jpg
new file mode 100644
index 0000000..af713f1
Binary files /dev/null and b/docs/images/chromplot2.jpg differ
diff --git a/docs/index.md b/docs/index.md
new file mode 100644
index 0000000..ad1bc63
--- /dev/null
+++ b/docs/index.md
@@ -0,0 +1,46 @@
+---
+hide:
+ - navigation
+ - toc
+---
+# Welcome to the Hunt Lab Small RNA Pipeline Documentation
+
+This documentation is for version 2 of the pipeline, select 1.1.0 from the version dropdown for the version 1 documentation.
+
+This is a pipeline developed in [Vicky Hunt's Lab](https://www.vickyhuntlab.org/) at the [University of Bath](https://www.bath.ac.uk), UK to analyse small RNA data in a quick and consistent way between projects. It is deliberately designed to be highly configurable so it can be easily applied to a wide range of scenarios, though most of our work is on nematode and nematomorph worms.
+
+# Background
+
+Most of the original methodology was devised by Mona Suleiman, Dominika Lastik and Vicky Hunt for the papers below. This was then developed into a usable pipeline by Kieran Reynolds in 2022 and updated to version 2 in 2025.
+
+> Suleiman, M., Kounosu, A., Murcott, B. et al. piRNA-like small RNAs target transposable elements in a Clade IV parasitic nematode. Sci Rep 12, 10156 (2022). [https://doi.org/10.1038/s41598-022-14247-1](https://doi.org/10.1038/s41598-022-14247-1)
+
+> Lastik, D., Kounosu, A., Dayi, M. et al. Small non-coding RNAs have predicted roles in reproductive biology and transposable element regulation in the parasitic worm Strongyloides venezuelensis. Sci Rep 15, 20608 (2025). [https://doi.org/10.1038/s41598-025-01968-2](https://doi.org/10.1038/s41598-025-01968-2)
+
+# Citation
+
+This pipeline is mainly a combination of existing tools so, doesn't have a paper attached to it but can be referenced with the repository URL and version number. For example:
+
+> This analysis was performed using HuntLab-smallRNA v2.0.0 ([https://github.com/Vicky-Hunt-Lab/HuntLab-smallRNA](https://github.com/Vicky-Hunt-Lab/HuntLab-smallRNA)).
+
+For the internal steps please cite the papers of the tools used:
+
+**Trim**:
+
+> Martin, Marcel. (2011). CUTADAPT removes adapter sequences from high-throughput sequencing reads. EMBnet.journal. 17. 10.14806/ej.17.1.200.
+
+**Sort**:
+
+> Langmead, B., Salzberg, S. Fast gapped-read alignment with Bowtie 2. Nat Methods 9, 357-359 (2012). https://doi.org/10.1038/nmeth.1923
+
+**Unitas**:
+
+> Gebert, D., Hewel, C. & Rosenkranz, D. unitas: the universal tool for annotation of small RNAs. BMC Genomics 18, 644 (2017). https://doi.org/10.1186/s12864-017-4031-9
+
+**Targetid**:
+
+> Langmead, B., Salzberg, S. Fast gapped-read alignment with Bowtie 2. Nat Methods 9, 357-359 (2012). https://doi.org/10.1038/nmeth.1923
+
+**Enrich**:
+
+> Carlos P Cantalapiedra, Ana Hernández-Plaza, Ivica Letunic, Peer Bork, Jaime Huerta-Cepas, eggNOG-mapper v2: Functional Annotation, Orthology Assignments, and Domain Prediction at the Metagenomic Scale, Molecular Biology and Evolution, Volume 38, Issue 12, December 2021, Pages 5825-5829, https://doi.org/10.1093/molbev/msab293
diff --git a/docs/installation.md b/docs/installation.md
new file mode 100644
index 0000000..05d5d5e
--- /dev/null
+++ b/docs/installation.md
@@ -0,0 +1,125 @@
+# Installation
+
+There are currently three different ways to install the pipeline, to allow for flexibility. The apptainer method is recommended as it is the quickest and most reliable, but the other two work if that is not possible.
+
+## Installing using Apptainer or Docker
+A containerised version of the pipeline is provided to allow for easy installation. Apptainer is recommended, as it allows for running with relative paths, but it will run in Docker as well. If you do not have these installed, see their websites for installation instructions:
+
+- Apptainer: [https://apptainer.org/docs/admin/main/installation.html](https://apptainer.org/docs/admin/main/installation.html)
+- Docker: [https://docs.docker.com/](https://docs.docker.com/)
+
+### Installing with Apptainer
+Create a SIF file with the code by running the following:
+
+```
+$ apptainer pull hlsmallrna.sif docker://kieranr51/huntlab-smallrna
+```
+
+To test it works try displaying the help message with the following command:
+
+```
+$ apptainer run hlsmallrna.sif hlsmallrna --help
+```
+
+If this displays the help message, you are good to go. All subsequent commands can be run with:
+
+```
+$ apptainer run -B $HOME:$HOME hlsmallrna.sif [COMMAND]
+```
+
+!!! note "Binding Directories in Apptainer"
+ To allow access to your files to the pipeline in apptainer you will need to bind the directories for you input and output. The above command binds your home directory, but you may need to change this if your files are outside of that directory. See [https://apptainer.org/docs/user/main/bind_paths_and_mounts.html](https://apptainer.org/docs/user/main/bind_paths_and_mounts.html) for more information.
+
+### Installing with Docker
+
+First pull the container with:
+
+```
+$ docker pull kieranr51/huntlab-smallrna
+```
+
+To test it works try displaying the help message with the following command:
+
+```
+$ docker run kieranr51/huntlab-smallrna hlsmallrna --help
+```
+
+If this displays the help message, you are good to go. All subsequent commands can be run with:
+
+```
+$ docker run -v $HOME:$HOME kieranr51/huntlab-smallrna [COMMAND]
+```
+
+!!! note "Paths and Mounting Directories in Docker"
+ To allow access to your files to the pipeline in apptainer you will need to mount the directories for you input and output. The above command mounts your home directory, but you may need to change this if your files are outside of that directory. See [https://docs.docker.com/engine/storage/bind-mounts/](https://docs.docker.com/engine/storage/bind-mounts/) for more information.
+
+ Docker hides the current path you are working in from the program as well, so you will have to specify the absolute path to your input and output files to get it to work. E.g. if your config file is in /home/user/smallrna/config.yml you would have to run:
+
+ ```
+ $ docker run -v $HOME:$HOME kieranr51/huntlab-smallrna hlsmallrna -o /home/user/smallrna/hlsmallrna_output /home/user/smallrna/config.yml
+ ```
+
+## Installing using Mamba or Conda
+
+If you can't use the containers, using mamba or conda to set up dependencies is the next supported route. Mamba is preferred. If you do not have these installed see the following links:
+
+- Mamba: [https://mamba.readthedocs.io/en/latest/installation/mamba-installation.html](https://mamba.readthedocs.io/en/latest/installation/mamba-installation.html)
+- Conda: [https://docs.conda.io/projects/conda/en/stable/user-guide/install/index.html](https://docs.conda.io/projects/conda/en/stable/user-guide/install/index.html)
+
+Then download the latest version of HuntLab-smallRNA from [https://github.com/Vicky-Hunt-Lab/HuntLab-smallRNA/releases](https://github.com/Vicky-Hunt-Lab/HuntLab-smallRNA/releases). Extract the zipped file, then navigate into the directory. To install run one of the following commands:
+
+```
+$ mamba env create -f environment.yml
+```
+
+```
+$ conda env create -f environment.yml
+```
+
+This will create an environment called huntlab-smallrna with the pipeline and dependencies installed into it. To test activate it with one of these commands:
+
+```
+$ mamba activate huntlab-smallrna
+```
+
+```
+$ conda activate huntlab-smallrna
+```
+
+Note you will have to run the above command each time you start your computer, before you use this software. Then you should be able to display the help message with the command:
+
+```
+$ hlsmallrna --help
+```
+
+## Manual Installation
+If none of the above are possible, you can retrieve each dependency and manually install it. See the following list of dependencies to install:
+
+- Python 3 - [https://www.python.org/](https://www.python.org/)
+- Pyyaml - [https://pyyaml.org/](https://pyyaml.org/)
+- Pysam - [https://pysam.readthedocs.io/en/latest/api.html](https://pysam.readthedocs.io/en/latest/api.html)
+- Biopython - [https://biopython.org/](https://biopython.org/)
+- Gffutils - [https://daler.github.io/gffutils/](https://daler.github.io/gffutils/)
+- Scipy - [https://scipy.org/](https://scipy.org/)
+- Matplotlib - [https://matplotlib.org/](https://matplotlib.org/)
+- Cutadapt - [https://cutadapt.readthedocs.io/en/stable/](https://cutadapt.readthedocs.io/en/stable/)
+- Unitas - [https://www.smallrnagroup.uni-mainz.de/software.html](https://www.smallrnagroup.uni-mainz.de/software.html)
+- Bowtie2 - [https://bowtie-bio.sourceforge.net/bowtie2/index.shtml](https://bowtie-bio.sourceforge.net/bowtie2/index.shtml)
+- Samtools - [http://www.htslib.org/](http://www.htslib.org/)
+- Bedtools - [https://bedtools.readthedocs.io/en/latest/](https://bedtools.readthedocs.io/en/latest/)
+- Eggnog-mapper - [https://github.com/eggnogdb/eggnog-mapper](https://github.com/eggnogdb/eggnog-mapper)
+- R -[https://www.r-project.org/](https://www.r-project.org/)
+- topGO - [https://bioconductor.org/packages/3.22/bioc/html/topGO.html](https://bioconductor.org/packages/3.22/bioc/html/topGO.html)
+- pip - [https://packaging.python.org/en/latest/tutorials/installing-packages/](https://packaging.python.org/en/latest/tutorials/installing-packages/)
+
+Once all of these are installed you need to download and extract the latest version of HuntLab-smallRNA from [https://github.com/Vicky-Hunt-Lab/HuntLab-smallRNA/releases](https://github.com/Vicky-Hunt-Lab/HuntLab-smallRNA/releases). Navigate into the extracted directory and run the following command to install:
+
+```
+$ pip install .
+```
+
+Then you should be able to display the help message with the command:
+
+```
+$ hlsmallrna --help
+```
\ No newline at end of file
diff --git a/docs/other-tasks.md b/docs/other-tasks.md
new file mode 100644
index 0000000..b349ada
--- /dev/null
+++ b/docs/other-tasks.md
@@ -0,0 +1,205 @@
+# Other Tasks
+
+This page covers some other bioinformatics tasks you may want to combine with the pipeline, using both programs included with and separate from the pipeline.
+
+## Plotting Genome Location
+
+### Basic Plot
+
+A common task you might want to do is plotting the location of things within the genome, for example, genome features you are interested in or the positions of small RNAs you have aligned to the genome. To make this easy, a script is included in this environment called `build_coord_files`. This automatically generates tsv files for use with the R package `chromPlot`. This section explains how you can generate a plot with it.
+
+Firstly, you need at least two files: a full genome in a FASTA or FASTQ file and a set of locations you are interested in, either as an alignment output (in SAM or BAM format) or as feature annotations in GFF format. To start produce the genome coordinate file with the command: (assuming `genome.fasta` is your genome file)
+
+```
+$ build_coord_files -o genome_coord.tsv genome.fasta
+```
+
+By default, the script will figure out what type of output to produce by the file extension. But if it can’t, the user will need to provide a flag to tell it what type of input you are using (in this case `--fasta`). Next we produce the coordinates file for the item of interest, assuming it is in `alignment.sam` you can run:
+
+```
+$ build_coord_files -o alignment_coord.tsv alignment.sam
+```
+
+This will produce two files, one for coordinates on the sense strand and one for coordinates on the antisense strand. Finally, we can use these files to plot the result. If it isn’t already, install chromPlot with bioconducter using the following R command:
+
+```r
+> BiocManager::install("chromPlot")
+```
+
+Then you should be able to create a plot using the following R code (of two example chromosomes called "SRAE_chr1" and "SRAE_chr2"):
+
+```r
+chrom_coords <- read.csv("/path/to/genome_coord.tsv", sep = "\t")
+sense_pos <- read.csv("/path/to/sense_alignment_coord.tsv", sep = "\t")
+antisense_pos <- read.csv("/path/to/antisense_alignment_coord.tsv", sep = "\t")
+
+library(chromPlot)
+
+chromPlot(
+ gaps = chrom_coords, annot1 = sense_pos, annot2 = antisense_pos, plotRndchr = T,
+ chr = c("SRAE_chr1", "SRAE_chr2"), bin = 50000,
+ chrSide = c(-1, 1, -1, -1, 1, -1, -1, 1)
+)
+```
+
+The script breaks down as follows:
+
+- First three lines load in the data we just produced with `build_coord_files`
+- Following line loads the chromPlot library into R
+- The final line plots the data with chromPlot:
+- `gaps` defines what each chromosome should look like
+- `annot1` and `annot2` define the annotation data to show
+- `plotRndchr` disables the assumption in chromPlot that the chromosomes are labelled with numbers
+- `chr` selects which chromosomes to plot by specifying an array of names
+- `bin` sets the bin size to group the features into when plotting
+- `chrSide` defines which side of the chromosome to plot each feature, the only change we have made from the default is setting the second item to `1` as opposed to `-1`. This ensures `annot2` is plotted on the opposite side to `annot1`
+
+The resulting plot will look something like this:
+
+
+
+### Dealing with Multiple Scaffolds
+
+In more complex cases, one or more of your chromosomes may be split into multiple scaffolds. `build_coord_files` has a scaffold aware mode that can automatically try to merge these, but it will only work if your scaffold’s IDs are in the form of `_scaffold` e.g. `SRAE_chrX_scaffold1`. Once you have this, you can generate the data files with the `-c` flag set, as follows:
+
+```
+$ build_coord_files -c -o genome_coord.tsv genome.fasta
+$ build_coord_files -c --genome-coords genome_coord.tsv -o alignment_coord.tsv alignment.sam
+```
+
+Note you need to provide the genome coords file when generating files for what you are plotting, this allows for a correction of coordinates to be made from where scaffolds are merged. In addition to the regular files, this produces `scaffold_info_genome_coord.tsv`, that can be used to add scaffold dividers to your plots. To do this, the following R script can be run:
+
+```r
+chrom_coords <- read.csv("/path/to/genome_coord.tsv", sep = "\t")
+sense_pos <- read.csv("/path/to/sense_alignment_coord.tsv", sep = "\t")
+antisense_pos <- read.csv("/path/to/antisense_alignment_coord.tsv", sep = "\t")
+scaffolds <- read.csv("/path/to/scaffold_info_genome_coord.tsv", sep = "\t")
+
+library(chromPlot)
+
+chromPlot(
+ gaps = chrom_coords, annot1 = sense_pos, annot2 = antisense_pos, plotRndchr = T,
+ chr = c("SRAE_chr1", "SRAE_chr2"), bin = 50000, bands = scaffolds,
+ chrSide = c(-1, 1, -1, -1, 1, -1, -1, 1)
+)
+```
+
+Compared to the last script, this adds:
+
+- The fourth line, that loads in the scaffold data
+- The bands argument, that draws red lines at the end of each scaffold
+
+Note that, by default the bands are 20,000 bp wide, as that looks good on the plot. At that size they may not appear if you make the bins too much bigger. Therefore, if you plan to make the bins bigger, pass a larger value to the `--band-width` flag when creating the genome file. That makes the bands bigger, which results in a more visible line on the plot.
+
+A plot produced by this method will look something like:
+
+
+ 
+## Signature Identification
+
+TODO
+
+## AT/GC Richness
+
+Nucleotide richness and the presence of any conserved motifs is identified using [WebLogo](https://weblogo.threeplusone.com). Sequences need to be the same length and can be either uploaded using a fasta file or can be pasted in.
+
+## Upstream motif (e.g. piRNA)
+
+1. Map to genome with bowtie2, convert sam to bam with `samtools`, convert bamToBed with `bedtools`
+
+```
+$ bowtie2 -x genome.fasta -f -U sequences.fasta -S sequences_mapped.sam
+$ samtools view sequences_mapped.sam -o sequences_mapped.bam
+$ bamToBed -i sequences.bam > sequences.bed
+```
+
+2. Extract flanking sequences with flankBed in `bedtools` (need to use bedtools format genome file - generate with faidx and cut):
+
+```
+$ samtools faidx genome.fasta
+$ cut -f 1,2 genome.fa.fai > chrom.sizes
+$ flankBed -i sequences.bed -g chrom.sizes -b 100 > sequences_100nt_flank.bed
+```
+
+or for 60 upstream only:
+
+```
+$ flankBed -i sequences.bed -g chrom.sizes -l 60 -r 0 > sequnences_60nt.bed
+```
+
+3. Extract fasta sequences from bed fastFromBed with bedtools.
+
+```
+$ fastaFromBed -fi genome.fasta -bed sequences_100nt_Flank.bed -fo sequences_100nt_Flank_seqs.fasta
+```
+
+## Ping Pong signature
+
+To find ping-pong signature, Unitas with the option -pp can be used to find 5’ overlaps of mapped sequence reads to the genome and calculate a Z-score for the enrichment of 10 bp overlaps. `map.file` should be in SAM or ELAND3 format.
+
+```
+$ perl unitas.pl -pp -input map.file –species x –refseq reference.fasta
+```
+
+## Dicer Signature (StepRNA)
+
+To install stepRNA, see the documentation on [the GitHub page](https://github.com/Vicky-Hunt-Lab/stepRNA).
+
+stepRNA is able to identify a Dicer processing signature from a small RNA sequence dataset. It requires two FASTA files as input that have had their adapters already trimmed from the ends.
+
+```
+$ stepRNA --reference REFERENCE.fa --reads READS.fa
+```
+
+- READS.fa should contain all of the small RNA sequencing reads.
+- REFERENCE.fa is usually filtered to investigate small RNAs of interest e.g. 26G sRNAs.
+
+Identical reads in the READ and REFERENCE files can also be removed before searching for a Dicer signature, if desired, by using `-e/--remove_exact`.
+
+Importantly the FASTA headers must be unique - this can be done by stepRNA with `-u/--make_unique`.
+
+stepRNA will then generate:
+
+- BAM alignment files (with different combinations of overhang lengths)
+- Overhang length CSV
+- Passenger Read Length CSV
+
+These can then be used to plot the distribution and/or further explore the reads that have aligned.
+
+Note: If the sRNA lengths in the read file are all the same i.e. 21 nt long, stepRNA will not be able to find overhangs due to the scoring system. In this case the `-m/--min_score` should be set to an appropriate value; we recommend sRNA length - 7nt (e.g. 21 - 6 = 15).
+
+Use `stepRNA --help` to bring up the command line help, a description of the methods can be found in the [stepRNA publication](https://www.frontiersin.org/journals/bioinformatics/articles/10.3389/fbinf.2022.994871/full).
+
+```
+stepRNA [-h] -r REFERENCE -q READS [-n NAME] [-d DIRECTORY]
+ [-m MIN_SCORE] [-e] [-u] [-j] [-V]
+
+Align a reference RNA file to read sequences. Output will be a set of CSV files containing information about the length of the reads, number of reads aligned to a reference sequence and the length of overhangs of the alignment.
+Reference RNA file will be automatically indexed
+
+Optional Arguments:
+ -h, --help show this help message and exit
+ -n NAME, --name NAME Prefix for the output files
+ -d DIRECTORY, --directory DIRECTORY
+ Directory to store the output files
+ -m MIN_SCORE, --min_score MIN_SCORE
+ Minimum score to accept, default is the shortest read
+ length
+
+Required Arguments:
+ -r REFERENCE, --reference REFERENCE
+ Path to the reference sequences
+ -q READS, --reads READS
+ Path to the read sequences
+
+Flags:
+ -e, --remove_exact Remove exact read matches to the reference sequence
+ -u, --make_unique Make FASTA headers unique in reference and reads i.e.
+ >Read_1 >Read_2
+ -j, --write_json Write count dictionaries to a JSON file
+ -V, --version Print version number then exit.
+```
+
+## miRNA
+
+TODO
\ No newline at end of file
diff --git a/docs/walkthrough/config.md b/docs/walkthrough/config.md
new file mode 100644
index 0000000..b1f0c0e
--- /dev/null
+++ b/docs/walkthrough/config.md
@@ -0,0 +1,142 @@
+# The Config File
+
+HuntLab-smallRNA 2 uses an overhauled config system that is different to HuntLab-smallRNA 1. This is based on a YAML file and the idea that anything that affects the analysis should be in the config file, while anything that doesn't is a command line argument. Alongside this it has been made easier to run multiple stages in one command. An annotated version of a minimal config file to run all 4 parts of the analysis is as follows:
+
+```yaml
+# The small RNA FASTQ file to use as input to the pipeline
+smallRNA_fastq: smallrna.fastq
+
+# The trim section, as it is present trim will be run
+trim:
+ # Use the bulit-in adapaters for the Qiagen kit
+ kit: qiagen
+
+# The sort section, as it is present sort will be run
+sort:
+ # The genome to align to filter small RNA
+ genome: genome.fasta
+
+# The unitas section, as it is present unitas will be run
+unitas:
+ # files to pass to unitas refseq and the labels to give them
+ refseq:
+ - miRNA: test/miRNA.fasta
+ - piRNA: test/piRNA.fasta
+ - tRNA: test/tRNA.fasta
+ - TE: test/transposable_elements.fasta
+
+# The targetid section, as it is present targetid will be run
+targetid:
+ # Files containing the interesting targets to align to
+ target_files:
+ - test/file1.fasta
+ - test/file2.fasta
+```
+
+!!! note "Absolute and Reletive Paths"
+
+ When parsing the config file, the small RNA pipeline assumes that all paths in the config file are reletive to the directory the config file is in. So if you move the config file, things will break. If you need to move the config file, use absolute paths (that start with `/`) instead.
+
+To run this, you could save the file as `hlsmallrna_config.yml` and run:
+
+```
+$ hlsmallrna hlsmallrna_config.yml
+```
+
+!!! tip
+ Remember if you install via apptainer or docker you will need to prepend the command from the install page. For apptainer this command becomes:
+ ```
+ $ apptainer run -B $HOME:$HOME hlsmallrna.sif hlsmallrna hlsmallrna_config.yml
+ ```
+
+The results would then appear in the directory `hlsmallrna_output/`, presuming all the sequence data referenced is real. A more comprehensive file to do a similar analysis to the first one is:
+
+```yaml
+# The small RNA FASTQ file to use as input to the pipeline
+smallRNA_fastq: smallrna.fastq
+# CDS FASTA of the species of interest, can be used in sort, unitas and targetid
+cds: cds.fasta
+# unspliced transcriptome FASTA of the species of interest, can be used in unitas and targetid
+unspliced_transcriptome: unspliced.fasta
+
+
+# The trim section, as it is present trim will be run
+trim:
+ # The 5’ adapter sequence to trim
+ 5_prime: ACGTTTAG
+ # The 3’ adapter sequence to trim
+ 3_prime: CGTAGGAT
+ # The quality filter cutoff to use
+ min_quality: 20
+
+# The sort section, as it is present sort will be run
+sort:
+ # The genome to align to filter small RNA
+ genome: genome.fasta
+ # If True, also align to the CDS file above
+ align_to_cds: True
+ # Minimum length of small RNA sequence to keep
+ min_length: 15
+ # Maximum length of small RNA sequence to keep
+ max_length: 50
+ # Max number to allow when aligning to the genome and CDS
+ mismatches: 0
+
+# The unitas section, as it is present unitas will be run
+unitas:
+ # files to pass to unitas refseq and the labels to give them
+ refseq:
+ # add genes to unitas, special keyword to add using cds and unspliced_transcriptome specified earlier
+ - gene
+ - miRNA: test/miRNA.fasta
+ - piRNA: test/piRNA.fasta
+ - tRNA: test/tRNA.fasta
+ - TE: test/transposable_elements.fasta
+ # species name to pass to unitas
+ species: x
+
+# The targetid section, as it is present targetid will be run
+targetid:
+ # Minimum length of small RNA, higher number will speed up runtime, used to calculate bowtie2 seed length
+ min_seq_length: 5
+ # Number of mismatches to allow when aligning to the targets
+ mismatches: 0
+ # Files containing the interesting targets to align to
+ target_files:
+ - test/file1.fasta
+ - test/file2.fasta
+ # If present also enrich GO terms, KEGG pathways and Pfams of targets
+ enrich:
+ # Path to the eggnog-mapper data dir
+ eggnog_data_dir: /home/user/eggnog-mapper-data
+ # Target files to ignore during enrichment
+ exclude_files:
+ - test/file2.fasta
+```
+
+If you want to skip a step, just leave it out of the file. e.g. if you only want to run sort and unitas, you could use the file:
+
+```yaml
+smallRNA_fastq: smallrna.fastq
+sort:
+ genome: genome.fasta
+unitas:
+ refseq:
+ - miRNA: test/miRNA.fasta
+ - piRNA: test/piRNA.fasta
+ - tRNA: test/tRNA.fasta
+ - TE: test/transposable_elements.fasta
+```
+
+!!! tip "Command Line Arguments"
+ While most things are now specified in the config file, there are still a few CLI arguments that are as follows:
+
+ `-o`, `--ouput` - directory to write pipeline output to (default: hlsmallrna_output)
+
+ `-t`, `--threads` - number of threads to use when running analysis tools (default: 4)
+
+ `-k`, `--keep-files` - If set, don't delete intermediate files, useful for debugging or if you need an intermediate file later
+
+ `-v`, `--verbose` - makes the pipeline print out the output of the intermediate commands
+
+The following sections will explain each parameter in the config file in detail.
diff --git a/docs/walkthrough/index.md b/docs/walkthrough/index.md
new file mode 100644
index 0000000..2fdc160
--- /dev/null
+++ b/docs/walkthrough/index.md
@@ -0,0 +1,26 @@
+# Walkthrough
+
+This walkthrough covers all of the major tasks of the pipeline and how they can be performed, either together or separately and considerations for your analysis. To start, here is an overview of the analysis the pipeline can perform:
+
+## Trim
+
+This filters low quality ends, trims 5' adapters and then trims 3' adapters. Any reads without a detectable 3' adapter are removed, as they may not cover a full small RNA.
+
+## Sort
+
+Firstly, aligns the small RNA to the genome and removes any that don’t align to remove contamination. Then sorts into FASTQ files of each length within a specified range. Finally, plots a graph of frequency of each length and first base of small RNA to give an overview of your sample.
+
+## Unitas
+
+This step takes the files of each small RNA length from sort and runs unitas on them to classify the origins of the small RNA. Then plots a graph of the origins of each length of small RNA.
+
+## TargetID
+
+This aligns the reverse complement of the small RNA exactly to FASTAs of interesting targets (e.g. genes in the organism). These confident targets can then optionally be enriched for GO terms, KEGG pathways and Pfams using eggnog-mapper to allocate.
+
+!!! note "Pipeline Stage Ordering"
+ While you do not have to run all of these parts of the pipeline, they will always run in the order:
+
+ Trim -> Sort -> Unitas -> TargetID -> Enrich
+
+ If you want to do a different order you will need to write multiple config files and run the pipeline multiple times.
diff --git a/docs/walkthrough/prereq.md b/docs/walkthrough/prereq.md
new file mode 100644
index 0000000..9432cb7
--- /dev/null
+++ b/docs/walkthrough/prereq.md
@@ -0,0 +1,63 @@
+# Prerequisites and Considerations
+
+Before we get on to how to run the pipeline this page covers what files are needed to run the pipeline and things to think about to inform the parameters to set. Read each section you want to run carefully.
+
+## Trim
+The trimming in the pipeline is what we think is best practice for the small RNA data we use - remove low quality ends and trim any 5' adapters first, then trim the 3' adapter and remove any sequences that lack it - as they are likely incomplete or noise. If you want to do something different, you can do this and provide the already trimmed small RNAs to the pipeline to run the other analysis. To run this, you will need the following data:
+
+- Small RNA-seq reads in FASTQ or FASTA format
+- 3' adapter sequence from small RNA library prep
+- Optionally, a 5' adapter sequence and quality cutoff to use
+
+!!! tip "Small RNA Library Prep Kits"
+ For the Quiagen and NEB kits the best adapters to use are built into the pipeline. You just need to know the kit and if we have data in the forward or reverse direction.
+
+!!! note "Paired End Sequencing"
+ Most small RNA-seq provides a single sequence, but if you have paired end data, you should only provide the read that contains the sequence of the small RNA to the program (Usually the forward read), not it's reverse complement.
+
+## Sort
+
+This takes trimmed reads and aligns them to the genome and optionally a set of CDS sequences to catch small RNA derived from spliced genes. Then sorts them into length and plots length and first base to give you an overview of the small RNAs you have. It also outputs a summary table on the number and percentage of small RNA reads aligned and FASTQs of small RNAs of each class (length and first base) for easier downstream analysis. For this you'll need:
+
+- A FASTQ or FASTA of trimmed small RNA sequences (can be automatically retrieved if trim is run)
+- A genome assembly to align the reads to
+- Optionally, a set of CDS sequences to also align reads to
+- Optionally, A minimum and maximum length of small RNA to consider, otherwise defaults to 18 - 30 nt
+- Optionally, A number of mismatches to use when aligning, otherwise defaults to 0
+
+!!! tip "Mismatches"
+ By default, sort allows no mismatches or gaps, which is fine if you have a good reference genome that is of a close relative to the sample small RNA-seq was performed on and want to be conservative about what you retain. If you know your genome is likely to have lots of SNPs from the reference or are uncertain about how good you assembly is, you may want to up the amount of mismatches allowed to compensate for this.
+
+!!! warning "Downstream use of BLAST"
+ While it seems like a good to run small RNAs through BLAST downstream, if you do this you need to be careful of the parameters used, particularly word size and E value. This is because BLAST will not return any hits if your sequence is shorter than the word size set (often above 20nt) and the short size of small RNA may mean even close to exact marches have a high E value. Due to this a basic search may return no hits, even if huntlab-smallrna has found an exact match in your genome.
+
+## Unitas
+
+Unitas runs unitas on a number of files representing small RNA lengths and plots a graphs of the origin or class for each small RNA length. This requires the following files:
+
+- A directory of small RNAs of different lengths named in the form length.fastq - where X is the length of the small RNA. This can be automatically derived if sort has been run
+- One or both of:
+ - A species name in the form supported by unitas - see examples in [https://www.smallrnagroup.uni-mainz.de/software/unitas_documentation_1.7.0.pdf](https://www.smallrnagroup.uni-mainz.de/software/unitas_documentation_1.7.0.pdf)
+ - One or more files containing known sequences of a particular class (e.g. miRNA, rRNA, tRNA, Transposable Elements etc.)
+- Optionally a CDS and/or unspliced transcriptome files to supplement the above
+
+!!! tip "Labelling for Unitas"
+ Unitas expects files that have the class of each sequence labelled in their ID e.g. for TEs with ids te1 and te2, the file would look like:
+ ```
+ >TE|te1
+ ACGT...
+ >TE|te2
+ CGAT...
+ ```
+ Unless you have mixed files, the pipeline can label this for you, so don't worry about adding that yet. See the unitas section for information on this.
+
+## TargetID
+This stage aligns the reverse complement exactly to lists of potential targets, to find likely targets of a small RNA, outputting both and overall list and by class (e.g. 22Us). Then it can optionally enrich the targets GO terms, KEGG pathways and Pfams to help understand the tragets. For this you will need:
+
+- FASTQ file of small RNAs or a directory or sorted small RNAs - can be derived from trim or sort if run
+- One or more FASTA files containing sequences of potential targets (e.g. genes, TEs, host genes etc.)
+- Optionally, the number of mismatches to tolerate when aligning, if not 0
+- If you want to run the enrichment, an eggnog-mapper data directiory (see targetid section on how to create this)
+- Optionally, a list of target files to exclude from the enrichment
+
+Now you know what you need you can continue to the next page of the tutorial to understand the config file you need to write.
\ No newline at end of file
diff --git a/docs/walkthrough/sort.md b/docs/walkthrough/sort.md
new file mode 100644
index 0000000..a0b55b6
--- /dev/null
+++ b/docs/walkthrough/sort.md
@@ -0,0 +1,121 @@
+# Running Sort
+
+## Basic Configuration
+
+At minimum, sort requires a set of trimmed reads and a genome to align to from the species the small RNAs came from. If trim is in the same config file, the set to trimmed reads will be taken from the output of trim. Examples are as follows.
+
+Without running trim.
+```yaml
+smallRNA_fastq: smallrna_trimmed.fastq
+sort:
+ genome: genome.fasta
+```
+
+With running trim.
+```yaml
+smallRNA_fastq: smallrna_untrimmed.fastq
+trim:
+ ...
+sort:
+ genome: genome.fasta
+```
+
+Once this config is written, it can be run as usual using:
+
+```
+$ hlsmallrna config.yml
+```
+
+!!! tip
+ Remember if you install via apptainer or docker you will need to prepend the command from the install page. For apptainer this command becomes:
+ ```
+ $ apptainer run -B $HOME:$HOME hlsmallrna.sif hlsmallrna config.yml
+ ```
+
+## Optional Parameters
+
+### align_to_cds
+
+If true aligns reads that don't align to the genome to the specified CDS file, keeping any that align. This allows for small RNA that are fragments of genes and align accross introns to be kept as they may not align to the genome.
+
+```yaml
+cds: cds_sequences.fasta
+...
+sort:
+ ...
+ align_to_cds: True
+```
+
+### min_length
+
+Sets a minimum sequence length to look at (default is 18 nt). Anything shorter is filtered during sorting. e.g. to filter anything below 10 nt set the following:
+
+```yaml
+sort:
+ min_length: 10
+```
+
+### max_length
+
+Sets a maximum sequence length to look at (default is 30 nt). Anything longer is filtered during sorting. e.g. to filter anything above 50 nt set the following:
+
+```yaml
+sort:
+ max_length: 50
+```
+
+### mismatches
+
+Number of mismatches to allow in bowtie2 when aligning to the genome and CDS, defaults to 0. See the [Considerations](prereq.md) page for information about why you'd want to change this. The example below sets this to allow 1 mismatch:
+
+```yaml
+sort:
+ mismatches: 1
+```
+
+!!! warning
+ The small RNA pipeline supresses gaps by makeing them equal to 100 mismatches. If you set this over 100, some gaps may be included in the alignment.
+
+## Output Files
+
+### alignment_report.tsv
+
+This is a table showing an overview of how many small RNA reads successfully aligned to the genome and optionally the CDS. It contains the following metrics:
+
+- Total Reads
+- Total Mapped
+- Percentage Mapped
+- Total Unmapped
+- Percentage Unmapped
+- Mapped to Genome
+- Percentage Mapped to Genome
+- Unmapped to Genome
+- Percentage Unmapped to Genome
+- Mapped to CDS
+- Percentage Mapped to CDS (but not to genome)
+- Unmapped to CDS or genome
+- Percentage Unmapped to CDS or genome
+
+### counts.tab
+
+Read counts is TSV format for each small RNA sequence, designed as input to a differential expression analysis, if that wanted.
+
+### rna_length_report.csv
+
+Table of counts of small RNA for each length and first base, in a easily human readable format.
+
+### baseplot.png
+
+Plot of frequency of length and first base in the dataset.
+
+### baseplot_data.csv
+
+That length and first base data used in baseplot.png - normalised to percentages. Should be used if the user wants to replot that graph or feed into downstream analysis.
+
+### binned_length_rna
+
+Directory containing FASTQ files for each RNA length - named lengthX.fastq where X is the length of small RNA in the file.
+
+### binned_group_rna
+
+Directory containing FASTQ files for each RNA group - defined by length and first base - named {length}{first_base}.fastq (e.g. 22G.fastq) where {length} is the length of small RNA in the file and {first_base} is the letter corresponding to the first base.
\ No newline at end of file
diff --git a/docs/walkthrough/targetid.md b/docs/walkthrough/targetid.md
new file mode 100644
index 0000000..f120f82
--- /dev/null
+++ b/docs/walkthrough/targetid.md
@@ -0,0 +1,119 @@
+# Running TargetID
+
+## Basic Configuration
+
+To run TargetID, you need two things:
+
+1. A small RNA file - if sort is specified the output of that is used - if not `smallRNA_fastq` is used.
+
+2. One or more files containing potential targets of small RNAs (e.g. genes, transposons, lncRNA etc.) specified in the `targetid` section as `target_files`.
+
+For example:
+
+```yaml
+smallRNA_fastq: small_rna.fastq
+targetid:
+ target_files:
+ - targets/lncRNA.fasta
+ - targets/genes.fasta
+```
+
+Once this config is written, it can be run as usual using:
+
+```
+$ hlsmallrna config.yml
+```
+
+!!! tip
+ Remember if you install via apptainer or docker you will need to prepend the command from the install page. For apptainer this command becomes:
+ ```
+ $ apptainer run -B $HOME:$HOME hlsmallrna.sif hlsmallrna config.yml
+ ```
+
+## Optional Parameters
+
+### min_seq_length
+
+This is used to set the seed length in bowtie2 when aligning to targets (defaults to 5). Larger numbers makes it fasta but less accurate for sequences close to this length. For example, to set to 18 add the following:
+
+```yaml
+targetid:
+ min_seq_length: 18
+```
+
+### mismatches
+
+Number of mismatches to allow in bowtie2 when aligning to the target files, defaults to 0. See the [Considerations](prereq.md) page for information about why you'd want to change this. The example below sets this to allow 1 mismatch:
+
+```yaml
+targetid:
+ mismatches: 1
+```
+
+### enrich - eggnog_data_dir
+
+The pipeline can optionally enrich GO terms, KEGG pathways and PFams for the targets of the small RNA. This is a rough analysis based on `emapper.py --itype CDS` mode that uses `diamond blastx` to map, so is likely only a rough analysis. May be useful if you want to have a quick look at if there are any common annotations within the targets though should be redone probably later. To run this you need to add the eggnog_data_dir to the config file. This directory can be created with the command:
+
+```
+$ mkdir -p /home/user/path/to/eggnog-data
+$ download_eggnog_data.py -y --data_dir /home/user/path/to/eggnog-data
+```
+
+!!! tip
+ Remember if you install via apptainer or docker you will need to prepend the command from the install page. For apptainer this command becomes:
+ ```
+ $ mkdir -p /home/user/path/to/eggnog-data
+ $ apptainer run -B $HOME:$HOME hlsmallrna.sif download_eggnog_data.py -y --data_dir /home/user/path/to/eggnog-data
+ ```
+Then add the following to the config file
+
+```yaml
+targetid:
+ ...
+ enrich:
+ eggnog_data_dir: /home/user/path/to/eggnog-data
+```
+
+### enrich - exclude
+
+If you have run targetid with targets it doesn't make sense to map to the eggnog database (e.g. lncRNA), you can exclude them with the exclude key, but must ensure the path exactly matches the one in the target_files list. e.g.
+
+```yaml
+targetid:
+ target_files:
+ - targets/genes.fasta
+ - targets/lncRNA.fasta
+ enrich:
+ eggnog_data_dir: /home/user/path/to/eggnog-data
+ exclude:
+ - targets/lncRNA.fasta
+```
+
+## Output Files
+### rna_target_list.tsv
+
+List of small RNA and there targets identified from the files specified.
+
+### targets_by_group
+
+Same information as `rna_target_list.tsv`, but split by small RNA group (length and first base).
+
+## Output Files - Enrich
+
+These files are only produced if the `enrich` key is present.
+
+### eggnog_mapper
+
+Results of the eggnog-mapper run on the target files.
+
+### targets_enriched_go_terms.tsv
+
+Enriched GO terms in the targets of small RNA vs the rest of the target files.
+
+### targets_enriched_kegg_pathways.tsv
+
+Enriched KEGG pathways in the targets of small RNA vs the rest of the target files.
+
+### targets_enriched_pfams.tsv
+
+Enriched PFams in the targets of small RNA vs the rest of the target files.
\ No newline at end of file
diff --git a/docs/walkthrough/trim.md b/docs/walkthrough/trim.md
new file mode 100644
index 0000000..2483b88
--- /dev/null
+++ b/docs/walkthrough/trim.md
@@ -0,0 +1,64 @@
+# Running Trim
+
+## Basic Configuration
+
+Trim requires a file of untrimmed small RNAs and the adpaters to trim - either explicitly or through the kit used - to be specified. Examples of config files for possible usage is as follows.
+
+Trimming adapters with the 5' sequence ACGTTTAGAT and the 3' sequence GCTAGGATA:
+
+```yaml
+smallRNA_fastq: smallrna_untrimmed.fastq
+trim:
+ 5_prime: ACGTTTAGAT
+ 3_prime: GCTAGGATA
+```
+
+Trimming adapters from the qiagen kit in forward direction:
+
+```yaml
+smallRNA_fastq: smallrna_untrimmed.fastq
+trim:
+ kit: qiagen
+```
+
+Trimming adapters from the qiagen kit in reverse direction:
+
+```yaml
+smallRNA_fastq: smallrna_untrimmed.fastq
+trim:
+ kit: qiagen-rev
+```
+
+Currently the valid values for kit are:
+
+- `qiagen` - qiagen kit, forward direction
+- `qiagen-rev` - qiagen kit, reverse direction
+- `neb` - NEB kit, forward direction
+- `neb-rev` - NEB kit, reverse direction
+
+Once this config is written, it can be run as usual using:
+
+```
+$ hlsmallrna config.yml
+```
+
+!!! tip
+ Remember if you install via apptainer or docker you will need to prepend the command from the install page. For apptainer this command becomes:
+ ```
+ $ apptainer run -B $HOME:$HOME hlsmallrna.sif hlsmallrna config.yml
+ ```
+
+## Optional Parameters
+### min_quality
+This sets the quality cutoff in cutadapt when filtering low quality ends. Defaults to 20 - example below sets it to 15:
+
+```yaml
+trim:
+ ...
+ min_quality: 15
+```
+
+## Output Files
+### trimmed_reads.fq
+
+This is a FASTQ file containing all the trimmed reads that passed the filters of quaility and had a 3' adapter present.
\ No newline at end of file
diff --git a/docs/walkthrough/unitas.md b/docs/walkthrough/unitas.md
new file mode 100644
index 0000000..58d8708
--- /dev/null
+++ b/docs/walkthrough/unitas.md
@@ -0,0 +1,148 @@
+# Running Unitas
+
+## Basic Configuration
+
+Unitas requires two parameters to work - a directory containing length sorted small RNA FASTQ and a set of refseq FASTA files to use a reference for the different types of small RNA or features that they may come from.
+
+The former (directory containing length sorted small RNA FASTQ) can be specified in two different ways -
+1. Directly though the `size_sorted_fastqs` parameter e.g.
+
+```yaml
+size_sorted_fastqs: sorted_small_rna
+unitas:
+ ...
+```
+
+2. Through having both `smallRNA_fastq` and `sort` specified, in which case it uses the output of `sort`. e.g.
+
+```yaml
+smallRNA_fastq: trimmed_smallrna.fastq
+sort:
+ ...
+unitas:
+ ...
+```
+
+The latter (set of refseq FASTA files) are specified as a list and can be written in three different forms:
+
+1. If the FASTA file is not already labelled for unitas in the IDs (i.e. ID values do not start with `[class_name]|`) the class name for the pipeline to add can be specified as a key and the file as a value. For example to add a file of miRNA called `microRNA.fasta` that should be labelled `miRNA` specify the following:
+
+```yaml
+unitas:
+ refseq:
+ - miRNA: microRNA.fasta
+```
+
+2. If the file already has unitas labels in the IDs (i.e. ID values start with `[class_name]|`) you only need to specify the path to the FASTA file. e.g. if you have a labeled set of tRNA in a file called `transferRNA.fasta` specify the following:
+
+```yaml
+unitas:
+ refseq:
+ - transferRNA.fasta
+```
+
+3. If you want to add gene features from `cds` and/or `unspliced_transcriptome` specify the special string `gene`. e.g.
+
+```yaml
+cds: cds.fasta
+unspliced_transcriptome: unspliced.fasta
+unitas:
+ refseq:
+ - gene
+```
+
+These can be mixed and matched, so a full example could be:
+
+```yaml
+smallRNA_fastq: trimmed_smallrna.fastq
+sort:
+ ...
+unitas:
+ refseq:
+ - gene
+ - transferRNA.fasta
+ - miRNA: microRNA.fasta
+```
+
+Once this config is written, it can be run as usual using:
+
+```
+$ hlsmallrna config.yml
+```
+
+!!! tip
+ Remember if you install via apptainer or docker you will need to prepend the command from the install page. For apptainer this command becomes:
+ ```
+ $ apptainer run -B $HOME:$HOME hlsmallrna.sif hlsmallrna config.yml
+ ```
+
+## Optional Parameters
+### species
+
+If unitas allows, you can set a species to have unitas download existing data on it's small RNAs from the internet to help with classification. For example:
+
+```yaml
+unitas:
+ species: caenorabditis_elegans
+```
+
+## Utility Programs
+
+Two utility programs to help with creating input for unitas are contained in the small RNA pipeline.
+
+### label_for_unitas
+
+This is a simple program that adds a class label that unitas can read to a FASTA file of sequences. It is recommened to not use this directly most of the time and instead use the key-value labelling syntax shown above. However, if you need a unitas labelled file, you can run the following:
+
+```
+$ label_for_unitas [class] [FASTA_file]
+```
+
+e.g. to label `transposons.fasta` with TE and write the result to `transposons_labelled.fasta`, you could run:
+
+```
+$ label_for_unitas TE transposons.fasta -o transposons_labelled.fasta
+```
+
+!!! tip
+ Remember if you install via apptainer or docker you will need to prepend the command from the install page. For apptainer this command becomes:
+ ```
+ $ apptainer run -B $HOME:$HOME hlsmallrna.sif label_for_unitas TE transposons.fasta -o transposons_labelled.fasta
+ ```
+
+
+### extract_nc
+
+If you want to test what small RNA come form noncoding reigons - outside of genes - specifically you can use this tool to extract those from the genome. It also adds the unitas label noncoding to the ID of the sequences for you.
+
+!!! warning
+ Most of the time you do not want to do this on the first run - as Unitas will likely double classify anything that aligns to something like a transposon. It is advised to take anything unitas has not classified on the first run and only run that on the noncoding if you want the information.
+
+It can be run as follows:
+
+```
+$ extract_nc [genome_FASTA] [annotation_GFF] -o [output_FASTA]
+```
+
+!!! tip
+ Remember if you install via apptainer or docker you will need to prepend the command from the install page. For apptainer this command becomes:
+ ```
+ $ apptainer run -B $HOME:$HOME hlsmallrna.sif extract_nc [genome_FASTA] [annotation_GFF] -o [output_FASTA]
+ ```
+
+## Output Files
+### unitas
+
+Directory containing the raw output of unitas for each length file run. Can be inspected to extract details for particular small RNAs.
+
+### unitas_summary.csv
+
+Human-readable unitas classification summary by length. Conatins the classifications for each RNA length on a row, allowing for inspection of all the details at once.
+
+### unitasGraph.png
+
+Graph of unitas classification by length - uses only top level unitas classes.
+
+### unitas_graph_data.csv
+
+Data for `unitasGraph.png`, in case it needs to be replotted, misses the more granular classes from `unitas_summary.csv`.
\ No newline at end of file
diff --git a/environment.yml b/environment.yml
index 427d8da..07298eb 100644
--- a/environment.yml
+++ b/environment.yml
@@ -3,22 +3,23 @@ channels:
- bioconda
- conda-forge
dependencies:
- - python=3.8
- - toml=0.10.*
- - pysam=0.19.*
- - biopython=1.*
- - gffutils=0.11.*
- - matplotlib=3.5.*
- - cutadapt=4.*
- - fastqc=0.11.*
- - bbmap=38.*
- - unitas=1.6.1
- - perl-archive-extract=0.88
- - bioconda/label/main::perl-lwp-simple=6.39
- - bowtie2=2.4.*
- - samtools=1.*
- - bedtools=2.*
+ - python
+ - pyyaml
+ - pysam
+ - biopython
+ - gffutils
+ - matplotlib
+ - cutadapt
+ - unitas
+ - perl-archive-extract
+ - bioconda/label/main::perl-lwp-simple
+ - bowtie2
+ - samtools
+ - bedtools
+ - eggnog-mapper
+ - scipy
+ - r
+ - bioconductor-topgo
- pip
- pip:
- .
-
diff --git a/example_config.yml b/example_config.yml
new file mode 100644
index 0000000..e3a2f40
--- /dev/null
+++ b/example_config.yml
@@ -0,0 +1,37 @@
+smallRNA_fastq: smallrna.fastq
+# size_sorted_fastqs: binned_reads
+cds: cds.fasta
+unspliced_transcriptome: unspliced.fasta
+trim:
+ # input: smallrna.fastq
+ kit: qiagen
+ # 5_prime: ACGTTTAG
+ # 3_prime: CGTAGGAT
+ min_quality: 20
+# new behaviour: output by first base as well e.g. file of 26G
+sort:
+ # input: output/trimmed_reads.fq
+ genome: genome.fasta
+ align_to_cds: True
+ min_length: 5
+ max_length: 50
+ mismatches: 0
+unitas:
+ # input: output/binned_rna
+ refseq:
+ - gene
+ - miRNA: test/miRNA.fasta
+ - piRNA: test/piRNA.fasta
+ - tRNA: test/tRNA.fasta
+ - TE: test/transposable_elements.fasta
+ species: x
+targetid:
+ min_seq_length: 5
+ target_files:
+ - test/file1.fasta
+ - test/file2.fasta
+ mismatches: 0
+ enrich:
+ eggnog_data_dir: /home/user/eggnog-mapper-data
+ exclude_files:
+ - test/file2.fasta
\ No newline at end of file
diff --git a/hlsmallrna/__init__.py b/hlsmallrna/__init__.py
index 1eed4a0..4e8ea52 100644
--- a/hlsmallrna/__init__.py
+++ b/hlsmallrna/__init__.py
@@ -1,4 +1,4 @@
-# Copyright 2022 Vicky Hunt Lab Members
+# Copyright 2022 - 2025 Vicky Hunt Lab Members
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
@@ -11,10 +11,11 @@
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
-from .__main__ import main, main_ssoverlap
+from .__main__ import main
+from .label_for_unitas import label_for_unitas_cli
def climain():
main()
-def ssoverlap_main():
- main_ssoverlap()
\ No newline at end of file
+def labelforunitas_main():
+ label_for_unitas_cli()
\ No newline at end of file
diff --git a/hlsmallrna/__main__.py b/hlsmallrna/__main__.py
index 28d5ed5..46dd3ce 100644
--- a/hlsmallrna/__main__.py
+++ b/hlsmallrna/__main__.py
@@ -1,5 +1,4 @@
-#!/usr/bin/env python3
-# Copyright 2022 Vicky Hunt Lab Members
+# Copyright 2022 - 2025 Vicky Hunt Lab Members
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
@@ -13,438 +12,156 @@
# See the License for the specific language governing permissions and
# limitations under the License.
-from ast import Num
-import os.path
+import os
import glob
import shutil
-from math import inf
+from functools import partial
from multiprocessing import Pool
from argparse import ArgumentParser
-from Bio import SeqIO, SeqRecord
-from gffutils.iterators import DataIterator
-
+from .validate_config import load_and_validate_input
from .trim import run_trim
-from .fastqc import run_fastqc, cut_rna_below_cutoff
from .genome_align import align_to_genome, bin_rna_size, graph_length
-from .create_noncoding import extract_noncoding
-from .unitas import copy_and_label_file, run_unitas_annotation, merge_summary, graph_unitas_classification_type
-from .targetid import revcomp_input_file, find_targets, build_summary_files
-from .ss_overlap import samestrand_overlap
-
-from .config import get_config_key, mkdir_if_not_exists, load_config, do_log
-
-def validate_file(path_to_file, file_type):
- '''
- Examine a file and determine if it is of the requested file type
- '''
- if file_type == 'directory':
- return os.path.isdir(path_to_file)
- elif file_type == 'fasta':
- try:
- item = next(SeqIO.parse(path_to_file, 'fasta'))
- return type(item) == SeqRecord.SeqRecord
- except:
- return False
- elif file_type == 'fastq':
- try:
- item = next(SeqIO.parse(path_to_file, 'fastq'))
- return type(item) == SeqRecord.SeqRecord
- except:
- return False
- elif file_type == 'gff':
- iter = DataIterator(path_to_file)
- return len(iter.peek) > 0
-
+from .label_for_unitas import label_file_for_unitas
+from .unitas import find_min_max_lengths, run_unitas_annotation, merge_summary, graph_unitas_classification_type
+from .targetid import merge_and_revcomp_input_files, find_targets, build_summary_files, do_term_enrichment
-def process_command(small_rna, adapter, front, anywhere, cutoff, quiet):
+def run_trim_reads(trim_config, program_paths, output, threads, verbose, keep_intermediate):
'''
- Code to run when the user chooses the process command
+ Call run_trim with the correct values based on the config
'''
+ run_trim(
+ trim_config['input'], output, five_prime_adapter=trim_config['5_prime'], three_prime_adapter=trim_config['3_prime'], quality_cutoff=trim_config['min_quality'],
+ threads=threads, cutadapt_path=program_paths['cutadapt'], keep_intermediate=keep_intermediate, verbose=verbose
+ )
- if not validate_file(small_rna, 'fastq'):
- print(f'Error: expected a small RNA FASTQ with at least one sequence, got {small_rna}')
- return False
-
- do_log(quiet, '==> Starting command Process')
- small_rna_path = small_rna
-
- if adapter is not None or front is not None or anywhere is not None:
- run_trim(small_rna_path, adapter, front, anywhere)
- small_rna_path = os.path.join(get_config_key('general', 'output_directory'), 'trimmed_rna.fastq')
- else:
- do_log(quiet, '==> No adapter sequence provided, skipping trim step')
-
- if cutoff > 0:
- run_fastqc(small_rna_path, quiet=quiet)
- cut_rna_below_cutoff(small_rna_path, cutoff, quiet=quiet)
- else:
- # create the cut_sequences.fastq, even if fastqc wasn't run
- shutil.copy2(small_rna_path, os.path.join(get_config_key('general', 'output_directory'), 'cut_sequences.fastq'))
-
- do_log(quiet, '==> Completed command Process')
-
-def sort_command(genome, small_rna, cds, min_length, max_length, num_mismatches, disable_align, threads, quiet):
+def run_sort_reads(sort_config, program_paths, output, threads, verbose, keep_intermediate):
'''
- Code to run when the user chooses the sort command
+ Call align_to_genome with the correct values based on the config
'''
-
- if not disable_align and not validate_file(genome, 'fasta'):
- print(f'Error: expected a genome in FASTA format, got {genome}')
-
- # Add note if genome positional argument has been missed
- if genome is None:
- print('It looks like you missed the Genome FASTA argument, try adding one to the end of your command')
-
- return False
-
- if validate_file(small_rna, 'fastq'):
- small_rna_filetype = 'fastq'
- elif validate_file(small_rna, 'fasta'):
- small_rna_filetype = 'fasta'
- else:
- print(f'Error: expected a small RNA FASTQ or FASTA with at least one sequence, got {small_rna}')
- return False
-
- do_log(quiet, '==> Starting command Sort')
- if not disable_align:
- new_fastq = align_to_genome(genome, small_rna, cds, threads=threads, small_rna_filetype=small_rna_filetype, mismatches=num_mismatches, quiet=quiet)
+ if 'genome' in sort_config:
+ sort_input = align_to_genome(
+ sort_config['genome'], sort_config['input'], sort_config['cds'], program_paths, output,
+ verbose=verbose, keep_intermediate=keep_intermediate, threads=threads, mismatches=sort_config['mismatches']
+ )
else:
- new_fastq = small_rna
-
- table_file = bin_rna_size(new_fastq, min_length, max_length, quiet=quiet)
+ print('====> No genome specified, skipping align...')
+ sort_input = sort_config['input']
- graph_length(table_file)
+ table_file = bin_rna_size(sort_input, sort_config['min_length'], sort_config['max_length'], output, verbose=verbose)
+ graph_length(table_file, output)
- do_log(quiet, '==> Completed command Sort')
-
-def extractnc_command(genome, gff, quiet):
+def run_unitas_on_lenths(unitas_config, min_length, max_length, program_paths, output, threads, verbose, keep_intermediate):
'''
- Code to run when the user chooses to extract the noncoding mRNA region
+ Call label_for_unitas and unitas commands based on the config file
'''
+ REFSEQ_LABELED_DIR = os.path.join(output, 'unitas_refseq_labelled')
+ UNITAS_OUTPUT_DIR = os.path.join(output, 'unitas')
- if not validate_file(genome, 'fasta'):
- print(f'Error: expected a genome in FASTA format, got {genome}')
- return False
-
- if not validate_file(gff, 'gff'):
- print(f'Error: expected a GFF file with at least one feature, got {gff}')
- return False
-
- do_log(quiet, '==> Starting command extractNC')
+ if not os.path.exists(REFSEQ_LABELED_DIR):
+ os.makedirs(REFSEQ_LABELED_DIR)
- extract_noncoding(
- genome, gff,
- quiet=quiet,
- output=os.path.join(get_config_key('general', 'output_directory'), 'noncoding.fasta')
- )
+ if not os.path.exists(UNITAS_OUTPUT_DIR):
+ os.makedirs(UNITAS_OUTPUT_DIR)
- do_log(quiet, '==> Completed command extractNC')
+ # First sort out any refseq files specified
+ refseq_paths = []
-# Parallelism bits for unitas
-# initialize worker processes
-def init_worker(a, b, c, d):
- # declare scope of a new global variable
- global species_name, ref_seqs, quiet, UNITAS_OUTPUT
- # store argument in the global variable for this process
- species_name = a
- ref_seqs = b
- quiet = c
- UNITAS_OUTPUT = d
+ if 'refseq' in unitas_config:
+ for refseq in unitas_config['refseq']:
+ if type(refseq) == str:
+ refseq_paths.append(refseq)
+ elif type(refseq) == dict:
+ for k in refseq.keys():
+ labelled_file = os.path.join(REFSEQ_LABELED_DIR, f'{os.path.basename(refseq[k])}.labelled.fasta')
+ label_file_for_unitas(refseq[k], labelled_file, k)
+ refseq_paths.append(labelled_file)
-# easiest way to implement this quickly
-def unitas_threads(small_rna):
- run_unitas_annotation(small_rna, species_name, ref_seqs, quiet=quiet, unitas_output=UNITAS_OUTPUT)
+ # build and run unitas functions in parrellel
+ small_rna_files = glob.glob(os.path.join(unitas_config['inputs'], '*.fq')) + glob.glob(os.path.join(unitas_config['inputs'], '*.fastq'))
-def unitas_command(small_rna_path, species_name, ref_seqs, cds, unspliced_transcriptome, threads, quiet):
- '''
- Code to run when the user chooses the unitas command
- '''
-
- if not validate_file(small_rna_path, 'directory'):
- print(f'Error: expected a directory containing small RNA FASTQ of varing lengths, got {small_rna_path}')
- return False
-
- if cds is not None:
- new_path = copy_and_label_file(cds, 'CDS', 'Gene')
- if ref_seqs is None:
- ref_seqs = [new_path]
- else:
- ref_seqs.append(new_path)
-
- if unspliced_transcriptome is not None:
- new_path = copy_and_label_file(unspliced_transcriptome, 'UnsplicedTranscriptome', 'Gene')
- if ref_seqs is None:
- ref_seqs = [new_path]
- else:
- ref_seqs.append(new_path)
+ unitas_mappable_function = partial(run_unitas_annotation, unitas_config['species'], refseq_paths, program_paths, UNITAS_OUTPUT_DIR, verbose)
+ with Pool(threads) as p:
+ p.map(unitas_mappable_function, small_rna_files)
- if species_name == 'x' and ref_seqs is None:
- print(f'Error: expected at least one of, a non x species name (-s) or at least one reference file (-r)')
- return False
+ unitas_length_results = merge_summary(output)
+ graph_unitas_classification_type(unitas_length_results, min_length, max_length, output)
- do_log(quiet, '==> Starting command Unitas')
- UNITAS_OUTPUT = os.path.join(get_config_key('general', 'output_directory'), 'unitas')
+ if not keep_intermediate:
+ shutil.rmtree(os.path.join(output, 'unitas_refseq_labelled'))
- mkdir_if_not_exists(UNITAS_OUTPUT)
-
- small_rna_list = glob.glob(os.path.join(small_rna_path, '*.fastq'))
-
- with Pool(threads, initializer=init_worker, initargs=(species_name, ref_seqs, quiet, UNITAS_OUTPUT,)) as p:
- p.map(unitas_threads, small_rna_list)
-
- table_path = merge_summary()
- graph_unitas_classification_type(table_path)
- do_log(quiet, '==> Completed command Unitas')
-
-def targetid_command(small_rna, targets, min_seq_length, mismatches_allowed, threads, quiet):
+def run_targetid(targetid_config, program_paths, output, threads, verbose, keep_intermediate):
'''
- Code to run when the user chooses the targetid command
+ Call the targetid functions based on the settings in the config file
'''
+ if os.path.isdir(targetid_config['inputs']):
+ small_rna_input = glob.glob(os.path.join(targetid_config['inputs'], '*.fq')) + glob.glob(os.path.join(targetid_config['inputs'], '*.fastq'))
+ else:
+ small_rna_input = [targetid_config['inputs']]
- if not validate_file(small_rna, 'fastq'):
- print(f'Error: expected a small RNA FASTQ with at least one sequence, got {small_rna}')
- return False
-
- if targets is None:
- print(f'Error: expected at least one target file (-t)')
- return False
-
- do_log(quiet, '==> Starting TargetID command')
-
- if targets is None:
- print('Error: You need to supply at least one target file with -t')
+ revcomp_file = merge_and_revcomp_input_files(small_rna_input, output)
+ target_sam_files = find_targets(
+ revcomp_file, targetid_config['target_files'], output, program_paths, threads=threads,
+ min_seq_length=targetid_config['min_seq_length'], mismatches_allowed=targetid_config['mismatches'],
+ verbose=verbose
+ )
+ targets = build_summary_files(target_sam_files, output, program_paths, verbose=verbose)
- revcomp_file = revcomp_input_file(small_rna, quiet=quiet)
- sam_files = find_targets(revcomp_file, targets, threads=threads, min_seq_length=min_seq_length, mismatches_allowed=mismatches_allowed, quiet=quiet)
- build_summary_files(sam_files, quiet=quiet)
+ if 'enrich' in targetid_config:
+ do_term_enrichment(
+ targetid_config['enrich']['included_files'], targets, targetid_config['enrich']['eggnog_data'], output, program_paths,
+ threads=threads, verbose=verbose
+ )
- do_log(quiet, '==> Ending TargetID command')
+ if not keep_intermediate:
+ shutil.rmtree(os.path.join(output, 'bowtie_indexes'))
+ shutil.rmtree(os.path.join(output, 'target_alignments'))
+ os.remove(os.path.join(output, 'revcomp_rna.fastq'))
def main():
- parser = ArgumentParser(description='Pipeline to process small RNAs')
- parser.add_argument('-q', '--quiet', help='Supress output from intermediate commands', action='count', default=0)
- parser.add_argument('-C', '--path-to-config', help='Path to the TOML format config file to use', default='config.toml')
+ parser = ArgumentParser(description='Pipeline to process small RNAs - version 2')
- subparsers = parser.add_subparsers(dest='command')
+ parser.add_argument('config_file', help='Path to YAML config file')
- parser_process = subparsers.add_parser('process', help='Preprocessing for the RNA')
- parser_process.add_argument('-a', '--adapter', help='Sequence of the adapter to remove from the 3\' end')
- parser_process.add_argument('-g', '--front', help='Sequence of the adapter to remove from the 5\' end')
- parser_process.add_argument('-b', '--anywhere', help='Sequence of the adapters to remove from both ends')
- parser_process.add_argument('-c', '--cutoff', help='Quality cutoff to trin RNA sequences at', default=20, type=int)
- parser_process.add_argument('small_rna', help='Path to FASTQ containing the small RNA')
-
- parser_sort = subparsers.add_parser('sort', help='Find RNAs that align to a genome and sort them by length')
- parser_sort.add_argument('-d', '--cds', help='Optional CDS region, also align this to the CDS region as well as the genome')
- parser_sort.add_argument('-l', '--min-length', help='Minimum length to bin', type=int, default=-inf)
- parser_sort.add_argument('-x', '--max-length', help='Maximum length to bin', type=int, default=inf)
- parser_sort.add_argument('-m', '--ref-mismatches', type=int, default=None, help='Number of mismatches to use in bowtie2, None for default behaviour')
- parser_sort.add_argument('--disable-alignment', action='store_true', help='Skip the alignment to the reference genome step')
- parser_sort.add_argument('small_rna', help='Path to FASTQ containing the small RNA')
- parser_sort.add_argument('genome', nargs='?', default=None, help='Genome to align against')
-
- parser_extractnc = subparsers.add_parser('extractnc', help='Extarct the noncoding region from a fasta with a GFF file')
- parser_extractnc.add_argument('genome', help='FASTA containing the genome to extract from')
- parser_extractnc.add_argument('gff_file', help='GFF file containing annotations of CDS and mRNA regions')
-
- parser_unitas = subparsers.add_parser('unitas', help='Run unitas on split files and merge results')
- parser_unitas.add_argument('-d', '--cds', help='Optional CDS region, passed to unitas')
- parser_unitas.add_argument('-u', '--unspliced-transcriptome', help='Optional, unspliced transcriptome, passed to unitas')
- parser_unitas.add_argument('-r', '--refseq', help='References for use with unitas', nargs='*', default=None)
- parser_unitas.add_argument('-s', '--species', help='Species to set in unitas arguments', default='x')
- parser_unitas.add_argument('path_to_rnas', help='Path to the folder with varying length RNAs in')
-
- parser_targetid = subparsers.add_parser('targetid', help='Align small RNA to a number of genome features to find out what is targeted')
- parser_targetid.add_argument('-m', '--min-seq-length', help='Minimum sequence length to properly align', default=5)
- parser_targetid.add_argument('-t', '--target-files', help='Files containing genome features that could be targeted', nargs='+')
- parser_targetid.add_argument('--num-mismatches', help='Number of mismatches to allow in the alignment, defaults to 0', default=0, type=int)
- parser_targetid.add_argument('small_rna', help='Path to the FASTQ containing the small RNA to find targets of')
-
- parser_all = subparsers.add_parser('all', help='Run process, sort and unitas one after the other')
- parser_all.add_argument('-a', '--adapter', help='Sequence of the adapter to remove from the 3\' end')
- parser_all.add_argument('-d', '--cds', help='Optional CDS region, also align this to the CDS region as well as the genome')
- parser_all.add_argument('-g', '--front', help='Sequence of the adapter to remove from the 5\' end')
- parser_all.add_argument('-b', '--anywhere', help='Sequence of the adapters to remove from both ends')
- parser_all.add_argument('-c', '--cutoff', help='Quality cutoff to trin RNA sequences at', default=20, type=int)
- parser_all.add_argument('-l', '--min-length', help='Minimum length to bin', type=int, default=-inf)
- parser_all.add_argument('-x', '--max-length', help='Maximum length to bin', type=int, default=inf)
- parser_all.add_argument('-r', '--refseq', help='References for use with unitas', nargs='*', default=None)
- parser_all.add_argument('-s', '--species', help='Species to set in unitas arguments', default='x')
- parser_all.add_argument('-u', '--unspliced-transcriptome', help='Optional, unspliced transcriptome, passed to unitas')
- parser_all.add_argument('-m', '--ref-mismatches', type=int, default=None, help='Number of mismatches to use in bowtie2 when aligning to the genome, None for default behaviour')
- parser_all.add_argument('--disable-alignment', action='store_true', help='Skip the alignment to the reference genome step')
- parser_all.add_argument('small_rna', help='Path to FASTQ containing the small RNA')
- parser_all.add_argument('genome', nargs='?', default=None, help='Genome to align against')
+ parser.add_argument('-o', '--output', help='Directory to write output to (will be created if it doesn\'t exist, default: hlsmallrna_output)', default='hlsmallrna_output')
+ parser.add_argument('-t', '--threads', help='Number of threads to run the pipeline on (default: 4)', type=int, default=4)
+ parser.add_argument('-k', '--keep-files', help='If set, keep the intermediate files, to help debug the application', action='store_true')
+ parser.add_argument('-v', '--verbose', help='If set, print the output for the intermediate commands', action='store_true')
args = parser.parse_args()
- do_log(args.quiet, f"Output directory set to {os.path.abspath(get_config_key('general', 'output_directory'))}")
- if os.path.exists(get_config_key('general', 'output_directory')):
- do_log(args.quiet, 'WARNING: output directory already exists, this run may use previously generated data, if you didn\'t mean this, try again with a different output directory')
-
- load_config(args.path_to_config, quiet=args.quiet)
-
- def get_command_args(name):
- arguments = vars(args)
-
- if name in arguments and arguments[name] is not None:
- return arguments[name]
- else:
- try:
- return get_config_key('command', name)
- except KeyError:
- return None
-
- mkdir_if_not_exists(get_config_key('general', 'output_directory'))
- num_threads = get_config_key('general', 'threads')
-
- if args.command == 'process':
- process_command(
- get_command_args('small_rna'),
- get_command_args('adapter'),
- get_command_args('front'),
- get_command_args('anywhere'),
- get_command_args('cutoff'),
- get_command_args('quiet')
- )
-
- if args.command == 'sort':
- sort_command(
- get_command_args('genome'),
- get_command_args('small_rna'),
- get_command_args('cds'),
- get_command_args('min_length'),
- get_command_args('max_length'),
- get_command_args('ref_mismatches'),
- get_command_args('disable_alignment'),
- num_threads,
- get_command_args('quiet')
- )
-
- if args.command == 'extractnc':
- extractnc_command(
- get_command_args('genome'),
- get_command_args('gff_file'),
- get_command_args('quiet')
- )
-
- if args.command == 'unitas':
- unitas_command(
- get_command_args('path_to_rnas'),
- get_command_args('species'),
- get_command_args('refseq'),
- get_command_args('cds'),
- get_command_args('unspliced_transcriptome'),
- num_threads,
- get_command_args('quiet')
- )
-
- if args.command == 'targetid':
- targetid_command(
- get_command_args('small_rna'),
- get_command_args('target_files'),
- get_command_args('min_seq_length'),
- get_command_args('num_mismatches'),
- num_threads,
- get_command_args('quiet')
- )
-
- if args.command == 'all':
- out_code = process_command(
- get_command_args('small_rna'),
- get_command_args('adapter'),
- get_command_args('front'),
- get_command_args('anywhere'),
- get_command_args('cutoff'),
- get_command_args('quiet')
- )
-
- if out_code is not None:
- return
+ if not os.path.exists(args.output):
+ os.makedirs(args.output)
- out_code = sort_command(
- get_command_args('genome'),
- os.path.join(get_config_key('general', 'output_directory'), 'cut_sequences.fastq'),
- get_command_args('cds'),
- get_command_args('min_length'),
- get_command_args('max_length'),
- get_command_args('ref_mismatches'),
- get_command_args('disable_alignment'),
- num_threads,
- get_command_args('quiet')
- )
+ shutil.copy2(args.config_file, os.path.join(args.output, 'config.yml'))
- if out_code is not None:
- return
-
- unitas_command(
- os.path.join(get_config_key('general', 'output_directory'), 'binned_rna'),
- get_command_args('species'),
- get_command_args('refseq'),
- get_command_args('cds'),
- get_command_args('unspliced_transcriptome'),
- num_threads,
- get_command_args('quiet')
- )
+ print('==> Validating config and converting files...')
+ success, err_msg, config = load_and_validate_input(args.config_file, args.output)
-def main_ssoverlap():
- '''
- Main function form the overlap_ss script
- '''
+ if not success:
+ print(f'Error: {err_msg}')
+ exit(1)
- parser = ArgumentParser(description='Looks for overlaps of two classes of RNA on the same strand of DNA')
-
- parser.add_argument('genome', help='FASTA file containing the genome of the organism')
- parser.add_argument('rna_file_1', help='File containing the RNA to use as a base')
- parser.add_argument('rna_file_2', help='File containing the RNA to look for overlaps with')
-
- parser.add_argument('-q', '--quiet', help='Supress output from intermediate commands', action='count', default=0)
-
- args = parser.parse_args()
+ if 'trim' in config:
+ print('==> Trimming Reads...')
+ run_trim_reads(config['trim'], config['program_paths'], args.output, args.threads, args.verbose, args.keep_files)
- OUTPUT_DIR = os.path.join(get_config_key('general', 'output_directory'), 'samestrand_overlap')
-
- mkdir_if_not_exists(OUTPUT_DIR)
-
- do_log(args.quiet, '==> Converting RNA to FASTA format')
-
- global longest_rna
- longest_rna = 0
- for rna_file in [args.rna_file_1, args.rna_file_2]:
- if (
- rna_file.endswith('.fa')
- or rna_file.endswith('.fasta')
- or rna_file.endswith('.fna')
- or rna_file.endswith('.ffn')
- or rna_file.endswith('.faa')
- or rna_file.endswith('.frn')
- ):
- filetype = 'fasta'
- elif rna_file.endswith('.fq') or rna_file.endswith('.fastq'):
- filetype = 'fastq'
-
- def convert_to_fasta(seq):
- global longest_rna
- if len(seq) > longest_rna:
- longest_rna = len(seq)
-
- return seq
-
- seqs = SeqIO.parse(rna_file, filetype)
- seqs = map(convert_to_fasta, seqs)
- SeqIO.write(seqs, os.path.join(OUTPUT_DIR, os.path.basename(rna_file) + '.fasta'), 'fasta')
-
- samestrand_overlap(
- args.genome,
- os.path.join(OUTPUT_DIR, os.path.basename(args.rna_file_1) + '.fasta'),
- os.path.join(OUTPUT_DIR, os.path.basename(args.rna_file_2) + '.fasta'),
- longest_rna,
- quiet=args.quiet
- )
+ if 'sort' in config:
+ print('==> Aligning to genome and sorting reads by length...')
+ run_sort_reads(config['sort'], config['program_paths'], args.output, args.threads, args.verbose, args.keep_files)
+
+ if 'unitas' in config:
+ if 'sort' in config:
+ min_length = config['sort']['min_length']
+ max_length = config['sort']['max_length']
+ else:
+ min_length, max_length = find_min_max_lengths(config['unitas']['inputs'])
+ print('==> Running Unitas in parellel on each sequence length...')
+ run_unitas_on_lenths(config['unitas'], min_length, max_length, config['program_paths'], args.output, args.threads, args.verbose, args.keep_files)
-if __name__ == '__main__':
- main()
\ No newline at end of file
+ if 'targetid' in config:
+ print('==> Finding targets of small RNA...')
+ run_targetid(config['targetid'], config['program_paths'], args.output, args.threads, args.verbose, args.keep_files)
\ No newline at end of file
diff --git a/bin/build_coord_files b/hlsmallrna/bin/build_coord_files.py
similarity index 99%
rename from bin/build_coord_files
rename to hlsmallrna/bin/build_coord_files.py
index fc4dae4..672bb1f 100755
--- a/bin/build_coord_files
+++ b/hlsmallrna/bin/build_coord_files.py
@@ -1,5 +1,5 @@
#!/usr/bin/env python3
-# Copyright 2022 Vicky Hunt Lab Members
+# Copyright 2022 - 2025 Vicky Hunt Lab Members
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
@@ -280,7 +280,7 @@ def process_fa_out(path_to_fa_out, output, feature, scaffold_index=None):
writer.writerow(['Chrom', 'Name', 'Start', 'End'])
writer.writerows(antisense_targets)
-if __name__ == '__main__':
+def main():
parser = ArgumentParser(description='Convert common bioinformatics file formats into coordinate file to plot')
parser.add_argument('-a', '--fasta', action='store_true', help='Treat input as a FASTA file')
diff --git a/hlsmallrna/create_noncoding.py b/hlsmallrna/bin/extract_nc.py
old mode 100644
new mode 100755
similarity index 86%
rename from hlsmallrna/create_noncoding.py
rename to hlsmallrna/bin/extract_nc.py
index 86d880a..56b9979
--- a/hlsmallrna/create_noncoding.py
+++ b/hlsmallrna/bin/extract_nc.py
@@ -1,4 +1,5 @@
-# Copyright 2022 Vicky Hunt Lab Members
+#!/usr/bin/env python3
+# Copyright 2022 - 2025 Vicky Hunt Lab Members
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
@@ -16,8 +17,7 @@
from sys import argv
from collections import defaultdict
-
-from .config import do_log
+from argparse import ArgumentParser
def get_pairs(array):
'''
@@ -115,14 +115,14 @@ def merge_cds(coding_reigon):
return n_merged
-def extract_noncoding(genome, gff_path, quiet=0, output='result.fasta'):
+def extract_noncoding(genome, gff_path, output='result.fasta'):
'''
Extract the noncoding reigon from the genome, basied on a GFF file
'''
gff_iter = DataIterator(gff_path)
genome_data = SeqIO.parse(genome, 'fasta')
- do_log(quiet, '====> Calculating coordinates')
+ print('====> Calculating coordinates')
coordinates = defaultdict(lambda: [])
mRNAs = defaultdict(lambda: [])
@@ -141,13 +141,13 @@ def extract_noncoding(genome, gff_path, quiet=0, output='result.fasta'):
if item[2] == 'CDS':
coding_reigon[item[0] + item[6]].append([int(item[3]), int(item[4])])
- do_log(quiet, '====> Merging and validating coordinates')
+ print('====> Merging and validating coordinates')
cds_merged = merge_cds(coding_reigon)
mRNA_merged = merge_cds(mRNAs)
validate_gff(mRNAs, coding_reigon)
- do_log(quiet, f'Merged {cds_merged} coding reigons and {mRNA_merged} mRNAs')
+ print(f'Merged {cds_merged} coding reigons and {mRNA_merged} mRNAs')
for key in mRNAs.keys():
for item in mRNAs[key]:
@@ -159,12 +159,21 @@ def extract_noncoding(genome, gff_path, quiet=0, output='result.fasta'):
coordinates[key].append(item[0])
coordinates[key].append(item[1])
- do_log(quiet, '====> Extarcting fragments')
+ print('====> Extarcting fragments')
for key in coordinates.keys():
coordinates[key].sort()
SeqIO.write(extract_fragments(genome_data, coordinates, mRNA_start, mRNA_end), output, 'fasta')
-if __name__ == '__main__':
- extract_noncoding(argv[1], argv[2])
\ No newline at end of file
+def main():
+ parser = ArgumentParser('')
+
+ parser.add_argument('genome', help='FASTA containing the genome to extract from')
+ parser.add_argument('gff_file', help='GFF file containing annotations of CDS and mRNA regions')
+
+ parser.add_argument('-o', '--output', help='FASTA file to write output to', default='result.fasta')
+
+ args = parser.parse_args()
+
+ extract_noncoding(args.genome, args.gff_file, output=args.output)
\ No newline at end of file
diff --git a/bin/revcomp_rna b/hlsmallrna/bin/revcomp_rna.py
old mode 100644
new mode 100755
similarity index 57%
rename from bin/revcomp_rna
rename to hlsmallrna/bin/revcomp_rna.py
index 2651b39..75eed45
--- a/bin/revcomp_rna
+++ b/hlsmallrna/bin/revcomp_rna.py
@@ -1,5 +1,5 @@
#!/usr/bin/env python3
-# Copyright 2022 Vicky Hunt Lab Members
+# Copyright 2022 - 2025 Vicky Hunt Lab Members
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
@@ -17,16 +17,6 @@
from Bio import SeqIO
-parser = ArgumentParser(description='Simple script to reverse complement a FASTA file')
-
-parser.add_argument('-i', '--input', help='FASTA file to reverse complement')
-parser.add_argument('-o', '--output', help='FASTA file to write to')
-
-args = parser.parse_args()
-
-if args.input is None or args.output is None:
- raise Exception('Both -i and -o need to be set to run this')
-
def reverse_complement(seq):
revcomp = seq.reverse_complement()
revcomp.id = seq.id
@@ -34,7 +24,18 @@ def reverse_complement(seq):
return revcomp
-seqs = SeqIO.parse(args.input, 'fasta')
-seqs = map(reverse_complement, seqs)
+def main():
+ parser = ArgumentParser(description='Simple script to reverse complement a FASTA file')
+
+ parser.add_argument('-i', '--input', help='FASTA file to reverse complement')
+ parser.add_argument('-o', '--output', help='FASTA file to write to')
+
+ args = parser.parse_args()
+
+ if args.input is None or args.output is None:
+ raise Exception('Both -i and -o need to be set to run this')
+
+ seqs = SeqIO.parse(args.input, 'fasta')
+ seqs = map(reverse_complement, seqs)
-SeqIO.write(seqs, args.output, 'fasta')
\ No newline at end of file
+ SeqIO.write(seqs, args.output, 'fasta')
diff --git a/hlsmallrna/config.py b/hlsmallrna/config.py
deleted file mode 100644
index aa28cb3..0000000
--- a/hlsmallrna/config.py
+++ /dev/null
@@ -1,123 +0,0 @@
-import toml
-
-import os
-
-from multiprocessing import cpu_count
-
-# global object to load config into, it works, but isn't nice
-# or good practice
-CONFIG = {}
-
-DEFALT_CONFIG = {
- 'general': {
- 'threads': cpu_count(),
- 'output_directory': './output'
- },
- 'cli-tools': {
- 'trim': {
- 'trim_pass_threads': True,
- 'path_to_trim': 'cutadapt',
- 'trim_type': '',
- 'trim_params': []
- },
- 'fastqc': {
- 'fastqc_pass_threads': True,
- 'path_to_fastqc': 'fastqc',
- 'fastqc_params': []
- },
- 'bbmap': {
- 'bbmap_pass_threads': True,
- 'path_to_bbmap': 'bbmap.sh',
- 'bbmap_align_params': [],
- 'bbmap_index_params': []
- },
- 'unitas': {
- 'unitas_pass_threads': True,
- 'path_to_unitas': 'unitas.pl',
- 'unitas_params': []
- },
- 'bowtie2': {
- 'bowtie2_pass_threads': True,
- 'path_to_bowtie2': 'bowtie2',
- 'bowtie2_params': [],
- 'path_to_bowtie2_build': 'bowtie2-build',
- 'bowtie2_build_params': []
- },
- 'samtools': {
- 'path_to_samtools': 'samtools',
- 'samtools_view_params': ['-h', '-F', '256', '-F', '4'],
- 'samtools_fastq_params': []
- },
- 'bedtools': {
- 'path_to_bedtools': 'bedtools',
- 'bedtools_bamtofastq_params': []
- }
- }
-}
-
-def load_config(path_to_config, quiet=0):
- '''
- Load the config file into the global CONFIG varible and validate the result
- '''
- global CONFIG
-
- try:
- CONFIG = toml.load(path_to_config)
-
- result = validate_toml(CONFIG, DEFALT_CONFIG)
-
- if result is not None:
- raise Exception(result)
- except FileNotFoundError:
- do_log(quiet, f'No {path_to_config} found, running with defaults...')
-
-def validate_toml(in_dict, compare):
- '''
- Check the type in a TOML dict (in_dict) match the types of another dict,
- validates the types in a loaded config file
- '''
- for key in compare:
- if key in in_dict:
- if type(in_dict[key]) != type(compare[key]):
- print('Error in config file!')
-
- return key + ' should be ' + str(type(compare[key])) + ' but found ' + str(type(in_dict[key]))
-
- if type(in_dict[key]) == dict:
- result = validate_toml(in_dict[key], compare[key])
-
- if result is not None:
- return result
-
-def get_config_key(*keys):
- '''
- Retrieve a config key from the config file if it exists or the
- default config if it doesn't
- '''
- try:
- last_item = CONFIG
-
- for key in keys:
- last_item = last_item[key]
- except KeyError:
- last_item = DEFALT_CONFIG
-
- for key in keys:
- last_item = last_item[key]
-
- return last_item
-
-def mkdir_if_not_exists(path):
- '''
- Shorthand to make directories only if they don't exist already
- without crashing
- '''
- if not os.path.exists(path):
- os.mkdir(path)
-
-def do_log(quiet, *print_args, **print_kwds):
- '''
- Function that prints if quiet isn't set
- '''
- if quiet < 2:
- print(*print_args, **print_kwds)
\ No newline at end of file
diff --git a/hlsmallrna/fastqc.py b/hlsmallrna/fastqc.py
deleted file mode 100644
index 1bb674f..0000000
--- a/hlsmallrna/fastqc.py
+++ /dev/null
@@ -1,101 +0,0 @@
-# Copyright 2022 Vicky Hunt Lab Members
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-# http://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software
-# distributed under the License is distributed on an "AS IS" BASIS,
-# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-# See the License for the specific language governing permissions and
-# limitations under the License.
-import os.path
-import zipfile
-import glob
-import csv
-
-from subprocess import run
-from io import StringIO
-from math import inf
-
-from Bio import SeqIO
-
-from .config import get_config_key, mkdir_if_not_exists, do_log
-
-def run_fastqc(path_to_fastq, quiet=0):
- '''
- Call fastqc and trim off any bases with a lower quatile below 20
- '''
- outdir = os.path.join(get_config_key('general', 'output_directory'), 'fastqc')
- mkdir_if_not_exists(outdir)
-
- command = [
- get_config_key('cli-tools', 'fastqc', 'path_to_fastqc'),
- '--outdir', outdir,
- path_to_fastq
- ]
-
- command = command + get_config_key('cli-tools', 'fastqc', 'fastqc_params')
-
- if get_config_key('cli-tools', 'fastqc', 'fastqc_pass_threads'):
- threads = get_config_key('general', 'threads')
-
- command.append('-t')
- command.append(str(threads))
-
- do_log(quiet, '====> Running FastQC')
- result = run(command, capture_output=(quiet != 0))
- result.check_returncode()
-
- return result.returncode == 0
-
-def cut_rna_below_cutoff(fastq_path, cutoff, quiet=0):
- '''
- Read the fastqc result and automatically trim off any base that is below the threshold
- '''
-
- do_log(quiet, '====> Extracting results from FastQC Output')
-
- # Unzip file produced by fastqc
- path_to_zip = os.path.join(get_config_key('general', 'output_directory'), 'fastqc')
- zip_file = glob.glob(os.path.join(path_to_zip, '*.zip'))[0]
-
- with zipfile.ZipFile(zip_file, 'r') as f:
- f.extractall(path=path_to_zip)
-
- # parse in quality score data, cut out part of file
- with open('.' + os.path.join(zip_file.strip('.zip'), 'fastqc_data.txt')) as f:
- data = f.readlines()
-
- start_line = None
- end_line = None
- for i, line in enumerate(data):
- if line.startswith('>>Per base sequence quality'):
- start_line = i + 2
-
- if start_line is not None and line.startswith('>>END_MODULE'):
- end_line = i
- break
-
- # use csv reader for simplicity
- cutoff_base = inf
-
- reader = reversed(list(csv.reader(StringIO(''.join(data[start_line:end_line])), delimiter='\t')))
- for line in reader:
- if float(line[3]) < cutoff:
- if cutoff_base == int(line[0]) + 1:
- cutoff_base = int(line[0])
- elif cutoff_base == inf:
- if float(line[3]) < cutoff:
- cutoff_base = int(line[0])
- else:
- cutoff_base = None
-
- # cut sequences from 3' end
- cut_func = lambda x: x[:cutoff_base]
-
- sequences = SeqIO.parse(fastq_path, 'fastq')
- cut_seqs = map(cut_func, sequences)
- SeqIO.write(cut_seqs, os.path.join(get_config_key('general', 'output_directory'), 'cut_sequences.fastq'), 'fastq')
\ No newline at end of file
diff --git a/hlsmallrna/genome_align.py b/hlsmallrna/genome_align.py
index 3fc04c1..bc6ae29 100644
--- a/hlsmallrna/genome_align.py
+++ b/hlsmallrna/genome_align.py
@@ -1,4 +1,4 @@
-# Copyright 2022 Vicky Hunt Lab Members
+# Copyright 2022 - 2025 Vicky Hunt Lab Members
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
@@ -11,20 +11,21 @@
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
-from collections import defaultdict
+from shutil import copy2
+
import os
import os.path
import csv
+import shutil
from subprocess import run
+from collections import defaultdict
import numpy as np
from Bio import SeqIO, SeqRecord, Seq
from matplotlib import pyplot as plt
-from .config import get_config_key, mkdir_if_not_exists, do_log
-
def into_seqrecord(seqs):
'''
Convert a set of Biopython sequences into SeqRecord objects
@@ -35,130 +36,58 @@ def into_seqrecord(seqs):
yield record
-def make_fastqs_unique(fastq1, fastq2, output):
+def merge_genome_cds_hits(results_seq, cds_seq):
'''
- Take 2 FastQ files, combine and remove duplicates from the second
+ Merge FASTQ hits from the genome and CDS
'''
- unique_sequences = []
- cds_sequences = defaultdict(lambda: 0)
-
- for sequence in SeqIO.parse(fastq1, 'fastq'):
- unique_sequences.append(str(sequence.seq))
+ for seq in SeqIO.parse(results_seq, 'fastq'):
+ yield seq.seq
- try:
- for sequence in SeqIO.parse(fastq2, 'fastq'):
- cds_sequences[str(sequence.seq)] += 1
- except FileNotFoundError:
- pass
-
- for key in cds_sequences.keys():
- if key not in unique_sequences:
- for i in range(cds_sequences[key]):
- unique_sequences.append(key)
+ if os.path.exists(cds_seq):
+ for seq in SeqIO.parse(cds_seq, 'fastq'):
+ yield seq.seq
- SeqIO.write(into_seqrecord(unique_sequences), output, 'fastq')
-
-def make_fastq_overlap_only(fastq1, fastq2, output):
- '''
- Create a fastq with only sequences that appear in both fastq files
- '''
- fastq1_seqs = []
-
- for sequence in SeqIO.parse(fastq1, 'fastq'):
- fastq1_seqs.append(str(sequence.seq))
-
- try:
- fastq2_seqs = defaultdict(lambda: 0)
-
- for sequence in SeqIO.parse(fastq2, 'fastq'):
- fastq2_seqs[str(sequence.seq)] += 1
-
- for seq in fastq2_seqs.keys():
- if seq not in fastq1_seqs:
- for i in range(fastq2_seqs[seq]):
- fastq1_seqs.append(seq)
-
- result_seqs = fastq1_seqs
- except FileNotFoundError as e:
- result_seqs = fastq1_seqs
-
- SeqIO.write(into_seqrecord(result_seqs), output, 'fastq')
-
-def remove_symbols_from_header(fasta):
- '''
- Remove @s from the FASTA header before doing the alignment
- '''
- for seq in SeqIO.parse(fasta, 'fasta'):
- seq.id = seq.id.replace('@', '_')
-
- yield seq
-
-def align_to_genome(genome, small_rnas, cds, quiet=0, threads=4, small_rna_filetype='fastq', mismatches=None):
+def align_to_genome(genome, small_rnas, cds, program_paths, output, verbose=False, keep_intermediate=False, threads=4, mismatches=0):
'''
Align the small RNAs to the genome and filter out any that are unsuccessful
'''
- INTERMEDIATE_SAM = os.path.join(get_config_key('general', 'output_directory'), 'mapped_sequences.sam')
- INTERMEDIATE_BAM = os.path.join(get_config_key('general', 'output_directory'), 'mapped_sequences.bam')
- UNMAPPED_BAM = os.path.join(get_config_key('general', 'output_directory'), 'unmapped_sequences.bam')
- RESULT_FASTQ = os.path.join(get_config_key('general', 'output_directory'), 'genome_mapped_sequences.fastq')
- RESULT_UNMAPPED_FASTQ = os.path.join(get_config_key('general', 'output_directory'), 'genome_unmapped_sequences.fastq')
- INDEX_DIRECTORY = os.path.join(get_config_key('general', 'output_directory'), 'bbmap_index')
-
- CDS_INTERMEDIATE_SAM = os.path.join(get_config_key('general', 'output_directory'), 'cds_sequences.sam')
- CDS_INTERMEDIATE_BAM = os.path.join(get_config_key('general', 'output_directory'), 'cds_sequences.bam')
- CDS_UNMAPPED_BAM = os.path.join(get_config_key('general', 'output_directory'), 'cds_unmapped_sequences.bam')
- CDS_UNMAPPED_FASTQ = os.path.join(get_config_key('general', 'output_directory'), 'cds_unmapped_sequences.fastq')
- CDS_RESULT_FASTQ = os.path.join(get_config_key('general', 'output_directory'), 'cds_sequences.fastq')
-
- FINAL_FASTQ = os.path.join(get_config_key('general', 'output_directory'), 'mapped_sequences.fastq')
- FINAL_UNMAPPED_FASTQ = os.path.join(get_config_key('general', 'output_directory'), 'unmapped_sequences.fastq')
+ INTERMEDIATE_SAM = os.path.join(output, 'mapped_sequences.sam')
+ INTERMEDIATE_BAM = os.path.join(output, 'mapped_sequences.bam')
+ UNMAPPED_BAM = os.path.join(output, 'unmapped_sequences.bam')
+ RESULT_FASTQ = os.path.join(output, 'genome_mapped_sequences.fastq')
+ RESULT_UNMAPPED_FASTQ = os.path.join(output, 'genome_unmapped_sequences.fastq')
+ INDEX_DIRECTORY = os.path.join(output, 'bowtie2_index')
- mkdir_if_not_exists(INDEX_DIRECTORY)
+ CDS_INTERMEDIATE_SAM = os.path.join(output, 'cds_sequences.sam')
+ CDS_INTERMEDIATE_BAM = os.path.join(output, 'cds_sequences.bam')
+ CDS_UNMAPPED_BAM = os.path.join(output, 'cds_unmapped_sequences.bam')
+ CDS_UNMAPPED_FASTQ = os.path.join(output, 'cds_unmapped_sequences.fastq')
+ CDS_RESULT_FASTQ = os.path.join(output, 'cds_sequences.fastq')
- if small_rna_filetype == 'fasta':
- NEW_SMALL_RNAS = os.path.join(get_config_key('general', 'output_directory'), 'corrected_headers.fasta')
- SeqIO.write(remove_symbols_from_header(small_rnas), NEW_SMALL_RNAS, 'fasta')
+ FINAL_FASTQ = os.path.join(output, 'mapped_sequences.fastq')
+ FINAL_UNMAPPED_FASTQ = os.path.join(output, 'unmapped_sequences.fastq')
- small_rnas = NEW_SMALL_RNAS
+ if not os.path.exists(INDEX_DIRECTORY):
+ os.makedirs(INDEX_DIRECTORY)
- bbmap_build_index = [
- get_config_key('cli-tools', 'bowtie2', 'path_to_bowtie2_build'),
+ bowtie2_build_index = [
+ program_paths['bowtie2-build'],
'--threads', str(threads),
genome,
os.path.join(INDEX_DIRECTORY, 'genome_index')
]
- bbmap_build_index = bbmap_build_index + get_config_key('cli-tools', 'bowtie2', 'bowtie2_build_params')
-
cds_build_index = [
- get_config_key('cli-tools', 'bowtie2', 'path_to_bowtie2_build'),
+ program_paths['bowtie2-build'],
'--threads', str(threads),
cds,
os.path.join(INDEX_DIRECTORY, 'cds_index')
]
- cds_build_index = cds_build_index + get_config_key('cli-tools', 'bowtie2', 'bowtie2_build_params')
-
- if mismatches is not None:
- bbmap_align_reads = [
- get_config_key('cli-tools', 'bowtie2', 'path_to_bowtie2'),
- '--threads', str(threads),
- '-L', '18',
- '--no-1mm-upfront',
- '--score-min', 'L,' + str(-mismatches) + ',0',
- '--end-to-end',
- '--mp', '1,1',
- '--ignore-quals',
- '--rdg', '9,1',
- '--rfg', '9,1',
- '-x', os.path.join(INDEX_DIRECTORY, 'genome_index'),
- '-U', small_rnas,
- '-S', INTERMEDIATE_SAM
- ]
- elif mismatches == 0:
- bbmap_align_reads = [
- get_config_key('cli-tools', 'bowtie2', 'path_to_bowtie2'),
+ if mismatches == 0:
+ bowtie2_align_reads = [
+ program_paths['bowtie2'],
'--threads', str(threads),
'-L', '18',
'--no-1mm-upfront',
@@ -170,39 +99,24 @@ def align_to_genome(genome, small_rnas, cds, quiet=0, threads=4, small_rna_filet
'-S', INTERMEDIATE_SAM
]
else:
- bbmap_align_reads = [
- get_config_key('cli-tools', 'bowtie2', 'path_to_bowtie2'),
- '--threads', str(threads),
- '-L', '18',
- '-x', os.path.join(INDEX_DIRECTORY, 'genome_index'),
- '-U', small_rnas,
- '-S', INTERMEDIATE_SAM
- ]
-
- if small_rna_filetype == 'fasta':
- bbmap_align_reads.append('-f')
-
- bbmap_align_reads = bbmap_align_reads + get_config_key('cli-tools', 'bowtie2', 'bowtie2_params')
-
- if mismatches is not None:
- cds_align_reads = [
- get_config_key('cli-tools', 'bowtie2', 'path_to_bowtie2'),
+ bowtie2_align_reads = [
+ program_paths['bowtie2'],
'--threads', str(threads),
'-L', '18',
- '--no-1mm-upfront',
'--score-min', 'L,' + str(-mismatches) + ',0',
'--end-to-end',
'--mp', '1,1',
'--ignore-quals',
- '--rdg', '9,1',
- '--rfg', '9,1',
- '-x', os.path.join(INDEX_DIRECTORY, 'cds_index'),
+ '--rdg', '100,100',
+ '--rfg', '100,100',
+ '-x', os.path.join(INDEX_DIRECTORY, 'genome_index'),
'-U', small_rnas,
- '-S', CDS_INTERMEDIATE_SAM
+ '-S', INTERMEDIATE_SAM
]
- elif mismatches == 0:
+
+ if mismatches == 0:
cds_align_reads = [
- get_config_key('cli-tools', 'bowtie2', 'path_to_bowtie2'),
+ program_paths['bowtie2'],
'--threads', str(threads),
'-L', '18',
'--no-1mm-upfront',
@@ -210,38 +124,38 @@ def align_to_genome(genome, small_rnas, cds, quiet=0, threads=4, small_rna_filet
'--end-to-end',
'-M', '0',
'-x', os.path.join(INDEX_DIRECTORY, 'cds_index'),
- '-U', small_rnas,
+ '-U', RESULT_UNMAPPED_FASTQ,
'-S', CDS_INTERMEDIATE_SAM
]
else:
cds_align_reads = [
- get_config_key('cli-tools', 'bowtie2', 'path_to_bowtie2'),
+ program_paths['bowtie2'],
'--threads', str(threads),
'-L', '18',
+ '--score-min', 'L,' + str(-mismatches) + ',0',
+ '--end-to-end',
+ '--mp', '1,1',
+ '--ignore-quals',
+ '--rdg', '100,100',
+ '--rfg', '100,100',
'-x', os.path.join(INDEX_DIRECTORY, 'cds_index'),
- '-U', small_rnas,
+ '-U', RESULT_UNMAPPED_FASTQ,
'-S', CDS_INTERMEDIATE_SAM
]
- if small_rna_filetype == 'fasta':
- cds_align_reads.append('-f')
-
- cds_align_reads = cds_align_reads + get_config_key('cli-tools', 'bowtie2', 'bowtie2_params')
-
samtools_view = [
- get_config_key('cli-tools', 'samtools', 'path_to_samtools'),
+ program_paths['samtools'],
'view',
'-h',
'-b',
'-F', '4',
+ '-F', '256',
'-o', INTERMEDIATE_BAM,
- INTERMEDIATE_SAM
+ INTERMEDIATE_SAM,
]
- samtools_view = samtools_view + get_config_key('cli-tools', 'samtools', 'samtools_view_params')
-
bedtools_bamtofastq_command = [
- get_config_key('cli-tools', 'samtools', 'path_to_samtools'),
+ program_paths['samtools'],
'fastq',
INTERMEDIATE_BAM,
'-o', RESULT_FASTQ,
@@ -249,19 +163,18 @@ def align_to_genome(genome, small_rnas, cds, quiet=0, threads=4, small_rna_filet
]
cds_samtools_view = [
- get_config_key('cli-tools', 'samtools', 'path_to_samtools'),
+ program_paths['samtools'],
'view',
'-h',
'-b',
'-F', '4',
+ '-F', '256',
'-o', CDS_INTERMEDIATE_BAM,
CDS_INTERMEDIATE_SAM
]
- cds_samtools_view = cds_samtools_view + get_config_key('cli-tools', 'samtools', 'samtools_view_params')
-
cds_bedtools_bamtofastq_command = [
- get_config_key('cli-tools', 'samtools', 'path_to_samtools'),
+ program_paths['samtools'],
'fastq',
CDS_INTERMEDIATE_BAM,
'-o', CDS_RESULT_FASTQ,
@@ -269,7 +182,7 @@ def align_to_genome(genome, small_rnas, cds, quiet=0, threads=4, small_rna_filet
]
unmapped_samtools_view = [
- get_config_key('cli-tools', 'samtools', 'path_to_samtools'),
+ program_paths['samtools'],
'view',
'-h',
'-b',
@@ -279,7 +192,7 @@ def align_to_genome(genome, small_rnas, cds, quiet=0, threads=4, small_rna_filet
]
unmapped_bedtools_bamtofastq_command = [
- get_config_key('cli-tools', 'samtools', 'path_to_samtools'),
+ program_paths['samtools'],
'fastq',
UNMAPPED_BAM,
'-o', RESULT_UNMAPPED_FASTQ,
@@ -287,7 +200,7 @@ def align_to_genome(genome, small_rnas, cds, quiet=0, threads=4, small_rna_filet
]
cds_unmapped_samtools_view = [
- get_config_key('cli-tools', 'samtools', 'path_to_samtools'),
+ program_paths['samtools'],
'view',
'-h',
'-b',
@@ -297,48 +210,67 @@ def align_to_genome(genome, small_rnas, cds, quiet=0, threads=4, small_rna_filet
]
cds_unmapped_bedtools_bamtofastq_command = [
- get_config_key('cli-tools', 'samtools', 'path_to_samtools'),
+ program_paths['samtools'],
'fastq',
CDS_UNMAPPED_BAM,
'-o', CDS_UNMAPPED_FASTQ,
'-0', CDS_UNMAPPED_FASTQ
]
- do_log(quiet, '====> Building BBMap Index')
- run(bbmap_build_index, capture_output=(quiet != 0))
+ print('====> Building bowtie2 Index...')
+ run(bowtie2_build_index, capture_output=not verbose)
if cds is not None:
- run(cds_build_index, capture_output=(quiet != 0))
+ run(cds_build_index, capture_output=not verbose)
+
+ print('====> Aligning reads to the genome...')
+ run(bowtie2_align_reads, capture_output=not verbose)
- do_log(quiet, '====> Aligning reads to the genome')
- run(bbmap_align_reads, capture_output=(quiet != 0))
+ run(samtools_view, capture_output=not verbose)
+ run(bedtools_bamtofastq_command, capture_output=not verbose)
+ run(unmapped_samtools_view, capture_output=not verbose)
+ run(unmapped_bedtools_bamtofastq_command, capture_output=not verbose)
if cds is not None:
- run(cds_align_reads, capture_output=(quiet != 0))
+ run(cds_align_reads, capture_output=not verbose)
- do_log(quiet, '====> Converting to FASTQ')
- run(samtools_view, capture_output=(quiet != 0))
- run(bedtools_bamtofastq_command, capture_output=(quiet != 0))
- run(unmapped_samtools_view, capture_output=(quiet != 0))
- run(unmapped_bedtools_bamtofastq_command, capture_output=(quiet != 0))
+ print('====> Converting to FASTQ...')
if cds is not None:
- run(cds_samtools_view, capture_output=(quiet != 0))
- run(cds_bedtools_bamtofastq_command, capture_output=(quiet != 0))
- run(cds_unmapped_samtools_view, capture_output=(quiet != 0))
- run(cds_unmapped_bedtools_bamtofastq_command, capture_output=(quiet != 0))
+ run(cds_samtools_view, capture_output=not verbose)
+ run(cds_bedtools_bamtofastq_command, capture_output=not verbose)
+ run(cds_unmapped_samtools_view, capture_output=not verbose)
+ run(cds_unmapped_bedtools_bamtofastq_command, capture_output=not verbose)
+
+ copy2(CDS_UNMAPPED_FASTQ, FINAL_UNMAPPED_FASTQ)
+ else:
+ copy2(RESULT_UNMAPPED_FASTQ, FINAL_UNMAPPED_FASTQ)
- make_fastqs_unique(RESULT_FASTQ, CDS_RESULT_FASTQ, FINAL_FASTQ)
- make_fastq_overlap_only(RESULT_UNMAPPED_FASTQ, CDS_UNMAPPED_FASTQ, FINAL_UNMAPPED_FASTQ)
+ SeqIO.write(into_seqrecord(merge_genome_cds_hits(RESULT_FASTQ, CDS_RESULT_FASTQ)), FINAL_FASTQ, 'fastq')
+ create_stats_table(small_rnas, output)
- create_stats_table(small_rnas, get_config_key('general', 'output_directory'), small_rna_filetype=small_rna_filetype)
+ if not keep_intermediate:
+ shutil.rmtree(INDEX_DIRECTORY)
+
+ os.remove(INTERMEDIATE_SAM)
+ os.remove(INTERMEDIATE_BAM)
+ os.remove(UNMAPPED_BAM)
+ os.remove(RESULT_FASTQ)
+ os.remove(RESULT_UNMAPPED_FASTQ)
+
+ if cds is not None:
+ os.remove(CDS_INTERMEDIATE_SAM)
+ os.remove(CDS_INTERMEDIATE_BAM)
+ os.remove(CDS_UNMAPPED_BAM)
+ os.remove(CDS_UNMAPPED_FASTQ)
+ os.remove(CDS_RESULT_FASTQ)
return FINAL_FASTQ
-def create_stats_table(smallrna, output_dir, small_rna_filetype='fastq'):
+def create_stats_table(smallrna, output_dir):
'''
Count sequences to create the statistics file
'''
input_reads = 0
- for read in SeqIO.parse(smallrna, small_rna_filetype):
+ for read in SeqIO.parse(smallrna, 'fastq'):
input_reads += 1
overall_mapped = 0
@@ -361,7 +293,7 @@ def create_stats_table(smallrna, output_dir, small_rna_filetype='fastq'):
cds_unmapped = input_reads - cds_mapped
genome_unmapped = input_reads - genome_mapped
- with open(os.path.join(output_dir, 'report.tsv'), 'w') as f:
+ with open(os.path.join(output_dir, 'alignment_report.tsv'), 'w') as f:
writer = csv.writer(f, delimiter='\t')
writer.writerow(['Total Reads', input_reads])
@@ -374,68 +306,81 @@ def create_stats_table(smallrna, output_dir, small_rna_filetype='fastq'):
writer.writerow(['Unmapped to Genome', genome_unmapped])
writer.writerow(['Percentage Unmapped to Genome', genome_unmapped / input_reads * 100])
writer.writerow(['Mapped to CDS', cds_mapped])
- writer.writerow(['Percentage Mapped to CDS', cds_mapped / input_reads * 100])
- writer.writerow(['Unmapped to CDS', cds_unmapped])
- writer.writerow(['Percentage Unmapped to CDS', cds_unmapped / input_reads * 100])
+ writer.writerow(['Percentage Mapped to CDS (but not to genome)', cds_mapped / input_reads * 100])
+ writer.writerow(['Unmapped to CDS or genome', cds_unmapped])
+ writer.writerow(['Percentage Unmapped to CDS or genome', cds_unmapped / input_reads * 100])
-def bin_rna_size(rna_file, min_length, max_length, quiet=0):
+def bin_rna_size(rna_file, min_length, max_length, output, verbose=False):
'''
- Bin the RNAs in the RNA file into new files by length
+ Bin the RNAs in the RNA file into new files by length and first base
'''
+ BIN_LENGTH_DIRECTORY = os.path.join(output, 'binned_length_rna')
+ if not os.path.exists(BIN_LENGTH_DIRECTORY):
+ os.makedirs(BIN_LENGTH_DIRECTORY)
- BINS_DIRECTORY = os.path.join(get_config_key('general', 'output_directory'), 'binned_rna')
- mkdir_if_not_exists(BINS_DIRECTORY)
+ BIN_BASE_DIRECTORY = os.path.join(output, 'binned_group_rna')
+ if not os.path.exists(BIN_BASE_DIRECTORY):
+ os.makedirs(BIN_BASE_DIRECTORY)
- do_log(quiet, '====> Sorting RNA into arrays by length')
- rnas = sorted(SeqIO.parse(rna_file, 'fastq'), key=lambda x: len(x))
+ print('====> Sorting RNA into arrays by length and first base...')
- table_path = os.path.join(get_config_key('general', 'output_directory'), 'rna_length_report.csv')
+ length_firstbase = defaultdict(lambda: defaultdict(lambda: []))
+ sequence_counts = defaultdict(lambda: 0)
- last_length = 0
- current_rnas = []
- start_base_count = {'A': 0, 'C': 0, 'G': 0, 'T': 0}
- with open(table_path, 'w') as table_file:
- writer = csv.writer(table_file)
- writer.writerow(['RNA Length', 'A', 'C', 'G', 'U', 'All Frequency'])
+ if os.path.exists(os.path.join(output, 'unmapped_sequences.fastq')):
+ unmapped = 0
+ for seq in SeqIO.parse(os.path.join(output, 'unmapped_sequences.fastq'), 'fastq'):
+ unmapped += 1
- for rna_seq in rnas:
- if len(rna_seq) != last_length:
- if len(current_rnas) > 0 and last_length >= min_length and last_length <= max_length:
- filename = os.path.join(BINS_DIRECTORY, 'length' + str(last_length) + '.fastq')
- with open(filename, 'a') as f:
- SeqIO.write(current_rnas, f, 'fastq')
+ sequence_counts['__not_aligned'] = unmapped
+ else:
+ sequence_counts['__not_aligned'] = 0
+
+ for seq in SeqIO.parse(rna_file, 'fastq'):
+ first_base = seq.seq[0].replace('T', 'U')
+ if len(seq) >= min_length and len(seq) <= max_length:
+ sequence_counts[str(seq.seq)] += 1
+ length_firstbase[len(seq)][first_base].append(seq)
- writer.writerow([
- last_length,
- start_base_count['A'], start_base_count['C'],
- start_base_count['G'], start_base_count['T'],
- len(current_rnas)
- ])
+ for length in length_firstbase:
+ length_rnas = [x for xs in length_firstbase[length].values() for x in xs]
+ SeqIO.write(length_rnas, os.path.join(BIN_LENGTH_DIRECTORY, f'length{length}.fastq'), 'fastq')
- current_rnas = []
- start_base_count = {'A': 0, 'C': 0, 'G': 0, 'T': 0}
+ for first_base in length_firstbase[length]:
+ SeqIO.write(length_firstbase[length][first_base], os.path.join(BIN_BASE_DIRECTORY, f'{length}{first_base}.fastq'), 'fastq')
- last_length = len(rna_seq)
+ table_path = os.path.join(output, 'rna_length_report.csv')
+ counts_path = os.path.join(output, 'counts.tab')
+ with open(counts_path, 'w') as f:
+ writer = csv.writer(f, delimiter='\t')
+ for sequence in sequence_counts:
+ writer.writerow([sequence, sequence_counts[sequence]])
- current_rnas.append(rna_seq)
- try:
- start_base_count[rna_seq[0]] += 1
- except KeyError:
- do_log(quiet, f'Warning : An RNA started with ambiguous nt {rna_seq[0]}')
+ with open(table_path, 'w') as table_file:
+ writer = csv.writer(table_file)
+ writer.writerow(['RNA Length', 'A', 'C', 'G', 'U', 'N', 'All Frequency'])
+
+ for length in range(min_length, max_length + 1):
+ length_counts = [
+ len(length_firstbase[length]['A']), len(length_firstbase[length]['C']),
+ len(length_firstbase[length]['G']), len(length_firstbase[length]['U']),
+ len(length_firstbase[length]['N'])
+ ]
+
+ writer.writerow([length] + length_counts + [sum(length_counts)])
return table_path
-def graph_length(path_to_table):
+def graph_length(path_to_table, output):
'''
Create the graph showing the length and starting base of all the RNA
'''
-
with open(path_to_table) as f:
reader = csv.reader(f)
next(reader)
labels = []
- bases = {'A': [], 'C': [], 'G': [], 'T': []}
+ bases = {'A': [], 'C': [], 'G': [], 'T': [], 'N': []}
for line in reader:
labels.append(str(line[0]))
@@ -444,12 +389,13 @@ def graph_length(path_to_table):
bases['C'].append(int(line[2]))
bases['G'].append(int(line[3]))
bases['T'].append(int(line[4]))
+ bases['N'].append(int(line[5]))
- total_rna = sum(bases['A']) + sum(bases['C']) + sum(bases['G']) + sum(bases['T'])
+ total_rna = sum(bases['A']) + sum(bases['C']) + sum(bases['G']) + sum(bases['T']) + sum(bases['N'])
for key in bases.keys():
bases[key] = np.array(bases[key]) / total_rna * 100
- with open(os.path.join(get_config_key('general', 'output_directory'), 'baseplot_data.csv'), 'w') as f:
+ with open(os.path.join(output, 'baseplot_data.csv'), 'w') as f:
writer = csv.writer(f)
writer.writerow(['RNA Length'] + labels)
@@ -457,6 +403,7 @@ def graph_length(path_to_table):
writer.writerow(['T'] + list(bases['T']))
writer.writerow(['C'] + list(bases['C']))
writer.writerow(['G'] + list(bases['G']))
+ writer.writerow(['N'] + list(bases['N']))
fig, ax = plt.subplots(figsize=(11.69, 8.27))
@@ -464,11 +411,13 @@ def graph_length(path_to_table):
ax.bar(labels, bases['C'], bottom=bases['A'], label='C', align='edge', width=0.8)
ax.bar(labels, bases['G'], bottom=bases['A'] + bases['C'], label='G', align='edge', width=0.8)
ax.bar(labels, bases['T'], bottom=bases['A'] + bases['C'] + bases['G'], label='U', align='edge', width=0.8)
+ ax.bar(labels, bases['N'], bottom=bases['A'] + bases['C'] + bases['G'] + bases['T'], label='N', align='edge', width=0.8)
+ ax.set_xticks(list(range(0, len(labels))))
ax.set_xticklabels(labels, fontsize=7)
ax.set_xlabel('Length of small RNA')
ax.set_ylabel('Percentage of small RNA')
plt.legend(title='First Nucliotide')
- plt.savefig(os.path.join(get_config_key('general', 'output_directory'), 'baseplot.png'))
\ No newline at end of file
+ plt.savefig(os.path.join(output, 'baseplot.png'))
\ No newline at end of file
diff --git a/hlsmallrna/kegg_pathway_names.py b/hlsmallrna/kegg_pathway_names.py
new file mode 100644
index 0000000..c4be98b
--- /dev/null
+++ b/hlsmallrna/kegg_pathway_names.py
@@ -0,0 +1,573 @@
+from collections import defaultdict
+
+KEGG_PATHWAY_NAMES = defaultdict(lambda: 'NONE', {
+ "map01100": "Metabolic pathways",
+ "map01110": "Biosynthesis of secondary metabolites",
+ "map01120": "Microbial metabolism in diverse environments",
+ "map01200": "Carbon metabolism",
+ "map01210": "2-Oxocarboxylic acid metabolism",
+ "map01212": "Fatty acid metabolism",
+ "map01230": "Biosynthesis of amino acids",
+ "map01232": "Nucleotide metabolism",
+ "map01250": "Biosynthesis of nucleotide sugars",
+ "map01240": "Biosynthesis of cofactors",
+ "map01220": "Degradation of aromatic compounds",
+ "map00010": "Glycolysis / Gluconeogenesis",
+ "map00020": "Citrate cycle (TCA cycle)",
+ "map00030": "Pentose phosphate pathway",
+ "map00040": "Pentose and glucuronate interconversions",
+ "map00051": "Fructose and mannose metabolism",
+ "map00052": "Galactose metabolism",
+ "map00053": "Ascorbate and aldarate metabolism",
+ "map00500": "Starch and sucrose metabolism",
+ "map00520": "Amino sugar and nucleotide sugar metabolism",
+ "map00620": "Pyruvate metabolism",
+ "map00630": "Glyoxylate and dicarboxylate metabolism",
+ "map00640": "Propanoate metabolism",
+ "map00650": "Butanoate metabolism",
+ "map00660": "C5-Branched dibasic acid metabolism",
+ "map00562": "Inositol phosphate metabolism",
+ "map00190": "Oxidative phosphorylation",
+ "map00195": "Photosynthesis",
+ "map00196": "Photosynthesis - antenna proteins",
+ "map00710": "Carbon fixation in photosynthetic organisms",
+ "map00720": "Carbon fixation pathways in prokaryotes",
+ "map00680": "Methane metabolism",
+ "map00910": "Nitrogen metabolism",
+ "map00920": "Sulfur metabolism",
+ "map00061": "Fatty acid biosynthesis",
+ "map00062": "Fatty acid elongation",
+ "map00071": "Fatty acid degradation",
+ "map00073": "Cutin, suberine and wax biosynthesis",
+ "map00100": "Steroid biosynthesis",
+ "map00120": "Primary bile acid biosynthesis",
+ "map00121": "Secondary bile acid biosynthesis",
+ "map00140": "Steroid hormone biosynthesis",
+ "map00561": "Glycerolipid metabolism",
+ "map00564": "Glycerophospholipid metabolism",
+ "map00565": "Ether lipid metabolism",
+ "map00600": "Sphingolipid metabolism",
+ "map00590": "Arachidonic acid metabolism",
+ "map00591": "Linoleic acid metabolism",
+ "map00592": "alpha-Linolenic acid metabolism",
+ "map01040": "Biosynthesis of unsaturated fatty acids",
+ "map00230": "Purine metabolism",
+ "map00240": "Pyrimidine metabolism",
+ "map00250": "Alanine, aspartate and glutamate metabolism",
+ "map00260": "Glycine, serine and threonine metabolism",
+ "map00270": "Cysteine and methionine metabolism",
+ "map00280": "Valine, leucine and isoleucine degradation",
+ "map00290": "Valine, leucine and isoleucine biosynthesis",
+ "map00300": "Lysine biosynthesis",
+ "map00310": "Lysine degradation",
+ "map00220": "Arginine biosynthesis",
+ "map00330": "Arginine and proline metabolism",
+ "map00340": "Histidine metabolism",
+ "map00350": "Tyrosine metabolism",
+ "map00360": "Phenylalanine metabolism",
+ "map00380": "Tryptophan metabolism",
+ "map00400": "Phenylalanine, tyrosine and tryptophan biosynthesis",
+ "map00410": "beta-Alanine metabolism",
+ "map00430": "Taurine and hypotaurine metabolism",
+ "map00440": "Phosphonate and phosphinate metabolism",
+ "map00450": "Selenocompound metabolism",
+ "map00460": "Cyanoamino acid metabolism",
+ "map00470": "D-Amino acid metabolism",
+ "map00480": "Glutathione metabolism",
+ "map00510": "N-Glycan biosynthesis",
+ "map00513": "Various types of N-glycan biosynthesis",
+ "map00512": "Mucin type O-glycan biosynthesis",
+ "map00515": "Mannose type O-glycan biosynthesis",
+ "map00514": "Other types of O-glycan biosynthesis",
+ "map00532": "Glycosaminoglycan biosynthesis - chondroitin sulfate / dermatan sulfate",
+ "map00534": "Glycosaminoglycan biosynthesis - heparan sulfate / heparin",
+ "map00533": "Glycosaminoglycan biosynthesis - keratan sulfate",
+ "map00531": "Glycosaminoglycan degradation",
+ "map00563": "Glycosylphosphatidylinositol (GPI)-anchor biosynthesis",
+ "map00601": "Glycosphingolipid biosynthesis - lacto and neolacto series",
+ "map00603": "Glycosphingolipid biosynthesis - globo and isoglobo series",
+ "map00604": "Glycosphingolipid biosynthesis - ganglio series",
+ "map00511": "Other glycan degradation",
+ "map00540": "Lipopolysaccharide biosynthesis",
+ "map00542": "O-Antigen repeat unit biosynthesis",
+ "map00541": "O-Antigen nucleotide sugar biosynthesis",
+ "map00550": "Peptidoglycan biosynthesis",
+ "map00552": "Teichoic acid biosynthesis",
+ "map00571": "Lipoarabinomannan (LAM) biosynthesis",
+ "map00572": "Arabinogalactan biosynthesis - Mycobacterium",
+ "map00543": "Exopolysaccharide biosynthesis",
+ "map00730": "Thiamine metabolism",
+ "map00740": "Riboflavin metabolism",
+ "map00750": "Vitamin B6 metabolism",
+ "map00760": "Nicotinate and nicotinamide metabolism",
+ "map00770": "Pantothenate and CoA biosynthesis",
+ "map00780": "Biotin metabolism",
+ "map00785": "Lipoic acid metabolism",
+ "map00790": "Folate biosynthesis",
+ "map00670": "One carbon pool by folate",
+ "map00830": "Retinol metabolism",
+ "map00860": "Porphyrin metabolism",
+ "map00130": "Ubiquinone and other terpenoid-quinone biosynthesis",
+ "map00900": "Terpenoid backbone biosynthesis",
+ "map00902": "Monoterpenoid biosynthesis",
+ "map00909": "Sesquiterpenoid and triterpenoid biosynthesis",
+ "map00904": "Diterpenoid biosynthesis",
+ "map00906": "Carotenoid biosynthesis",
+ "map00905": "Brassinosteroid biosynthesis",
+ "map00981": "Insect hormone biosynthesis",
+ "map00908": "Zeatin biosynthesis",
+ "map00903": "Limonene degradation",
+ "map00907": "Pinene, camphor and geraniol degradation",
+ "map01052": "Type I polyketide structures",
+ "map00522": "Biosynthesis of 12-, 14- and 16-membered macrolides",
+ "map01051": "Biosynthesis of ansamycins",
+ "map01059": "Biosynthesis of enediyne antibiotics",
+ "map01056": "Biosynthesis of type II polyketide backbone",
+ "map01057": "Biosynthesis of type II polyketide products",
+ "map00253": "Tetracycline biosynthesis",
+ "map00523": "Polyketide sugar unit biosynthesis",
+ "map01054": "Nonribosomal peptide structures",
+ "map01053": "Biosynthesis of siderophore group nonribosomal peptides",
+ "map01055": "Biosynthesis of vancomycin group antibiotics",
+ "map00940": "Phenylpropanoid biosynthesis",
+ "map00945": "Stilbenoid, diarylheptanoid and gingerol biosynthesis",
+ "map00941": "Flavonoid biosynthesis",
+ "map00944": "Flavone and flavonol biosynthesis",
+ "map00942": "Anthocyanin biosynthesis",
+ "map00943": "Isoflavonoid biosynthesis",
+ "map00946": "Degradation of flavonoids",
+ "map00901": "Indole alkaloid biosynthesis",
+ "map00403": "Indole diterpene alkaloid biosynthesis",
+ "map00950": "Isoquinoline alkaloid biosynthesis",
+ "map00960": "Tropane, piperidine and pyridine alkaloid biosynthesis",
+ "map00996": "Biosynthesis of various alkaloids",
+ "map00232": "Caffeine metabolism",
+ "map00965": "Betalain biosynthesis",
+ "map00966": "Glucosinolate biosynthesis",
+ "map00402": "Benzoxazinoid biosynthesis",
+ "map00311": "Penicillin and cephalosporin biosynthesis",
+ "map00332": "Carbapenem biosynthesis",
+ "map00261": "Monobactam biosynthesis",
+ "map00331": "Clavulanic acid biosynthesis",
+ "map00521": "Streptomycin biosynthesis",
+ "map00524": "Neomycin, kanamycin and gentamicin biosynthesis",
+ "map00525": "Acarbose and validamycin biosynthesis",
+ "map00401": "Novobiocin biosynthesis",
+ "map00404": "Staurosporine biosynthesis",
+ "map00405": "Phenazine biosynthesis",
+ "map00333": "Prodigiosin biosynthesis",
+ "map00254": "Aflatoxin biosynthesis",
+ "map00998": "Biosynthesis of various antibiotics",
+ "map00999": "Biosynthesis of various plant secondary metabolites",
+ "map00997": "Biosynthesis of various other secondary metabolites",
+ "map00362": "Benzoate degradation",
+ "map00627": "Aminobenzoate degradation",
+ "map00364": "Fluorobenzoate degradation",
+ "map00625": "Chloroalkane and chloroalkene degradation",
+ "map00361": "Chlorocyclohexane and chlorobenzene degradation",
+ "map00623": "Toluene degradation",
+ "map00622": "Xylene degradation",
+ "map00633": "Nitrotoluene degradation",
+ "map00642": "Ethylbenzene degradation",
+ "map00643": "Styrene degradation",
+ "map00791": "Atrazine degradation",
+ "map00930": "Caprolactam degradation",
+ "map00363": "Bisphenol degradation",
+ "map00621": "Dioxin degradation",
+ "map00626": "Naphthalene degradation",
+ "map00624": "Polycyclic aromatic hydrocarbon degradation",
+ "map00365": "Furfural degradation",
+ "map00984": "Steroid degradation",
+ "map00980": "Metabolism of xenobiotics by cytochrome P450",
+ "map00982": "Drug metabolism - cytochrome P450",
+ "map00983": "Drug metabolism - other enzymes",
+ "map01010": "Overview of biosynthetic pathways",
+ "map01060": "Biosynthesis of plant secondary metabolites",
+ "map01061": "Biosynthesis of phenylpropanoids",
+ "map01062": "Biosynthesis of terpenoids and steroids",
+ "map01063": "Biosynthesis of alkaloids derived from shikimate pathway",
+ "map01064": "Biosynthesis of alkaloids derived from ornithine, lysine and nicotinic acid",
+ "map01065": "Biosynthesis of alkaloids derived from histidine and purine",
+ "map01066": "Biosynthesis of alkaloids derived from terpenoid and polyketide",
+ "map01070": "Biosynthesis of plant hormones",
+ "map03020": "RNA polymerase",
+ "map03022": "Basal transcription factors",
+ "map03040": "Spliceosome",
+ "map03010": "Ribosome",
+ "map00970": "Aminoacyl-tRNA biosynthesis",
+ "map03013": "Nucleocytoplasmic transport",
+ "map03015": "mRNA surveillance pathway",
+ "map03008": "Ribosome biogenesis in eukaryotes",
+ "map03060": "Protein export",
+ "map04141": "Protein processing in endoplasmic reticulum",
+ "map04130": "SNARE interactions in vesicular transport",
+ "map04120": "Ubiquitin mediated proteolysis",
+ "map04122": "Sulfur relay system",
+ "map03050": "Proteasome",
+ "map03018": "RNA degradation",
+ "map03030": "DNA replication",
+ "map03410": "Base excision repair",
+ "map03420": "Nucleotide excision repair",
+ "map03430": "Mismatch repair",
+ "map03440": "Homologous recombination",
+ "map03450": "Non-homologous end-joining",
+ "map03460": "Fanconi anemia pathway",
+ "map03082": "ATP-dependent chromatin remodeling",
+ "map03083": "Polycomb repressive complex",
+ "map03230": "Viral genome structure",
+ "map03240": "Viral replication",
+ "map03250": "Viral life cycle - HIV-1",
+ "map03260": "Virion - Human immunodeficiency virus",
+ "map03271": "Virion - Rotavirus",
+ "map03262": "Virion - Coronavirus",
+ "map03261": "Virion - Influenza virus",
+ "map03263": "Virion - Enterovirus",
+ "map03264": "Virion - Flavivirus",
+ "map03265": "Virion - Ebolavirus and Lyssavirus",
+ "map03266": "Virion - Herpesvirus",
+ "map03268": "Virion - Orthopoxvirus",
+ "map03267": "Virion - Adenovirus",
+ "map03259": "Virion - Bacteriophages T4 and T7",
+ "map03258": "Virion - Bacteriophage lambda",
+ "map02010": "ABC transporters",
+ "map02060": "Phosphotransferase system (PTS)",
+ "map03070": "Bacterial secretion system",
+ "map02020": "Two-component system",
+ "map04010": "MAPK signaling pathway",
+ "map04013": "MAPK signaling pathway - fly",
+ "map04016": "MAPK signaling pathway - plant",
+ "map04011": "MAPK signaling pathway - yeast",
+ "map04012": "ErbB signaling pathway",
+ "map04014": "Ras signaling pathway",
+ "map04015": "Rap1 signaling pathway",
+ "map04310": "Wnt signaling pathway",
+ "map04330": "Notch signaling pathway",
+ "map04340": "Hedgehog signaling pathway",
+ "map04341": "Hedgehog signaling pathway - fly",
+ "map04350": "TGF-beta signaling pathway",
+ "map04390": "Hippo signaling pathway",
+ "map04391": "Hippo signaling pathway - fly",
+ "map04392": "Hippo signaling pathway - multiple species",
+ "map04370": "VEGF signaling pathway",
+ "map04371": "Apelin signaling pathway",
+ "map04630": "JAK-STAT signaling pathway",
+ "map04064": "NF-kappa B signaling pathway",
+ "map04668": "TNF signaling pathway",
+ "map04066": "HIF-1 signaling pathway",
+ "map04068": "FoxO signaling pathway",
+ "map04020": "Calcium signaling pathway",
+ "map04070": "Phosphatidylinositol signaling system",
+ "map04072": "Phospholipase D signaling pathway",
+ "map04071": "Sphingolipid signaling pathway",
+ "map04024": "cAMP signaling pathway",
+ "map04022": "cGMP-PKG signaling pathway",
+ "map04151": "PI3K-Akt signaling pathway",
+ "map04152": "AMPK signaling pathway",
+ "map04150": "mTOR signaling pathway",
+ "map04075": "Plant hormone signal transduction",
+ "map04080": "Neuroactive ligand-receptor interaction",
+ "map04060": "Cytokine-cytokine receptor interaction",
+ "map04061": "Viral protein interaction with cytokine and cytokine receptor",
+ "map04512": "ECM-receptor interaction",
+ "map04514": "Cell adhesion molecules",
+ "map04144": "Endocytosis",
+ "map04145": "Phagosome",
+ "map04142": "Lysosome",
+ "map04146": "Peroxisome",
+ "map04140": "Autophagy - animal",
+ "map04138": "Autophagy - yeast",
+ "map04136": "Autophagy - other",
+ "map04137": "Mitophagy - animal",
+ "map04139": "Mitophagy - yeast",
+ "map04148": "Efferocytosis",
+ "map04110": "Cell cycle",
+ "map04111": "Cell cycle - yeast",
+ "map04112": "Cell cycle - Caulobacter",
+ "map04113": "Meiosis - yeast",
+ "map04114": "Oocyte meiosis",
+ "map04210": "Apoptosis",
+ "map04214": "Apoptosis - fly",
+ "map04215": "Apoptosis - multiple species",
+ "map04216": "Ferroptosis",
+ "map04217": "Necroptosis",
+ "map04115": "p53 signaling pathway",
+ "map04218": "Cellular senescence",
+ "map04510": "Focal adhesion",
+ "map04520": "Adherens junction",
+ "map04530": "Tight junction",
+ "map04540": "Gap junction",
+ "map04550": "Signaling pathways regulating pluripotency of stem cells",
+ "map02024": "Quorum sensing",
+ "map05111": "Biofilm formation - Vibrio cholerae",
+ "map02025": "Biofilm formation - Pseudomonas aeruginosa",
+ "map02026": "Biofilm formation - Escherichia coli",
+ "map02030": "Bacterial chemotaxis",
+ "map02040": "Flagellar assembly",
+ "map04814": "Motor proteins",
+ "map04810": "Regulation of actin cytoskeleton",
+ "map04640": "Hematopoietic cell lineage",
+ "map04610": "Complement and coagulation cascades",
+ "map04611": "Platelet activation",
+ "map04613": "Neutrophil extracellular trap formation",
+ "map04620": "Toll-like receptor signaling pathway",
+ "map04624": "Toll and Imd signaling pathway",
+ "map04621": "NOD-like receptor signaling pathway",
+ "map04622": "RIG-I-like receptor signaling pathway",
+ "map04623": "Cytosolic DNA-sensing pathway",
+ "map04625": "C-type lectin receptor signaling pathway",
+ "map04650": "Natural killer cell mediated cytotoxicity",
+ "map04612": "Antigen processing and presentation",
+ "map04660": "T cell receptor signaling pathway",
+ "map04658": "Th1 and Th2 cell differentiation",
+ "map04659": "Th17 cell differentiation",
+ "map04657": "IL-17 signaling pathway",
+ "map04662": "B cell receptor signaling pathway",
+ "map04664": "Fc epsilon RI signaling pathway",
+ "map04666": "Fc gamma R-mediated phagocytosis",
+ "map04670": "Leukocyte transendothelial migration",
+ "map04672": "Intestinal immune network for IgA production",
+ "map04062": "Chemokine signaling pathway",
+ "map04911": "Insulin secretion",
+ "map04910": "Insulin signaling pathway",
+ "map04922": "Glucagon signaling pathway",
+ "map04923": "Regulation of lipolysis in adipocytes",
+ "map04920": "Adipocytokine signaling pathway",
+ "map03320": "PPAR signaling pathway",
+ "map04929": "GnRH secretion",
+ "map04912": "GnRH signaling pathway",
+ "map04913": "Ovarian steroidogenesis",
+ "map04915": "Estrogen signaling pathway",
+ "map04914": "Progesterone-mediated oocyte maturation",
+ "map04917": "Prolactin signaling pathway",
+ "map04921": "Oxytocin signaling pathway",
+ "map04926": "Relaxin signaling pathway",
+ "map04935": "Growth hormone synthesis, secretion and action",
+ "map04918": "Thyroid hormone synthesis",
+ "map04919": "Thyroid hormone signaling pathway",
+ "map04928": "Parathyroid hormone synthesis, secretion and action",
+ "map04916": "Melanogenesis",
+ "map04924": "Renin secretion",
+ "map04614": "Renin-angiotensin system",
+ "map04925": "Aldosterone synthesis and secretion",
+ "map04927": "Cortisol synthesis and secretion",
+ "map04260": "Cardiac muscle contraction",
+ "map04261": "Adrenergic signaling in cardiomyocytes",
+ "map04270": "Vascular smooth muscle contraction",
+ "map04970": "Salivary secretion",
+ "map04971": "Gastric acid secretion",
+ "map04972": "Pancreatic secretion",
+ "map04976": "Bile secretion",
+ "map04973": "Carbohydrate digestion and absorption",
+ "map04974": "Protein digestion and absorption",
+ "map04975": "Fat digestion and absorption",
+ "map04979": "Cholesterol metabolism",
+ "map04977": "Vitamin digestion and absorption",
+ "map04978": "Mineral absorption",
+ "map04962": "Vasopressin-regulated water reabsorption",
+ "map04960": "Aldosterone-regulated sodium reabsorption",
+ "map04961": "Endocrine and other factor-regulated calcium reabsorption",
+ "map04964": "Proximal tubule bicarbonate reclamation",
+ "map04966": "Collecting duct acid secretion",
+ "map04724": "Glutamatergic synapse",
+ "map04727": "GABAergic synapse",
+ "map04725": "Cholinergic synapse",
+ "map04728": "Dopaminergic synapse",
+ "map04726": "Serotonergic synapse",
+ "map04720": "Long-term potentiation",
+ "map04730": "Long-term depression",
+ "map04723": "Retrograde endocannabinoid signaling",
+ "map04721": "Synaptic vesicle cycle",
+ "map04722": "Neurotrophin signaling pathway",
+ "map04744": "Phototransduction",
+ "map04745": "Phototransduction - fly",
+ "map04740": "Olfactory transduction",
+ "map04742": "Taste transduction",
+ "map04750": "Inflammatory mediator regulation of TRP channels",
+ "map04320": "Dorso-ventral axis formation",
+ "map04360": "Axon guidance",
+ "map04361": "Axon regeneration",
+ "map04380": "Osteoclast differentiation",
+ "map04211": "Longevity regulating pathway",
+ "map04212": "Longevity regulating pathway - worm",
+ "map04213": "Longevity regulating pathway - multiple species",
+ "map04710": "Circadian rhythm",
+ "map04713": "Circadian entrainment",
+ "map04711": "Circadian rhythm - fly",
+ "map04712": "Circadian rhythm - plant",
+ "map04714": "Thermogenesis",
+ "map04626": "Plant-pathogen interaction",
+ "map05200": "Pathways in cancer",
+ "map05202": "Transcriptional misregulation in cancer",
+ "map05206": "MicroRNAs in cancer",
+ "map05205": "Proteoglycans in cancer",
+ "map05204": "Chemical carcinogenesis - DNA adducts",
+ "map05207": "Chemical carcinogenesis - receptor activation",
+ "map05208": "Chemical carcinogenesis - reactive oxygen species",
+ "map05203": "Viral carcinogenesis",
+ "map05230": "Central carbon metabolism in cancer",
+ "map05231": "Choline metabolism in cancer",
+ "map05235": "PD-L1 expression and PD-1 checkpoint pathway in cancer",
+ "map05210": "Colorectal cancer",
+ "map05212": "Pancreatic cancer",
+ "map05225": "Hepatocellular carcinoma",
+ "map05226": "Gastric cancer",
+ "map05214": "Glioma",
+ "map05216": "Thyroid cancer",
+ "map05221": "Acute myeloid leukemia",
+ "map05220": "Chronic myeloid leukemia",
+ "map05217": "Basal cell carcinoma",
+ "map05218": "Melanoma",
+ "map05211": "Renal cell carcinoma",
+ "map05219": "Bladder cancer",
+ "map05215": "Prostate cancer",
+ "map05213": "Endometrial cancer",
+ "map05224": "Breast cancer",
+ "map05222": "Small cell lung cancer",
+ "map05223": "Non-small cell lung cancer",
+ "map05166": "Human T-cell leukemia virus 1 infection",
+ "map05170": "Human immunodeficiency virus 1 infection",
+ "map05161": "Hepatitis B",
+ "map05160": "Hepatitis C",
+ "map05171": "Coronavirus disease - COVID-19",
+ "map05164": "Influenza A",
+ "map05162": "Measles",
+ "map05168": "Herpes simplex virus 1 infection",
+ "map05163": "Human cytomegalovirus infection",
+ "map05167": "Kaposi sarcoma-associated herpesvirus infection",
+ "map05169": "Epstein-Barr virus infection",
+ "map05165": "Human papillomavirus infection",
+ "map05110": "Vibrio cholerae infection",
+ "map05120": "Epithelial cell signaling in Helicobacter pylori infection",
+ "map05130": "Pathogenic Escherichia coli infection",
+ "map05132": "Salmonella infection",
+ "map05131": "Shigellosis",
+ "map05135": "Yersinia infection",
+ "map05133": "Pertussis",
+ "map05134": "Legionellosis",
+ "map05150": "Staphylococcus aureus infection",
+ "map05152": "Tuberculosis",
+ "map05100": "Bacterial invasion of epithelial cells",
+ "map05146": "Amoebiasis",
+ "map05144": "Malaria",
+ "map05145": "Toxoplasmosis",
+ "map05140": "Leishmaniasis",
+ "map05142": "Chagas disease",
+ "map05143": "African trypanosomiasis",
+ "map05310": "Asthma",
+ "map05322": "Systemic lupus erythematosus",
+ "map05323": "Rheumatoid arthritis",
+ "map05320": "Autoimmune thyroid disease",
+ "map05321": "Inflammatory bowel disease",
+ "map05330": "Allograft rejection",
+ "map05332": "Graft-versus-host disease",
+ "map05340": "Primary immunodeficiency",
+ "map05010": "Alzheimer disease",
+ "map05012": "Parkinson disease",
+ "map05014": "Amyotrophic lateral sclerosis",
+ "map05016": "Huntington disease",
+ "map05017": "Spinocerebellar ataxia",
+ "map05020": "Prion disease",
+ "map05022": "Pathways of neurodegeneration - multiple diseases",
+ "map05030": "Cocaine addiction",
+ "map05031": "Amphetamine addiction",
+ "map05032": "Morphine addiction",
+ "map05033": "Nicotine addiction",
+ "map05034": "Alcoholism",
+ "map05417": "Lipid and atherosclerosis",
+ "map05418": "Fluid shear stress and atherosclerosis",
+ "map05410": "Hypertrophic cardiomyopathy",
+ "map05412": "Arrhythmogenic right ventricular cardiomyopathy",
+ "map05414": "Dilated cardiomyopathy",
+ "map05415": "Diabetic cardiomyopathy",
+ "map05416": "Viral myocarditis",
+ "map04930": "Type II diabetes mellitus",
+ "map04940": "Type I diabetes mellitus",
+ "map04950": "Maturity onset diabetes of the young",
+ "map04936": "Alcoholic liver disease",
+ "map04932": "Non-alcoholic fatty liver disease",
+ "map04931": "Insulin resistance",
+ "map04933": "AGE-RAGE signaling pathway in diabetic complications",
+ "map04934": "Cushing syndrome",
+ "map01501": "beta-Lactam resistance",
+ "map01502": "Vancomycin resistance",
+ "map01503": "Cationic antimicrobial peptide (CAMP) resistance",
+ "map01521": "EGFR tyrosine kinase inhibitor resistance",
+ "map01524": "Platinum drug resistance",
+ "map01523": "Antifolate resistance",
+ "map01522": "Endocrine resistance",
+ "map07011": "Penicillins",
+ "map07012": "Cephalosporins - parenteral agents",
+ "map07013": "Cephalosporins - oral agents",
+ "map07021": "Aminoglycosides",
+ "map07019": "Tetracyclines",
+ "map07020": "Macrolides and ketolides",
+ "map07014": "Quinolones",
+ "map07023": "Rifamycins",
+ "map07026": "Antifungal agents",
+ "map07044": "Antiviral agents",
+ "map07053": "Anti-HIV agents",
+ "map07040": "Antineoplastics - alkylating agents",
+ "map07041": "Antineoplastics - antimetabolic agents",
+ "map07042": "Antineoplastics - agents from natural products",
+ "map07043": "Antineoplastics - hormones",
+ "map07045": "Antineoplastics - protein kinase inhibitors",
+ "map07032": "Hypnotics",
+ "map07030": "Anxiolytics",
+ "map07033": "Anticonvulsants",
+ "map07015": "Local analgesics",
+ "map07039": "Opioid analgesics",
+ "map07028": "Antipsychotics",
+ "map07029": "Antipsychotics - phenothiazines",
+ "map07031": "Antipsychotics - butyrophenones",
+ "map07027": "Antidepressants",
+ "map07056": "Agents for Alzheimer-type dementia",
+ "map07057": "Antiparkinsonian agents",
+ "map07055": "Sulfonamide derivatives - overview",
+ "map07016": "Sulfonamide derivatives - sulfa drugs",
+ "map07017": "Sulfonamide derivatives - diuretics",
+ "map07018": "Sulfonamide derivatives - hypoglycemic agents",
+ "map07037": "Antiarrhythmic drugs",
+ "map07038": "Antiulcer drugs",
+ "map07046": "Immunosuppressive agents",
+ "map07047": "Osteoporosis drugs",
+ "map07048": "Antimigraines",
+ "map07049": "Antithrombosis agents",
+ "map07050": "Antirheumatics - DMARDs and biological agents",
+ "map07051": "Antidiabetics",
+ "map07052": "Antidyslipidemic agents",
+ "map07054": "Antiglaucoma agents",
+ "map07220": "Cholinergic and anticholinergic drugs",
+ "map07215": "alpha-Adrenergic receptor agonists/antagonists",
+ "map07214": "beta-Adrenergic receptor agonists/antagonists",
+ "map07213": "Dopamine receptor agonists/antagonists",
+ "map07212": "Histamine H1 receptor antagonists",
+ "map07227": "Histamine H2/H3 receptor agonists/antagonists",
+ "map07211": "Serotonin receptor agonists/antagonists",
+ "map07228": "Eicosanoid receptor agonists/antagonists",
+ "map07224": "Opioid receptor agonists/antagonists",
+ "map07229": "Angiotensin receptor and endothelin receptor antagonists",
+ "map07225": "Glucocorticoid and mineralocorticoid receptor agonists/antagonists",
+ "map07226": "Progesterone, androgen and estrogen receptor agonists/antagonists",
+ "map07223": "Retinoic acid receptor (RAR) and retinoid X receptor (RXR) agonists/antagonists",
+ "map07222": "Peroxisome proliferator-activated receptor (PPAR) agonists",
+ "map07221": "Nicotinic cholinergic receptor antagonists",
+ "map07230": "GABA-A receptor agonists/antagonists",
+ "map07036": "Calcium channel blocking drugs",
+ "map07231": "Sodium channel blocking drugs",
+ "map07232": "Potassium channel blocking and opening drugs",
+ "map07235": "N-Methyl-D-aspartic acid (NMDA) receptor antagonists",
+ "map07233": "Ion transporter inhibitors",
+ "map07234": "Neurotransmitter transporter inhibitors",
+ "map07216": "Catecholamine transferase inhibitors",
+ "map07219": "Cyclooxygenase inhibitors",
+ "map07024": "HMG-CoA reductase inhibitors",
+ "map07217": "Renin-angiotensin system inhibitors",
+ "map07218": "HIV protease inhibitors",
+ "map07025": "Quinolines",
+ "map07034": "Eicosanoids",
+ "map07035": "Prostaglandins",
+ "map07110": "Benzoic acid family",
+ "map07112": "1,2-Diphenyl substitution family",
+ "map07114": "Naphthalene family",
+ "map07117": "Benzodiazepine family"
+})
\ No newline at end of file
diff --git a/bin/label_for_unitas b/hlsmallrna/label_for_unitas.py
similarity index 71%
rename from bin/label_for_unitas
rename to hlsmallrna/label_for_unitas.py
index 71f029b..269a913 100755
--- a/bin/label_for_unitas
+++ b/hlsmallrna/label_for_unitas.py
@@ -1,5 +1,4 @@
-#!/usr/bin/env python3
-# Copyright 2022 Vicky Hunt Lab Members
+# Copyright 2022 - 2025 Vicky Hunt Lab Members
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
@@ -18,12 +17,25 @@
from Bio import SeqIO
def prepend_label_to_id(label, seqs):
+ '''
+ Add the specified label to the start of the sequence ID
+ '''
for seq in seqs:
seq.id = label + '|' + seq.id
yield seq
-if __name__ == '__main__':
+def label_file_for_unitas(input, output, label):
+ '''
+ Prepend a specified label to all sequences in a FASTA file
+ '''
+ SeqIO.write(
+ prepend_label_to_id(label, SeqIO.parse(input, 'fasta')),
+ output,
+ 'fasta'
+ )
+
+def label_for_unitas_cli():
parser = ArgumentParser(description='Program to automatically prepend unitas labels to all the sequences in a FASTA file')
parser.add_argument('-o', '--output', help='File to output to, defaults to labelled.fasta', default='labeled.fasta')
@@ -33,7 +45,4 @@ def prepend_label_to_id(label, seqs):
args = parser.parse_args()
- data = SeqIO.parse(args.file_path, 'fasta')
- labelled_data = prepend_label_to_id(args.label, data)
-
- SeqIO.write(labelled_data, args.output, 'fasta')
\ No newline at end of file
+ label_file_for_unitas(args.file_path, args.output, args.label)
\ No newline at end of file
diff --git a/hlsmallrna/ss_overlap.py b/hlsmallrna/ss_overlap.py
deleted file mode 100644
index 47c0a45..0000000
--- a/hlsmallrna/ss_overlap.py
+++ /dev/null
@@ -1,45 +0,0 @@
-# Copyright 2022 Vicky Hunt Lab Members
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-# http://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software
-# distributed under the License is distributed on an "AS IS" BASIS,
-# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-# See the License for the specific language governing permissions and
-# limitations under the License.
-import os
-
-from subprocess import run
-
-from .config import get_config_key, mkdir_if_not_exists, do_log
-
-def samestrand_overlap(genome_file, rna_file_1, rna_file_2, longest_rna, quiet=0):
- '''
- Run the samestand overlap script
- '''
-
- do_log(quiet, '==> Running same strand overlap script')
-
- CWD = os.getcwd()
- OUTPUT_DIR = os.path.join(CWD, get_config_key('general', 'output_directory'), 'samestrand_overlap')
-
- mkdir_if_not_exists(OUTPUT_DIR)
- os.chdir(OUTPUT_DIR)
-
- command = [
- 'overlap_ss.sh',
- os.path.join(CWD, genome_file),
- os.path.join(CWD, rna_file_1),
- os.path.join(CWD, rna_file_2),
- str(longest_rna)
- ]
-
- run(command, capture_output=(quiet != 0))
-
- os.chdir(CWD)
-
- do_log(quiet, '==> Finished same strand overlap script')
\ No newline at end of file
diff --git a/hlsmallrna/targetid.py b/hlsmallrna/targetid.py
index 9120bf4..d8a8495 100644
--- a/hlsmallrna/targetid.py
+++ b/hlsmallrna/targetid.py
@@ -1,4 +1,4 @@
-# Copyright 2022 Vicky Hunt Lab Members
+# Copyright 2022 - 2025 Vicky Hunt Lab Members
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
@@ -16,12 +16,54 @@
import csv
from subprocess import run
+from collections import defaultdict, Counter
import pysam
from Bio import SeqIO
+from Bio.Seq import Seq
+from scipy.stats import fisher_exact
-from .config import get_config_key, mkdir_if_not_exists, do_log
+from .kegg_pathway_names import KEGG_PATHWAY_NAMES
+
+GO_ENRICHMENT_SCRIPT = '''
+interesting_genes <- c({{INTERSTING_GENE_LIST}})
+
+library(topGO)
+
+all_data <- readMappings("{{PATH_TO_GO_TXT}}")
+universe <- names(all_data)
+duplicated_list <- factor(as.integer(universe %in% interesting_genes))
+names(duplicated_list) <- universe
+
+go_ids <- c()
+go_signif <- c()
+go_names <- c()
+
+for (ont in c('BP', 'CC', 'MF')) {
+ myGOdata <- new("topGOdata", description="My project", ontology=ont, allGenes=duplicated_list, annot = annFUN.gene2GO, gene2GO = all_data)
+ significant <- sigGenes(myGOdata)
+ resultFisher <- runTest(myGOdata, algorithm="classic", statistic="fisher")
+ mysummary <- summary(attributes(resultFisher)$score <= 0.01)
+
+ if (length(mysummary) > 2) {
+ numsignif <- as.integer(mysummary[[3]])
+
+ allRes <- GenTable(myGOdata, classicFisher = resultFisher, orderBy = "resultFisher", ranksOf = "classicFisher", topNodes = numsignif)
+ go_ids <- c(go_ids, allRes$GO.ID)
+ go_names <- c(go_names, allRes$Term)
+ go_signif <- c(go_signif, allRes$classicFisher)
+ }
+}
+
+lines <- c()
+for (i in 1:length(go_ids)) {
+ lines <- append(lines, c(paste(go_ids[i], go_names[i], go_signif[i], sep = "\t")))
+}
+fileConn <- file("{{PATH_TO_OUTPUT_TSV}}")
+writeLines(lines, fileConn)
+close(fileConn)
+'''
def check_ids_unique(path_to_file, format):
'''
@@ -36,45 +78,62 @@ def check_ids_unique(path_to_file, format):
return True
-def revcomp_input_file(smallRNA, quiet=0):
+def seqio_merge_multiple(files, filetype):
'''
- Create a file containing the reverse complement of a file of small RNA
+ Merge SeqIO.parse over multiple files
'''
+ for filepath in files:
+ for seq in SeqIO.parse(filepath, filetype):
+ yield seq
- do_log(quiet, '====> Reverse complimenting RNA...')
- PATH_TO_REVCOMP = os.path.join(get_config_key('general', 'output_directory'), 'revcomp_rna.fastq')
+def seqio_revcomp_multiple(fastq_files):
+ '''
+ Merge SeqIO.parse over multiple files and reverse complement sequences
+ '''
+ seen_ids = set()
+
+ for seq in seqio_merge_multiple(fastq_files, 'fastq'):
+ if seq.id in seen_ids:
+ print(f'Error: duplicate ID {seq.id} found in your small RNA files, make IDs unique and try again')
+ exit(1)
+ else:
+ seen_ids.add(seq.id)
- def do_revcomp(seq):
revcomp = seq.reverse_complement()
revcomp.id = seq.id
- return revcomp
+ yield revcomp
+
+def merge_and_revcomp_input_files(smallRNA_files, output):
+ '''
+ Create a file containing the reverse complement of a file of small RNA
+ '''
- if not check_ids_unique(smallRNA, 'fastq'):
- print('Error: duplicate ID found in your small RNA file, make IDs unique and try again')
- exit(1)
+ print('====> Reverse complimenting RNA...')
+ PATH_TO_REVCOMP = os.path.join(output, 'revcomp_rna.fastq')
- seqs = SeqIO.parse(smallRNA, 'fastq')
- revcomps = map(do_revcomp, seqs)
- SeqIO.write(revcomps, PATH_TO_REVCOMP, 'fastq')
+ SeqIO.write(seqio_revcomp_multiple(smallRNA_files), PATH_TO_REVCOMP, 'fastq')
return PATH_TO_REVCOMP
-
-def find_targets(smallRNA, possible_target_list, threads=4, min_seq_length=2, mismatches_allowed=0, quiet=0):
+def find_targets(smallRNA, possible_target_list, output, program_paths, threads=4, min_seq_length=2, mismatches_allowed=0, verbose=False):
'''
Align the small RNA against the lists of possible targets with bowtie2 and
analyse the output
'''
CWD = os.getcwd()
- INDEX_DIR = os.path.join(get_config_key('general', 'output_directory'), 'bowtie_indexes')
- SAM_DIR = os.path.join(CWD, get_config_key('general', 'output_directory'), 'target_alignments')
+ INDEX_DIR = os.path.join(output, 'bowtie_indexes')
+ SAM_DIR = os.path.join(CWD, output, 'target_alignments')
sam_files = []
- mkdir_if_not_exists(SAM_DIR)
- mkdir_if_not_exists(INDEX_DIR)
+ if not os.path.exists(SAM_DIR):
+ os.makedirs(SAM_DIR)
+
+ if not os.path.exists(INDEX_DIR):
+ os.makedirs(INDEX_DIR)
+
os.chdir(INDEX_DIR)
for target in possible_target_list:
@@ -82,37 +141,34 @@ def find_targets(smallRNA, possible_target_list, threads=4, min_seq_length=2, mi
print(f'Error: duplicate ID found in your target file: {target}. Make IDs unique and try again')
exit(1)
- do_log(quiet, f'====> Builing index for {target}')
+ print(f'====> Building index for {target}...')
INDEX_NAME = os.path.basename(target) + '_index'
bowtie2_build_command = [
- get_config_key('cli-tools', 'bowtie2', 'path_to_bowtie2_build'),
+ program_paths['bowtie2-build'],
'--threads', str(threads),
os.path.join(CWD, target),
INDEX_NAME
]
- bowtie2_build_command = bowtie2_build_command + get_config_key('cli-tools', 'bowtie2', 'bowtie2_build_params')
-
- run(bowtie2_build_command, capture_output=(quiet != 0))
+ run(bowtie2_build_command, capture_output=not verbose)
- do_log(quiet, f'====> Aligning small RNA against {target}')
+ print(f'====> Aligning small RNA against {target}...')
sam_files.append(os.path.join(SAM_DIR, INDEX_NAME + '.sam'))
if mismatches_allowed > 0:
bowtie2_align_command = [
- get_config_key('cli-tools', 'bowtie2', 'path_to_bowtie2'),
+ program_paths['bowtie2'],
'--threads', str(threads),
'-L', str(min_seq_length),
- '--no-1mm-upfront',
'--score-min', 'L,-' + str(mismatches_allowed) + ',0',
'--end-to-end',
'--norc',
'--mp', '1,1',
'--ignore-quals',
- '--rdg', '9,1',
- '--rfg', '9,1',
+ '--rdg', '100,100',
+ '--rfg', '100,100',
'-x', INDEX_NAME,
'-U', os.path.join(CWD, smallRNA),
'-S', sam_files[-1],
@@ -120,7 +176,7 @@ def find_targets(smallRNA, possible_target_list, threads=4, min_seq_length=2, mi
]
else:
bowtie2_align_command = [
- get_config_key('cli-tools', 'bowtie2', 'path_to_bowtie2'),
+ program_paths['bowtie2'],
'--threads', str(threads),
'-L', str(min_seq_length),
'--no-1mm-upfront',
@@ -134,33 +190,30 @@ def find_targets(smallRNA, possible_target_list, threads=4, min_seq_length=2, mi
'-a'
]
- bowtie2_align_command = bowtie2_align_command + get_config_key('cli-tools', 'bowtie2', 'bowtie2_params')
-
- run(bowtie2_align_command, capture_output=(quiet != 0))
+ run(bowtie2_align_command, capture_output=not verbose)
os.chdir(CWD)
return sam_files
-def build_summary_files(sam_files, quiet=0):
+def build_summary_files(sam_files, output, program_paths, verbose=False):
'''
Takes the alignment SAM and extracts the sequences into a summary file and
a set of FASTA files
'''
-
- do_log(quiet, '====> Removing unaligend sequences and building summary')
-
- with open(os.path.join(get_config_key('general', 'output_directory'), 'rna_target_list.tsv'), 'w') as tablefile:
- writer = csv.writer(tablefile, delimiter='\t')
- writer.writerow(['Query Sequence', 'Target File', 'Target Sequence', 'Start Coordinate', 'End Coordinate', 'Strand'])
+ print('====> Removing unaligend sequences and building summary files...')
- for filename in sam_files:
+ target_rows = defaultdict(lambda: [])
+ # Create a set of targets to analyse for enrichment later
+ target_set = set()
+ for filename in sam_files:
samfile = pysam.AlignmentFile(filename, 'r')
fileheader = samfile.header
for read in samfile:
if not read.is_unmapped:
query_name = read.query_name
+ query_sequence = Seq(read.query_sequence).reverse_complement()
r_name = fileheader.get_reference_name(read.reference_id)
start = read.reference_start
end = read.reference_end
@@ -170,34 +223,205 @@ def build_summary_files(sam_files, quiet=0):
else:
strand = 'sense'
- writer.writerow([query_name, os.path.abspath(filename), r_name, start, end, strand])
-
- try:
- remove_duplicate_command = [
- get_config_key('cli-tools', 'samtools', 'path_to_samtools'),
- 'view',
- '-o', filename + '.nodup.sam',
- filename
- ]
-
- remove_duplicate_command = remove_duplicate_command + get_config_key('cli-tools', 'samtools', 'samtools_view_params')
-
- result = run(remove_duplicate_command)
- result.check_returncode()
- filename = filename + '.nodup.sam'
-
- bam_to_fastq_command = [
- get_config_key('cli-tools', 'samtools', 'path_to_samtools'),
- 'fastq',
- filename
- ]
-
- bam_to_fastq_command = bam_to_fastq_command + get_config_key('cli-tools', 'samtools', 'samtools_fastq_params')
-
- result = run(bam_to_fastq_command, capture_output=True)
- result.check_returncode()
- with open(filename + '.fastq', 'wb') as f:
- f.write(result.stdout)
-
- except:
- print('SAMTOOLS convertion failed, skipping creating filtered and FASTQ files...')
\ No newline at end of file
+ length = len(query_sequence)
+ first_base = query_sequence[0].replace('T', 'U')
+ target_rows[f'{length}{first_base}'].append([query_name, query_sequence, os.path.abspath(filename), r_name, start, end, strand])
+ target_set.add(r_name)
+
+ TRAGETS_DIR = os.path.join(output, 'targets_by_group')
+
+ if not os.path.exists(TRAGETS_DIR):
+ os.makedirs(TRAGETS_DIR)
+
+ for target in target_rows:
+ with open(os.path.join(TRAGETS_DIR, f'{target}_targets.tsv'), 'w') as f:
+ writer = csv.writer(f, delimiter='\t')
+ writer.writerow(['Query Name', 'Query Sequence', 'Target File', 'Target Sequence', 'Start Coordinate', 'End Coordinate', 'Strand'])
+ writer.writerows(target_rows[target])
+
+ with open(os.path.join(output, 'rna_target_list.tsv'), 'w') as tablefile:
+ writer = csv.writer(tablefile, delimiter='\t')
+ writer.writerow(['Query Name', 'Query Sequence', 'Target File', 'Target Sequence', 'Start Coordinate', 'End Coordinate', 'Strand'])
+ for target in target_rows:
+ writer.writerows(target_rows[target])
+ try:
+ remove_duplicate_command = [
+ program_paths['samtools'],
+ 'view',
+ '-h',
+ '-F', '256',
+ '-F', '4',
+ '-o', filename + '.nodup.sam',
+ filename
+ ]
+
+ result = run(remove_duplicate_command, capture_output=not verbose)
+ result.check_returncode()
+ filename = filename + '.nodup.sam'
+
+ bam_to_fastq_command = [
+ program_paths['samtools'],
+ 'fastq',
+ filename
+ ]
+
+ result = run(bam_to_fastq_command, capture_output=True)
+ result.check_returncode()
+ with open(filename + '.fastq', 'wb') as f:
+ f.write(result.stdout)
+
+ except:
+ print('====> SAMTOOLS convertion failed, skipping creating filtered and FASTQ files...')
+
+ return target_set
+
+def do_term_enrichment(files_to_annotate, target_list, path_to_emapper_data, output, program_paths, threads=4, verbose=False):
+ '''
+ Use eggnog-mapper to annotate GO, KEGG and PFam with diamond blastx to provide a quick overview enrichment of terms
+ '''
+ EGGNOG_MAPPER_OUTPUT = os.path.join(output, 'eggnog_mapper')
+ COMBINED_TARGETS_FILE = os.path.join(output, 'eggnog_mapper', 'targets.fasta')
+
+ if not os.path.exists(EGGNOG_MAPPER_OUTPUT):
+ os.makedirs(EGGNOG_MAPPER_OUTPUT)
+
+ print('====> Merging target files...')
+ SeqIO.write(seqio_merge_multiple(files_to_annotate, 'fasta'), COMBINED_TARGETS_FILE, 'fasta')
+
+ print('====> Annotating with eggnog-mapper...')
+ eggnog_mapper_command = [
+ program_paths['eggnog-mapper'], '-i', COMBINED_TARGETS_FILE, '-o', 'eggnog_mapper', '--output_dir', EGGNOG_MAPPER_OUTPUT, '--cpu', str(threads),
+ '--data_dir', path_to_emapper_data, '--itype', 'CDS'
+ ]
+ run(eggnog_mapper_command, capture_output=not verbose)
+
+ print('====> Reading eggnog-mapper output...')
+ go_terms = {}
+
+ kegg_pathways = {}
+ kegg_counts = defaultdict(lambda: 0)
+ total_kegg_annot = 0
+
+ pfams = {}
+ pfam_counts = defaultdict(lambda: 0)
+ total_pfam_annot = 0
+ with open(os.path.join(EGGNOG_MAPPER_OUTPUT, 'eggnog_mapper.emapper.annotations')) as f:
+ reader = csv.reader(f, delimiter='\t')
+
+ for line in reader:
+ if not line[0].startswith('#'):
+ if line[9] == '-':
+ go_terms[line[0]] = set()
+ else:
+ go_terms[line[0]] = set(line[9].split(','))
+
+ kegg_pathways[line[0]] = set()
+ for kegg in line[12].split(','):
+ if 'map' in kegg:
+ kegg_pathways[line[0]].add(kegg)
+ kegg_counts[kegg] += 1
+ total_kegg_annot += 1
+
+ if line[-1] == '-':
+ pfams[line[0]] = set()
+ else:
+ pfam_set = set(line[-1].split(','))
+ pfams[line[0]] = pfam_set
+
+ total_pfam_annot += len(pfam_set)
+
+ for pfam in pfam_set:
+ pfam_counts[pfam] += 1
+
+ with open(os.path.join(EGGNOG_MAPPER_OUTPUT, 'go_map.txt'), 'w') as f:
+ writer = csv.writer(f, delimiter='\t')
+ for name in go_terms:
+ writer.writerow([name, ','.join(go_terms[name])])
+
+ with open(os.path.join(EGGNOG_MAPPER_OUTPUT, 'kegg_map.txt'), 'w') as f:
+ writer = csv.writer(f, delimiter='\t')
+ for name in kegg_pathways:
+ writer.writerow([name, ','.join(kegg_pathways[name])])
+
+ with open(os.path.join(EGGNOG_MAPPER_OUTPUT, 'pfam_map.txt'), 'w') as f:
+ writer = csv.writer(f, delimiter='\t')
+ for name in pfams:
+ writer.writerow([name, ','.join(pfams[name])])
+
+ genes_to_analyse_string = ''
+ for gene in target_list:
+ genes_to_analyse_string += f',"{gene}"'
+
+ print('====> Enriching GO terms using TopGO...')
+ customised_go_enrichment = GO_ENRICHMENT_SCRIPT.replace(
+ '{{INTERSTING_GENE_LIST}}', genes_to_analyse_string[1:], 1
+ ).replace(
+ '{{PATH_TO_GO_TXT}}', os.path.join(EGGNOG_MAPPER_OUTPUT, 'go_map.txt'), 1
+ ).replace(
+ '{{PATH_TO_OUTPUT_TSV}}', os.path.join(output, 'targets_enriched_go_terms.tsv'), 1
+ )
+
+ GO_ENRICHMENT_PATH = os.path.join(EGGNOG_MAPPER_OUTPUT, 'enrichGOTerms.R')
+ with open(GO_ENRICHMENT_PATH, 'w') as f:
+ f.write(customised_go_enrichment)
+
+ go_enrich_command = [
+ program_paths['rscript'], GO_ENRICHMENT_PATH
+ ]
+ run(go_enrich_command, capture_output=not verbose)
+
+ print('====> Enriching KEGG pathways with scipy fisher_exact...')
+ kegg_output_table = []
+
+ kegg_in_set = []
+ for interesting_gene in target_list:
+ if interesting_gene in kegg_pathways:
+ kegg_in_set += kegg_pathways[interesting_gene]
+
+ kegg_in_total = len(kegg_in_set)
+ kegg_counter = Counter(kegg_in_set)
+
+ for kegg in kegg_counter:
+ in_intersting = kegg_counter[kegg]
+ out_intesting = kegg_counts[kegg] - in_intersting
+
+ in_not_intersting = kegg_in_total - in_intersting
+ out_not_intersting = total_kegg_annot - out_intesting - in_intersting - in_not_intersting
+
+ fisher_table = [[in_intersting, out_intesting], [in_not_intersting, out_not_intersting]]
+ result = fisher_exact(fisher_table, alternative='greater')
+
+ if result.pvalue < 0.01:
+ kegg_output_table.append([kegg, KEGG_PATHWAY_NAMES[kegg], result.pvalue])
+
+ with open(os.path.join(output, 'targets_enriched_kegg_pathways.tsv'), 'w') as f:
+ writer = csv.writer(f, delimiter='\t')
+ writer.writerows(kegg_output_table)
+
+ print('====> Enriching PFams with scipy fisher_exact...')
+ pfam_output_table = []
+
+ pfam_in_set = []
+ for interesting_gene in target_list:
+ if interesting_gene in pfams:
+ pfam_in_set += pfams[interesting_gene]
+
+ pfam_in_total = len(pfam_in_set)
+ pfam_counter = Counter(pfam_in_set)
+
+ for pfam in pfam_counter:
+ in_intersting = pfam_counter[pfam]
+ out_interesting = pfam_counts[pfam] - in_intersting
+
+ in_not_intersting = pfam_in_total - in_intersting
+ out_not_intersting = total_pfam_annot - in_intersting - out_interesting - in_not_intersting
+
+ fisher_table = [[in_intersting, out_interesting], [in_not_intersting, out_not_intersting]]
+ result = fisher_exact(fisher_table, alternative='greater')
+
+ if result.pvalue < 0.01:
+ pfam_output_table.append([pfam, result.pvalue])
+
+ with open(os.path.join(output, 'targets_enriched_pfams.tsv'), 'w') as f:
+ writer = csv.writer(f, delimiter='\t')
+ writer.writerows(pfam_output_table)
\ No newline at end of file
diff --git a/hlsmallrna/trim.py b/hlsmallrna/trim.py
index 5c24ac9..8cebf21 100644
--- a/hlsmallrna/trim.py
+++ b/hlsmallrna/trim.py
@@ -1,4 +1,4 @@
-# Copyright 2022 Vicky Hunt Lab Members
+# Copyright 2022 - 2025 Vicky Hunt Lab Members
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
@@ -11,48 +11,57 @@
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
-import os.path
+import os
+import sys
+import shutil
from subprocess import run
-from .config import do_log, get_config_key
-
-def run_trim(file_to_trim, start_adapter, end_adapter, both_adapters, quiet=0):
+def run_trim(
+ file_to_trim, output, five_prime_adapter=None, three_prime_adapter=None, quality_cutoff=20, threads=4,
+ cutadapt_path='cutadapt', keep_intermediate=False, verbose=False
+ ):
'''
Call cutadapt to remove adapters
'''
+ current_input = file_to_trim
+ files_to_remove = []
- outfile = os.path.join(get_config_key('general', 'output_directory'), 'trimmed_rna.fastq')
- command = [
- get_config_key('cli-tools', 'trim', 'path_to_trim'),
- file_to_trim
- ]
-
- if start_adapter is not None:
- command.append('-a')
- command.append(start_adapter)
-
- if end_adapter is not None:
- command.append('-g')
- command.append(end_adapter)
-
- if both_adapters is not None:
- command.append('-b')
- command.append(both_adapters)
+ # quality filter
+ if quality_cutoff > 0:
+ print('====> Removing low quality reads...', file=sys.stderr)
+ cutadapt_quality = [
+ cutadapt_path, '-q', str(quality_cutoff), '-o', os.path.join(output, 'quality_filtered.fq'), '-j', str(threads), current_input
+ ]
+ run(cutadapt_quality, capture_output=not verbose)
- if get_config_key('cli-tools', 'trim', 'trim_pass_threads'):
- threads = get_config_key('general', 'threads')
+ current_input = os.path.join(output, 'quality_filtered.fq')
+ files_to_remove.append(current_input)
- command.append('-j')
- command.append(str(threads))
+ # trim five prime adapters
+ if five_prime_adapter is not None:
+ print('====> Removing complete 5\' adapters...', file=sys.stderr)
+ cutadapt_five_prime = [
+ cutadapt_path, '-g', five_prime_adapter, '-o', os.path.join(output, 'five_prime_removed.fq'), '-j', str(threads), current_input
+ ]
+ run(cutadapt_five_prime, capture_output=not verbose)
- command = command + get_config_key('cli-tools', 'trim', 'trim_params')
+ current_input = os.path.join(output, 'five_prime_removed.fq')
+ files_to_remove.append(current_input)
- do_log(quiet, '====> Removing adpters from ends')
- do_log(quiet, command)
- with open(outfile, 'w') as f:
- result = run(command, stdout=f)
+ # trim three prime adapters and remove reads without
+ if three_prime_adapter is not None:
+ print('====> Removing 3\' adapters and filtering reads...', file=sys.stderr)
+ cutadapt_three_prime = [
+ cutadapt_path, '-a', three_prime_adapter, '-m', '1', '--discard-untrimmed', '-o', os.path.join(output, 'three_prime_removed.fq'),
+ '-j', str(threads), current_input
+ ]
+ run(cutadapt_three_prime, capture_output=not verbose)
- result.check_returncode()
+ current_input = os.path.join(output, 'three_prime_removed.fq')
+ files_to_remove.append(current_input)
- return result.returncode == 0
\ No newline at end of file
+ shutil.copy2(current_input, os.path.join(output, 'trimmed_reads.fq'))
+ if not keep_intermediate:
+ for filename in files_to_remove:
+ os.remove(filename)
\ No newline at end of file
diff --git a/hlsmallrna/unitas.py b/hlsmallrna/unitas.py
index 3ee55b2..da23217 100644
--- a/hlsmallrna/unitas.py
+++ b/hlsmallrna/unitas.py
@@ -1,4 +1,4 @@
-# Copyright 2022 Vicky Hunt Lab Members
+# Copyright 2022 - 2025 Vicky Hunt Lab Members
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
@@ -17,6 +17,7 @@
import csv
import re
+from math import inf
from subprocess import run
from collections import defaultdict
@@ -25,36 +26,41 @@
from matplotlib import pyplot as plt
from Bio import SeqIO
-from .config import do_log, get_config_key
-
-def copy_and_label_file(ref_file, filename, label):
+def find_min_max_lengths(unitas_in):
'''
- Labeling function for the CDS and Unspliced Transcriptome files
+ Calculate the minimum and maximum lengths for the unitas report based on fastq filenames
'''
+ minimum = inf
+ maximum = -inf
- def rename_iter(path):
- for seq in SeqIO.parse(path, 'fasta'):
- seq.id = f'{label}|{seq.id}'
+ filenames = glob.glob(os.path.join(unitas_in, '*.fq')) + glob.glob(os.path.join(unitas_in, '*.fastq'))
+ for name in filenames:
+ length = int(re.findall(r'length(\d+)', name)[0])
- yield seq
+ if length > maximum:
+ maximum = length
+
+ if length < minimum:
+ minimum = length
-
- new_path = os.path.join(get_config_key('general', 'output_directory'), f'{filename}.fasta')
- SeqIO.write(rename_iter(ref_file), new_path, 'fasta')
+ if maximum == -inf:
+ maximum = 30
+ if minimum == inf:
+ minimum = 18
- return new_path
+ return minimum, maximum
-def run_unitas_annotation(small_rna, species_name, ref_files, quiet=0, unitas_output='.'):
+def run_unitas_annotation(species_name, ref_files, program_paths, output, verbose, small_rna):
'''
Run Unitas on the small RNAs, given a file and species name
'''
CWD = os.getcwd()
SMALL_RNA_PATH = os.path.join(CWD, small_rna)
- os.chdir(unitas_output)
+ os.chdir(output)
unitas_command = [
- get_config_key('cli-tools', 'unitas', 'path_to_unitas'),
+ program_paths['unitas'],
'-input', SMALL_RNA_PATH,
'-species', species_name
]
@@ -64,138 +70,121 @@ def run_unitas_annotation(small_rna, species_name, ref_files, quiet=0, unitas_ou
unitas_command.append('-refseq')
unitas_command.append(os.path.join(CWD, refs))
- unitas_command = unitas_command + get_config_key('cli-tools', 'unitas', 'unitas_params')
-
- if get_config_key('cli-tools', 'unitas', 'unitas_pass_threads'):
- threads = get_config_key('general', 'threads')
-
- unitas_command.append('-threads')
- unitas_command.append(str(threads))
-
- do_log(quiet, '==> Running first unitas pass')
- result = run(unitas_command, capture_output=(quiet != 0))
+ print(f'====> Running unitas on {os.path.basename(small_rna)}...')
+ result = run(unitas_command, capture_output=not verbose)
result.check_returncode()
os.chdir(CWD)
-def merge_summary():
+def merge_summary(output):
'''
Merge together the summery files into one CSV
'''
- new_file = []
+ numbered_names = defaultdict(lambda: 'Unnamed')
+ numbered_results = defaultdict(lambda: defaultdict(lambda: 0))
+ named_results = defaultdict(lambda: defaultdict(lambda: 0))
+ unitas_order = defaultdict(lambda: [])
+
+ for filename in glob.glob(os.path.join(output, 'unitas', '*/unitas.annotation_summary.txt')):
+ length_regex_result = re.findall(r'length(\d+)', filename)
+ if len(length_regex_result) < 1:
+ namestr = f'File Name : {filename}'
+ else:
+ length = int(length_regex_result[0])
+ namestr = f'RNA Length {length}'
+ numbered_names[length] = namestr
+
+ namestr = length
- numbered_results = {}
- named_results = {}
- for filename in glob.glob(os.path.join(get_config_key('general', 'output_directory'), 'unitas', '*/unitas.annotation_summary.txt')):
with open(filename) as f:
reader = csv.reader(f, delimiter='\t')
- length = re.findall(r'length(\d+)', filename)
- if len(length) < 1:
- namestr = f'File Name : {filename}'
- else:
- namestr = f'RNA Length {length[0]}'
-
- names = []
- values = []
for line in reader:
- names.append(line[0])
- values.append(line[1])
+ if type(namestr) == int:
+ numbered_results[namestr][line[0]] = float(line[1])
+ else:
+ named_results[namestr][line[0]] = float(line[1])
- if len(length) > 1:
- numbered_results[length[1]] = [namestr, names, values]
- else:
- named_results[namestr] = [names, values]
+ unitas_order[namestr].append(line[0])
- for result in sorted(numbered_results):
- new_file.append([numbered_results[result][0]])
- new_file.append(numbered_results[result][1])
- new_file.append(numbered_results[result][2])
+ combined_unitas_file = []
+ for result in sorted(numbered_results.keys()):
+ names = []
+ values = []
- for result in sorted(named_results):
- new_file.append([result])
- new_file.append(named_results[result][0])
- new_file.append(named_results[result][1])
+ for unitas_cat in unitas_order[result]:
+ names.append(unitas_cat)
+ values.append(numbered_results[result][unitas_cat])
- unitas_table = os.path.join(get_config_key('general', 'output_directory'), 'unitas_summary.csv')
- with open(unitas_table, 'w') as f:
- writer = csv.writer(f)
+ combined_unitas_file.append([numbered_names[result]])
+ combined_unitas_file.append(names)
+ combined_unitas_file.append(values)
+ combined_unitas_file.append([])
+
+ for result in sorted(named_results.keys()):
+ names = []
+ values = []
+
+ for unitas_cat in unitas_order[result]:
+ names.append(unitas_cat)
+ values.append(numbered_results[result][unitas_cat])
- writer.writerows(new_file)
+ combined_unitas_file.append([result])
+ combined_unitas_file.append(names)
+ combined_unitas_file.append(values)
+ combined_unitas_file.append([])
- return unitas_table
+ unitas_table = os.path.join(output, 'unitas_summary.csv')
+ with open(unitas_table, 'w') as f:
+ writer = csv.writer(f)
+ writer.writerows(combined_unitas_file)
+
+ return numbered_results
-def graph_unitas_classification_type(path_to_table):
+def graph_unitas_classification_type(numbered_results, min_length, max_length, output):
'''
Craete a graph of small RNA by length against unitas type
'''
- with open(path_to_table) as f:
- reader = csv.reader(f)
-
- fig, ax = plt.subplots(figsize=(11.69, 8.27))
-
- labels = set()
- values = defaultdict(lambda: defaultdict(lambda: 0))
-
- # Use skip_group to skip any group not named by length
- skip_group = False
- for i, line in enumerate(reader):
- if i % 3 == 0:
- if not line[0].startswith('RNA Length '):
- skip_group = True
- else:
- line_name = re.findall('\d+', line[0])[0]
- elif i % 3 == 1:
- if skip_group:
- continue
-
- line_labels = set()
- for l in line:
- if not l[0].isspace():
- line_labels.add(l)
-
- all_labels = line
-
- labels = labels | line_labels
- elif i % 3 == 2:
- if skip_group:
- skip_group = False
- continue
-
- for j, name in enumerate(all_labels):
- if name in line_labels:
- values[line_name][name] = float(line[j])
-
- with open(os.path.join(get_config_key('general', 'output_directory'), 'unitas_graph_data.csv'), 'w') as f:
- writer = csv.writer(f)
-
- rna_lengths = list(values.keys())
- writer.writerow(['RNA Length'] + rna_lengths)
-
- offsets = None
- for l in labels:
- counts = []
-
- for key in rna_lengths:
- counts.append(values[key][l])
-
- ax.bar(rna_lengths, counts, label=l, bottom=offsets)
+ fig, ax = plt.subplots(figsize=(11.69, 8.27))
+ final_cats = ['gene', '5_UTR', '3_UTR', 'miRNA', 'tRNA', 'rRNA', 'TE', 'no annotation', 'low_complexity']
+ rna_lengths = list(range(min_length, max_length + 1))
+
+ # figure out if there are any other toplevel categories we need to include on the graph
+ included_cats = set(final_cats)
+ for length in numbered_results:
+ for category in numbered_results[length]:
+ if not category.startswith(' '):
+ included_cats.add(category)
+
+ for label in included_cats:
+ if label not in final_cats:
+ final_cats.append(label)
+
+ # Plot each caegory on the graph and add to the plot data table
+ base_offset = np.array([0] * len(rna_lengths))
+ plot_data = [['RNA Length'] + rna_lengths]
+ for category in final_cats:
+ category_bars = []
+ for length in rna_lengths:
+ if category in numbered_results[length]:
+ category_bars.append(numbered_results[length][category])
+ else:
+ category_bars.append(0)
- writer.writerow([l] + counts)
+ ax.bar(rna_lengths, category_bars, bottom=base_offset, label=category, align='edge', width=0.8)
+ plot_data.append([category] + category_bars)
+ base_offset = base_offset + np.array(category_bars)
- if offsets is None:
- offsets = np.array(counts)
- else:
- offsets = offsets + np.array(counts)
+ with open(os.path.join(output, 'unitas_graph_data.csv'), 'w') as f:
+ writer = csv.writer(f)
+ writer.writerows(plot_data)
+ ax.set_xticks(rna_lengths)
ax.set_xticklabels(rna_lengths, fontsize=7)
ax.set_xlabel('Length of small RNA')
- ax.set_ylabel('')
+ ax.set_ylabel('Unitas RPM')
plt.legend(title='Unitas Category')
- plt.savefig(os.path.join(get_config_key('general', 'output_directory'), 'unitasGraph.png'))
-
-if __name__ == '__main__':
- graph_unitas_classification_type('output/unitas_summery.csv')
\ No newline at end of file
+ plt.savefig(os.path.join(output, 'unitasGraph.png'))
\ No newline at end of file
diff --git a/hlsmallrna/validate_config.py b/hlsmallrna/validate_config.py
new file mode 100644
index 0000000..918cde1
--- /dev/null
+++ b/hlsmallrna/validate_config.py
@@ -0,0 +1,314 @@
+# Copyright 2022 - 2025 Vicky Hunt Lab Members
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+# http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+import os
+import gzip
+import shutil
+
+import yaml
+
+from Bio import SeqIO
+
+KNOWN_ADAPTERS = {
+ 'qiagen': {
+ '5': 'GTTCAGAGTTCTACAGTCCGACGATC',
+ '3': 'AACTGTAGGCACCATCAAT'
+ },
+ 'qiagen-rev': {
+ '5': 'AACTGTAGGCACCATCAAT',
+ '3': 'GTTCAGAGTTCTACAGTCCGACGATC'
+ },
+ 'nebnext': {
+ '5': None,
+ '3': 'AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC'
+ },
+ 'nebnext-rev': {
+ '5': None,
+ '3': 'GATCGTCGGACTGTAGAACTCTGAACGTGTAGATCTCGGTGGTCGCCGTATCATT'
+ }
+}
+
+def make_path_reletive_to_config(path, config):
+ '''
+ Convert a path to be reletive to the config directory
+ '''
+ return os.path.join(os.path.dirname(config), path)
+
+def unzip_if_needed(filepath, output_dir, tmpfilename='temp.fq'):
+ '''
+ Check if filepath is gzipped
+ '''
+ with open(filepath, 'rb') as test_f:
+ is_gzipped = test_f.read(2) == b'\x1f\x8b'
+
+ if not is_gzipped:
+ return filepath
+ else:
+ with gzip.open(filepath, 'rb') as file_in:
+ with open(os.path.join(output_dir, tmpfilename), 'wb') as file_out:
+ shutil.copyfileobj(file_in, file_out)
+
+ return os.path.join(output_dir, tmpfilename)
+
+
+def fasta_to_fastq_record(seqs):
+ '''
+ remove @s from the sequence IDs of a biopython object, replace with _ and
+ add phread scores of 40
+ '''
+ for seq in seqs:
+ if '@' in seq.id:
+ seq.id = seq.id.replace('@', '_')
+
+ seq.letter_annotations["phred_quality"] = [40] * len(seq)
+
+ yield seq
+
+def convert_to_fastq_if_needed(filepath, output_dir, tmpfilename='temp.fa'):
+ '''
+ Convert a FASTA file provided as input to a FASTQ if needed
+ '''
+ # if it is a FASTA file the first character is >
+ with open(filepath, 'r') as f:
+ is_fasta = f.read(1) == '>'
+
+ if not is_fasta:
+ return filepath
+ else:
+ SeqIO.write(
+ fasta_to_fastq_record(SeqIO.parse(filepath, 'fasta')),
+ os.path.join(output_dir, tmpfilename), 'fastq'
+ )
+ return os.path.join(output_dir, tmpfilename)
+
+def load_and_validate_input(input_yaml_path, output_dir):
+ '''
+ Check for errors in the small RNA pipeline input YAML
+ '''
+ # TODO: check file types and existance
+ with open(input_yaml_path) as f:
+ input_yaml = yaml.load(f.read(), Loader=yaml.Loader)
+
+ smallRNA_fastq_present = 'smallRNA_fastq' in input_yaml
+ size_sorted_fastqs_present = 'size_sorted_fastqs' in input_yaml
+
+ # smallRNA_fastq and size_sorted_fastqs are mutaually exclusive in the input file
+ if (smallRNA_fastq_present and size_sorted_fastqs_present) or not (smallRNA_fastq_present or size_sorted_fastqs_present):
+ return False, 'Exactly one of smallRNA_fastq or size_sorted_fastqs must be specified in the config file', input_yaml
+
+ if smallRNA_fastq_present:
+ input_yaml['smallRNA_fastq'] = convert_to_fastq_if_needed(
+ unzip_if_needed(
+ make_path_reletive_to_config(input_yaml['smallRNA_fastq'], input_yaml_path),
+ output_dir, tmpfilename='smallRNA.fastq'
+ ),
+ output_dir, tmpfilename='smallRNA.fasta.fastq'
+ )
+ else:
+ input_yaml['size_sorted_fastqs'] = make_path_reletive_to_config(input_yaml['size_sorted_fastqs'], input_yaml_path)
+
+ # cds and unspliced_transcriptome set to None if not present
+ if 'cds' not in input_yaml:
+ input_yaml['cds'] = None
+ else:
+ input_yaml['cds'] = unzip_if_needed(make_path_reletive_to_config(input_yaml['cds'], input_yaml_path), output_dir, tmpfilename='cds.fasta')
+
+ if 'unspliced_transcriptome' not in input_yaml:
+ input_yaml['unspliced_transcriptome'] = None
+ else:
+ input_yaml['unspliced_transcriptome'] = unzip_if_needed(
+ make_path_reletive_to_config(input_yaml['unspliced_transcriptome'], input_yaml_path), output_dir, tmpfilename='unspliced.fasta'
+ )
+
+ # program paths are all optional, default to in path
+ if 'program_paths' not in input_yaml:
+ input_yaml['program_paths'] = {}
+
+ if 'unitas' not in input_yaml['program_paths']:
+ input_yaml['program_paths']['unitas'] = 'unitas.pl'
+
+ if 'bowtie2-build' not in input_yaml['program_paths']:
+ input_yaml['program_paths']['bowtie2-build'] = 'bowtie2-build'
+
+ if 'bowtie2' not in input_yaml['program_paths']:
+ input_yaml['program_paths']['bowtie2'] = 'bowtie2'
+
+ if 'cutadapt' not in input_yaml['program_paths']:
+ input_yaml['program_paths']['cutadapt'] = 'cutadapt'
+
+ if 'samtools' not in input_yaml['program_paths']:
+ input_yaml['program_paths']['samtools'] = 'samtools'
+
+ if 'bedtools' not in input_yaml['program_paths']:
+ input_yaml['program_paths']['bedtools'] = 'bedtools'
+
+ if 'eggnog-mapper' not in input_yaml['program_paths']:
+ input_yaml['program_paths']['eggnog-mapper'] = 'emapper.py'
+
+ if 'rscript' not in input_yaml['program_paths']:
+ input_yaml['program_paths']['rscript'] = 'Rscript'
+
+ # trim
+ if 'trim' in input_yaml:
+ # if trim is specified, require smallRNA_fastq
+ if not smallRNA_fastq_present:
+ return False, 'For trim and sort steps smallRNA_fastq must be specified, not size_sorted_fastqs', input_yaml
+ else:
+ input_yaml['trim']['input'] = input_yaml['smallRNA_fastq']
+
+ # kit is mually exclusive to 5_prime and 3_prime
+ if ('5_prime' in input_yaml['trim'] or '3_prime' in input_yaml['trim']) and 'kit' in input_yaml['trim']:
+ return False, 'If a kit is specified in the trim config, 5_prime and 3_prime sequences should not be', input_yaml
+
+ # kit must be one of qiagen, qiagen-rev, neb or neb-rev
+ if 'kit' in input_yaml['trim'] and input_yaml['trim']['kit'] not in KNOWN_ADAPTERS.keys():
+ return False, f'The kit specified in trim should be one of: {", ".join(KNOWN_ADAPTERS.keys())}', input_yaml
+ elif 'kit' in input_yaml['trim']:
+ input_yaml['trim']['5_prime'] = KNOWN_ADAPTERS[input_yaml['trim']['kit']]['5']
+ input_yaml['trim']['3_prime'] = KNOWN_ADAPTERS[input_yaml['trim']['kit']]['3']
+
+ del input_yaml['trim']['kit']
+
+ # min_quality set to 20 if it doesn't exist
+ if 'min_quality' not in input_yaml['trim']:
+ input_yaml['trim']['min_quality'] = 20
+
+ # sort
+ if 'sort' in input_yaml:
+ if not smallRNA_fastq_present:
+ return False, 'For trim and sort steps smallRNA_fastq must be specified, not size_sorted_fastqs', input_yaml
+ elif 'trim' in input_yaml:
+ input_yaml['sort']['input'] = os.path.join(output_dir, 'trimmed_reads.fq')
+ else:
+ input_yaml['sort']['input'] = input_yaml['smallRNA_fastq']
+
+ if 'genome' in input_yaml['sort']:
+ input_yaml['sort']['genome'] = unzip_if_needed(
+ make_path_reletive_to_config(input_yaml['sort']['genome'], input_yaml_path),
+ output_dir, tmpfilename='genome.fasta'
+ )
+
+ input_yaml['sort']['cds'] = input_yaml['cds']
+
+ # align_to_cds defaults to false, if true cds must be present
+ if 'align_to_cds' not in input_yaml['sort']:
+ input_yaml['sort']['align_to_cds'] = False
+ elif input_yaml['sort']['align_to_cds'] and input_yaml['sort']['cds'] is None:
+ return False, 'If align_to_cds is specified in the sort section, a cds file must be specified as well', input_yaml
+
+ # min_length defaults to 15
+ if 'min_length' not in input_yaml['sort']:
+ input_yaml['sort']['min_length'] = 15
+
+ # max_length defaults to 30
+ if 'max_length' not in input_yaml['sort']:
+ input_yaml['sort']['max_length'] = 30
+
+ # mismatches defaults to 0
+ if 'mismatches' not in input_yaml['sort']:
+ input_yaml['sort']['mismatches'] = 0
+
+ # unitas
+ if 'unitas' in input_yaml:
+ # if sort is not run size_sorted_fastqs must be specified
+ if 'sort' not in input_yaml and 'size_sorted_fastqs' not in input_yaml:
+ return False, 'To run unitas, either the sort step must be run first or size_sorted_fastqs provided instead of smallRNA_fastq', input_yaml
+ elif 'size_sorted_fastqs' in input_yaml:
+ input_yaml['unitas']['inputs'] = input_yaml['size_sorted_fastqs']
+ elif 'sort' in input_yaml:
+ input_yaml['unitas']['inputs'] = os.path.join(output_dir, 'binned_length_rna')
+
+ # at least one of refseq or species must be specified
+ if 'refseq' not in input_yaml['unitas'] and ('species' not in input_yaml['unitas'] or input_yaml['unitas']['species'] == 'x'):
+ return False, 'At least one refseq or non-x species must be provided in the unitas config, to reuse cds and unsplices transcriptiome, specify gene', input_yaml
+
+ if 'species' not in input_yaml['unitas']:
+ input_yaml['unitas']['species'] = 'x'
+
+ if 'refseq' in input_yaml['unitas'] and 'gene' in input_yaml['unitas']['refseq']:
+ gene_index = input_yaml['unitas']['refseq'].index('gene')
+ del input_yaml['unitas']['refseq'][gene_index]
+
+ if input_yaml['cds'] is None and input_yaml['unspliced_transcriptome'] is None:
+ return False, 'If gene is specified in unitas or targetid refseq or unspliced_transcriptome must also be specified', input_yaml
+
+ if input_yaml['cds'] is not None:
+ input_yaml['unitas']['refseq'].append({'gene': input_yaml['cds']})
+
+ if input_yaml['unspliced_transcriptome'] is not None:
+ input_yaml['unitas']['refseq'].append({'gene': input_yaml['unspliced_transcriptome']})
+
+ for i in range(len(input_yaml['unitas']['refseq'])):
+ for seqname in input_yaml['unitas']['refseq'][i]:
+ input_yaml['unitas']['refseq'][i][seqname] = make_path_reletive_to_config(input_yaml['unitas']['refseq'][i][seqname], input_yaml_path)
+
+ # targetid
+ if 'targetid' in input_yaml:
+ if 'sort' in input_yaml:
+ input_yaml['targetid']['inputs'] = os.path.join(output_dir, 'binned_length_rna')
+ elif 'size_sorted_fastqs' in input_yaml:
+ input_yaml['targetid']['inputs'] = input_yaml['size_sorted_fastqs']
+ elif 'trim' in input_yaml:
+ input_yaml['targetid']['inputs'] = os.path.join(output_dir, 'trimmed_reads.fq')
+ else:
+ input_yaml['targetid']['inputs'] = input_yaml['smallRNA_fastq']
+
+ # at least one target file specified
+ if 'target_files' not in input_yaml['targetid'] or len(input_yaml['targetid']['target_files']) < 1:
+ return False, 'To use the targetid command, you must specify at least one target file, to reuse cds and unsplices transcriptiome, specify gene', input_yaml
+
+ if 'gene' in input_yaml['targetid']['target_files']:
+ gene_index = input_yaml['targetid']['target_files'].index('gene')
+ del input_yaml['targetid']['target_files'][gene_index]
+
+ if input_yaml['cds'] is None and input_yaml['unspliced_transcriptome'] is None:
+ return False, 'If gene is specified in unitas or targetid refseq or unspliced_transcriptome must also be specified', input_yaml
+
+ if input_yaml['cds'] is not None:
+ input_yaml['targetid']['target_files'].append(input_yaml['cds'])
+
+ if input_yaml['unspliced_transcriptome'] is not None:
+ input_yaml['targetid']['target_files'].append(input_yaml['unspliced_transcriptome'])
+
+ for i in range(len(input_yaml['targetid']['target_files'])):
+ input_yaml['targetid']['target_files'][i] = make_path_reletive_to_config(input_yaml['targetid']['target_files'][i], input_yaml_path)
+
+ # min_seq_length defaults to 5
+ if 'min_seq_length' not in input_yaml['targetid']:
+ input_yaml['targetid']['min_seq_length'] = 5
+
+ # mismatches defaults to 0
+ if 'mismatches' not in input_yaml['targetid']:
+ input_yaml['targetid']['mismatches'] = 0
+
+ # optional functional enrichment part
+ if 'enrich' in input_yaml['targetid']:
+ if 'eggnog_data' not in input_yaml['targetid']['enrich']:
+ return False, 'If enrich is specified, at least eggnog_data must be specified in it, either add eggnog_data or remove enrich', input_yaml
+
+ # Create a list of incuded files based on target_files and exclude
+ if 'exclude' in input_yaml['targetid']['enrich']:
+ if type(input_yaml['targetid']['enrich']['exclude']) == str:
+ input_yaml['targetid']['enrich']['exclude'] = [input_yaml['targetid']['enrich']['exclude']]
+ elif type(input_yaml['targetid']['enrich']['exclude']) != list:
+ return False, f'In enrich, the exclude part should be a list of paths, not {type(input_yaml["targetid"]["enrich"]["exclude"])}', input_yaml
+
+ for i in range(len(input_yaml['targetid']['enrich']['exclude'])):
+ input_yaml['targetid']['enrich']['exclude'][i] = make_path_reletive_to_config(input_yaml['targetid']['enrich']['exclude'][i], input_yaml_path)
+
+ input_yaml['targetid']['enrich']['included_files'] = set(input_yaml['targetid']['target_files']) - set(input_yaml['targetid']['enrich']['exclude'])
+ else:
+ input_yaml['targetid']['enrich']['included_files'] = set(input_yaml['targetid']['target_files'])
+
+ return True, None, input_yaml
\ No newline at end of file
diff --git a/mkdocs.yml b/mkdocs.yml
new file mode 100644
index 0000000..937549a
--- /dev/null
+++ b/mkdocs.yml
@@ -0,0 +1,41 @@
+site_name: Hunt Lab Small RNA Pipeline
+site_url: https://Vicky-Hunt-Lab.github.io/HuntLab-smallRNA/
+repo_url: https://github.com/Vicky-Hunt-Lab/HuntLab-smallRNA
+repo_name: Vicky-Hunt-Lab/HuntLab-smallRNA
+copyright: Copyright ©