README.md

Combine multiple compressed FASTQ files, converting them to FASTA, using pass.

st pass file1.fastq.gz file2.fastq.gz -o output.fasta

Note: almost every command can read multiple input files and convert between formats, but pass does nothing other than reading and writing while other command perform certain actions.

Aside from ID and sequence, any variable/function such as the sequence length (seqlen) can be written to delimited text.

st pass input.fasta --to-tsv id,seq,seqlen

id1	ACG	3
id1	ACGTACGT	7
id1	ACGTA	5

st count -k path data/*.fastq.gz

data/sample1.fastq.gz	30601
data/sample2.fastq.gz	15702
data/sample3.fastq.gz	264965
data/sample4.fastq.gz	1120
data/sample5.fastq.gz	7021
(...)

In count, one or several categorical variables/functions can be specified with -k/--key.

The function bin(variable, interval) groups continuous numeric values into intervals

st count -k 'bin(gc_percent, 10)' sequences.fasta

(10, 20]	57
(20, 30]	2113
(30, 40]	11076
(40, 50]	7184
(50, 60]	12

st set -i 'seq_{num}' seqs.fasta > renamed.fasta

>seq_1
SEQUENCE
>seq_2
SEQUENCE
>seq_3
SEQUENCE
(...)

seqs.fasta with a 'group' annotation in the header:

>id1 group1
SEQUENCE1
>id2 group1
SEQUENCE2
>id3 group1
SEQUENCE2
>id4 group2
SEQUENCE1
>id5 group2
SEQUENCE1

st unique 'desc,seq' seqs.fasta > grouped_uniques.fasta

>id1 group1
SEQUENCE1
>id2 group1
SEQUENCE2
>id4 group2
SEQUENCE1

This filter command removes sequencing reads with more than one expected sequencing error (like USEARCH can do) or sequence length of <100 bp.

st filter 'exp_err < 1 && seqlen >= 100' reads.fastq > filtered.fastq

The unique command returns all unique sequences and annotates the number of records with the same sequence in the header:

st unique seq -a abund={n_duplicates} input.fasta > uniques.fasta

>id1 abund=3
TCTTTAATAACCTGATTAG
>id3 abund=1
GGAGGATCCGAGCG
(...)

It is also possible to de-replicate by multiple keys, e.g. by sequence, but grouped by a sample attribute in the header:

st unique 'seq,attr(sample)' input.fasta > uniques.fasta

>id1 sample=1
SEQUENCE1
>id3 sample=2
SEQUENCE2
>id10 sample=1
SEQUENCE3
>id11 sample=3
SEQUENCE4
(...)

These steps could be part of an amplicon pipeline that de-multiplexes multi-marker amplicons.

Three different seqtool commands are chained, whereby find searches for a set of primers, which are removed by trim and finally split distributes the sequences into different files named by the forward primer. The primer names and positions are passed to the next command in the form of header attributes (-a/--attr option).

primers.fasta

>prA
PRIMER
>prB
PRIMER

Command for searching/trimming

st find file:primers.fasta -a primer='{pattern_name}' -a end='{match_end}' sequences.fasta |
  st trim -e '{attr(end)}..' | 
  st split -o '{attr(primer)}'

>id1 primer=prA end=22
SEQUENCE
>id4 primer=prA end=21
SEQUENCE
(...)

>id2 primer=prB end=20
SEQUENCE
>id3 primer=prB end=22
SEQUENCE
(...)

>id5 primer=undefined end=undefined
UNTRIMMEDSEQUENCE
(...)

>id1
SEQUENCE
>id2
SEQUENCE
(...)

id  genus
seq1  Actinomyces
seq2  Amycolatopsis
(...)

Using -m/--meta to include genus.tsv as metadata source:

st set -m genus.tsv --desc '{meta(genus)}' input.fasta > with_genus.fasta

>seq1 Actinomyces
SEQUENCE
>seq2 Amycolatopsis
SEQUENCE
(...)

>id1
SEQUENCE
>id2
SEQUENCE
>id3
SEQUENCE
>id4
SEQUENCE

id1
id4

st filter -m id_list.txt 'has_meta()' input.fasta > subset.fasta

>id1
SEQUENCE
>id4
SEQUENCE