BCH709 Introduction to Bioinformatics: Seqkit Tutorial

Usage and Examples

Quick Guide

Technical details and guides for use

FASTA/Q format parsing and writing

Seqkit also supports reading and writing xz (.xz) and zstd (.zst) formats since v2.2.0. Bzip2 format is supported since v2.4.0.

SeqKit uses the author’s lightweight and high-performance bioinformatics package bio for FASTA/Q parsing, which has high performance close to the famous C lib klib (kseq.h).

Notes:

Input and output files

Seqkit accepts input data from standard input (STDIN) and plain or gzip-compressed files. Files can be given via positional arguments or the flag --infile-list. For example:

seqkit seq   a.fasta b.fasta

seqkit seq --infile-list file-list.txt

seqkit seq --infile-list <(find -name "*.fq.gz" directory)

Result are printed to standard ouput (STDOUT) by default, you can also specify the output file via the flag -o/--out-file. The file name extension .gz is automatically recognized. For example:

seqkit seq a.fasta -o a.fasta.gz

seqkit grep -f IDs.txt read_1.fq.gz -o dir/read_1.fq.gz

Seqkit utlizies the pgzip package to read and write gzip file, and the outputted gzip file would be slighty larger than files generated by GNU gzip.

Seqkit writes gzip files very fast, much faster than the multi-threaded pigz, so there’s no need to pipe the result to gzip/pigz.

Sequence formats and types

SeqKit seamlessly support FASTA and FASTQ format. Sequence format is automatically detected. All subcommands except for faidx and bam can handle both formats. And only when some commands (subseq, split, sort and shuffle) which utilise FASTA index to improve perfrmance for large files in two pass mode (by flag --two-pass), only FASTA format is supported.

Sequence type (DNA/RNA/Protein) is automatically detected by leading subsequences of the first sequences in file or STDIN. The length of the leading subsequences is configurable by global flag --alphabet-guess-seq-length with default value of 10000. If length of the sequences is less than that, whole sequences will be checked.

Sequence ID

By default, most softwares, including seqkit, take the leading non-space letters as sequence identifier (ID). For example,

FASTA header ID
>123456 gene name 123456
>longname longname
>gi|110645304|ref|NC_002516.2| Pseudomona gi|110645304|ref|NC_002516.2|

But for some sequences from NCBI, e.g. >gi|110645304|ref|NC_002516.2| Pseudomona, the ID is NC_002516.2. In this case, we could set sequence ID parsing regular expression by global flag --id-regexp "\|([^\|]+)\| " or just use flag --id-ncbi. If you want the gi number, then use --id-regexp "^gi\|([^\|]+)\|".

FASTA index

For some commands, including subseq, split, sort and shuffle, when input files are (plain or gzipped) FASTA files, FASTA index would be optional used for rapid access of sequences and reducing memory occupation.

ATTENTION: the .seqkit.fai file created by SeqKit is a little different from .fai file created by samtools. SeqKit uses full sequence head instead of just ID as key.

Parallelization of CPU intensive jobs

The validation of sequences bases and complement process of sequences are parallelized for large sequences.

Parsing of line-based files, including BED/GFF file and ID list file are also parallelized.

The pgzip package reads and write gzip files in parallel.

The Parallelization is implemented by multiple goroutines in golang which are similar to but much lighter weight than threads. The concurrency number is configurable with global flag -j or --threads (default value: 4).

Using four threads is fast enough for most commands where FASTA/Q reading and writing is the performance bottleneck, and using more threads will not increase the speed.

Few commands could benefit from multiple (>4) threads:

Memory occupation

Most of the subcommands do not read whole FASTA/Q records in to memory, including stat, fq2fa, fx2tab, tab2fx, grep, locate, replace, seq, sliding, subseq.

Note that when using subseq --gtf | --bed, if the GTF/BED files are too big, the memory usage will increase. You could use --chr to specify chromesomes and --feature to limit features.

Some subcommands could either read all records or read the files twice by flag -2 (--two-pass), including sample, split, shuffle and sort. They use FASTA index for rapid acccess of sequences and reducing memory occupation.

Reproducibility

Subcommands sample and shuffle use random function, random seed could be given by flag -s (--rand-seed). This makes sure that sampling result could be reproduced in different environments with same random seed.

seqkit

SeqKit -- a cross-platform and ultrafast toolkit for FASTA/Q file manipulation

Version: 2.5.0

Author: Wei Shen <shenwei356@gmail.com>

Documents  : http://bioinf.shenwei.me/seqkit
Source code: https://github.com/shenwei356/seqkit
Please cite: https://doi.org/10.1371/journal.pone.0163962


Seqkit utlizies the pgzip (https://github.com/klauspost/pgzip) package to
read and write gzip file, and the outputted gzip file would be slighty
larger than files generated by GNU gzip.

Seqkit writes gzip files very fast, much faster than the multi-threaded pigz,
therefore there's no need to pipe the result to gzip/pigz.

Seqkit also supports reading and writing xz (.xz) and zstd (.zst) formats since v2.2.0.
Bzip2 format is supported since v2.4.0.

Compression level:
  format   range   default  comment
  gzip     1-9     5        https://github.com/klauspost/pgzip sets 5 as the default value.
  xz       NA      NA       https://github.com/ulikunitz/xz does not support.
  zstd     1-4     2        roughly equals to zstd 1, 3, 7, 11, respectively.
  bzip     1-9     6        https://github.com/dsnet/compress

Usage:
  seqkit [command]

Available Commands:
  amplicon        extract amplicon (or specific region around it) via primer(s)
  bam             monitoring and online histograms of BAM record features
  common          find common sequences of multiple files by id/name/sequence
  concat          concatenate sequences with the same ID from multiple files
  convert         convert FASTQ quality encoding between Sanger, Solexa and Illumina
  duplicate       duplicate sequences N times
  fa2fq           retrieve corresponding FASTQ records by a FASTA file
  faidx           create FASTA index file and extract subsequence
  fish            look for short sequences in larger sequences using local alignment
  fq2fa           convert FASTQ to FASTA
  fx2tab          convert FASTA/Q to tabular format (and length, GC content, average quality...)
  genautocomplete generate shell autocompletion script (bash|zsh|fish|powershell)
  grep            search sequences by ID/name/sequence/sequence motifs, mismatch allowed
  head            print first N FASTA/Q records
  head-genome     print sequences of the first genome with common prefixes in name
  locate          locate subsequences/motifs, mismatch allowed
  merge-slides    merge sliding windows generated from seqkit sliding
  mutate          edit sequence (point mutation, insertion, deletion)
  pair            match up paired-end reads from two fastq files
  range           print FASTA/Q records in a range (start:end)
  rename          rename duplicated IDs
  replace         replace name/sequence by regular expression
  restart         reset start position for circular genome
  rmdup           remove duplicated sequences by ID/name/sequence
  sample          sample sequences by number or proportion
  sana            sanitize broken single line FASTQ files
  scat            real time recursive concatenation and streaming of fastx files
  seq             transform sequences (extract ID, filter by length, remove gaps, reverse complement...)
  shuffle         shuffle sequences
  sliding         extract subsequences in sliding windows
  sort            sort sequences by id/name/sequence/length
  split           split sequences into files by id/seq region/size/parts (mainly for FASTA)
  split2          split sequences into files by size/parts (FASTA, PE/SE FASTQ)
  stats           simple statistics of FASTA/Q files
  subseq          get subsequences by region/gtf/bed, including flanking sequences
  sum             compute message digest for all sequences in FASTA/Q files
  tab2fx          convert tabular format to FASTA/Q format
  translate       translate DNA/RNA to protein sequence (supporting ambiguous bases)
  version         print version information and check for update
  watch           monitoring and online histograms of sequence features

Flags:
      --alphabet-guess-seq-length int   length of sequence prefix of the first FASTA record based on
                                        which seqkit guesses the sequence type (0 for whole seq)
                                        (default 10000)
      --compress-level int              compression level for gzip, zstd, xz and bzip2. type "seqkit -h"
                                        for the range and default value for each format (default -1)
  -h, --help                            help for seqkit
      --id-ncbi                         FASTA head is NCBI-style, e.g. >gi|110645304|ref|NC_002516.2|
                                        Pseud...
      --id-regexp string                regular expression for parsing ID (default "^(\\S+)\\s?")
      --infile-list string              file of input files list (one file per line), if given, they are
                                        appended to files from cli arguments
  -w, --line-width int                  line width when outputting FASTA format (0 for no wrap) (default 60)
  -o, --out-file string                 out file ("-" for stdout, suffix .gz for gzipped out) (default "-")
      --quiet                           be quiet and do not show extra information
  -t, --seq-type string                 sequence type (dna|rna|protein|unlimit|auto) (for auto, it
                                        automatically detect by the first sequence) (default "auto")
  -j, --threads int                     number of CPUs. can also set with environment variable
                                        SEQKIT_THREADS) (default 4)

Use "seqkit [command] --help" for more information about a command.

Tutorial

SeqKit

Create and Installation

Create a new Conda environment named seqkit and install seqkit along with other bioinformatics tools:

conda create -n seqkit -c bioconda seqkit bedops csvtk clustalo

Activate Environment

Activate the seqkit environment:

conda activate seqkit

Export Environment

Export the current environment configuration to seqkit.yml:

conda env export > seqkit.yml

Deactivate Environment

Deactivate the current Conda environment:

conda deactivate

Remove Conda Environment

Remove the seqkit Conda environment entirely:

conda remove --name seqkit --all

Reinstall Environment

Recreate the seqkit environment from the exported seqkit.yml file:

conda env create -f seqkit.yml

Activate Environment

Activate the seqkit environment again:

conda activate seqkit

Use SeqKit

Display the help information for seqkit:

seqkit --help

Make Homework Folder

Set up a dedicated folder for your SeqKit projects:

cd ~/${USER}
mkdir -p ~/bch709/seqkit
cd ~/bch709/seqkit
pwd
seqkit

Datasets

Test FASTQ File

Download the test FASTQ file:

wget https://raw.githubusercontent.com/plantgenomicslab/BCH709/gh-pages/data/test.fastq.gz
wget https://raw.githubusercontent.com/plantgenomicslab/BCH709/refs/heads/gh-pages/data/contig.fa
wget https://raw.githubusercontent.com/plantgenomicslab/BCH709/refs/heads/gh-pages/data/seqs.fa
wget https://raw.githubusercontent.com/plantgenomicslab/BCH709/refs/heads/gh-pages/data/rna.fa
wget https://raw.githubusercontent.com/plantgenomicslab/BCH709/refs/heads/gh-pages/data/file1.fa
wget https://raw.githubusercontent.com/plantgenomicslab/BCH709/refs/heads/gh-pages/data/file2.fa

wget https://mirbase.org/download/hairpin.fa

Human Genome from UCSC

Download the human genome and related files:

wget --no-check-certificate https://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/latest/hg38.fa.gz -O hsa.fa.gz
wget --no-check-certificate https://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/genes/hg38.refGene.gtf.gz -O hsa.gtf.gz
ls -algh

Examples

1. Read and Print Sequences

2. Only Print Names

3. Only Print Sequences

Define the output line width with the -w flag (use 0 for no wrap):

seqkit seq hairpin.fa -s -w 0

4. Reverse Complement Sequence

seqkit seq contig.fa -r -p

5. Remove Gaps and Convert Case

seqkit seq -g -u seqs.fa

6. RNA to DNA Conversion

seqkit seq --rna2dna rna.fa

Subseq

Usage

Retrieve subsequences by region, GTF, or BED files, including flanking sequences.

Recommendation: Use plain FASTA files to utilize SeqKit’s FASTA index for efficiency.

Region Definition: 1-based indexing with custom design.

Examples:

1-based index    1 2 3 4 5 6 7 8 9 10
negative index    0-9-8-7-6-5-4-3-2-1
           seq    A C G T N a c g t n
           1:1    A
           2:4      C G T
         -4:-2                c g t
         -4:-1                c g t n
         -1:-1                      n
          2:-2      C G T N a c g t
          1:-1    A C G T N a c g t n

Usage Command

seqkit subseq [flags]

Flags

Examples

Recommendation: Use plain FASTA files to utilize SeqKit’s FASTA index.

  1. First 12 Bases:

     cat hairpin.fa | seqkit subseq -r 1:12

  2. Last 12 Bases:

     cat hairpin.fa | seqkit subseq -r -12:-1

  3. Excluding First and Last 12 Bases:

     cat hairpin.fa | seqkit subseq -r 13:-13

  4. Get Subsequence by GTF File:

    #### Human Genome Example

    Specify chromosomes and features:

     seqkit subseq --gtf hsa.gtf.gz --chr chr1 --feature CDS hsa.fa.gz > chr1.gtf.cds.fa

     seqkit stat chr1.gtf.cds.fa

  5. Remove Duplicated Sequences:

     seqkit subseq --gtf hsa.gtf.gz --chr chr1 --feature CDS hsa.fa.gz | seqkit rmdup > chr1.gtf.rmdup.fa

     seqkit stat chr1.gtf.rmdup.fa

Human Genome Statistics

seqkit stat hsa.fa.gz

file    format  type  num_seqs        sum_len  min_len       avg_len      max_len
hsa.fa  FASTA   DNA        194  3,099,750,718      970  15,978,096.5  248,956,422

Building FASTA Index:

seqkit faidx hsa.fa.gz -o hsa.fa.fai2

Sorting by Sequence Length:

seqkit sort --by-length --reverse --two-pass hsa.fa.gz > hsa.sorted.fa

[INFO] create and read FASTA index ...
[INFO] read sequence IDs and lengths from FASTA index ...
[INFO] 194 sequences loaded
[INFO] sorting ...
[INFO] output ...

seqkit fx2tab --length hsa.sorted.fa --name --only-id | cut -f 1,4 | more

Shuffling Sequences:

seqkit shuffle hsa.fa.gz --two-pass > hsa.shuffled.fa

[INFO] create and read FASTA index ...
[INFO] read sequence IDs from FASTA index ...
[INFO] 194 sequences loaded
[INFO] shuffle ...
[INFO] output ...

Splitting into Files with Single Sequence:

seqkit split --by-id hsa.fa.gz --two-pass

[INFO] split by ID. idRegexp: ^([^\s]+)\s?
[INFO] create and read FASTA index ...
[INFO] read sequence IDs from FASTA index ...
[INFO] 194 sequences loaded
[INFO] write 1 sequences to file: hsa.id_KI270743.1.fa
[INFO] write 1 sequences to file: hsa.id_KI270706.1.fa
[INFO] write 1 sequences to file: hsa.id_KI270717.1.fa
[INFO] write 1 sequences to file: hsa.id_KI270718.1.fa
[INFO] write 1 sequences to file: hsa.id_KI270468.1.fa
...

Sliding

Description

Slide sequences with support for circular genomes.

Usage Command

seqkit sliding [flags]

Flags

Examples

  1. General Use:

     echo -e ">seq\nACGTacgtNN" | seqkit sliding -s 3 -W 6

  2. Circular Genome:

     echo -e ">seq\nACGTacgtNN" | seqkit sliding -s 3 -W 6 -C

  3. Generate GC Content for Plotting:

     cat hairpin.fa | seqkit fx2tab | head -n 1 | seqkit tab2fx | seqkit sliding -s 5 -W 30 | seqkit fx2tab -n -g

Stat

Description

Provide simple statistics of FASTA/Q files.

Usage Command

seqkit stat [flags]

Examples

  1. General Use:

     seqkit stat *.f*{a,q}.gz

Fq2fa

Description

Convert FASTQ files to FASTA format.

Usage Command

seqkit fq2fa [flags]

Example

seqkit fq2fa test.fastq.gz -o test_.fa.gz
zcat test_.fa.gz
zcat test.fastq.gz

Fx2tab & Tab2fx

Fx2tab

Description: Convert FASTA/Q to tabular format, providing various information like sequence length, GC content, and GC skew.

Usage Command

seqkit fx2tab [flags]

Flags

Tab2fx

Description: Convert tabular format (first two/three columns) back to FASTA/Q format.

Usage Command

seqkit tab2fx [flags]

Flags

Examples

  1. Default Output:

     seqkit fx2tab hairpin.fa | head -n 2

  2. Print Sequence Length, GC Content, and Only Names:

     seqkit fx2tab hairpin.fa -l -g -n -i -H | head -n 4 | csvtk -t -C '&' pretty

  3. Use Fx2tab and Tab2fx in a Pipeline:

     cat hairpin.fa | seqkit fx2tab | seqkit tab2fx
 zcat test.fastq.gz | seqkit fx2tab | seqkit tab2fx

  4. Sort Sequences by Length:

     cat hairpin.fa | seqkit fx2tab -l | sort -t"$(echo -e '\t')" -n -k4,4 | seqkit tab2fx

    Alternatively:

     seqkit sort -l hairpin.fa

  5. Get First 1000 Sequences:

     seqkit fx2tab hairpin.fa | head -n 1000 | seqkit tab2fx
 seqkit fx2tab test.fastq.gz | head -n 1000 | seqkit tab2fx

Extension

After converting FASTA to tabular format with seqkit fx2tab, you can handle the data with CSV/TSV tools such as csvtk, a cross-platform, efficient, and practical CSV/TSV toolkit.

Grep

Description

Search sequences by patterns of names or sequence motifs.

Usage Command

seqkit grep [flags]

Flags

Examples

  1. Extract Human Hairpins (Names Starting with hsa):

     cat hairpin.fa | seqkit grep -r -p ^hsa

  2. Remove Human and Mice Hairpins:

     cat hairpin.fa | seqkit grep -r -p ^hsa -p ^mmu -v

  3. Extract New Entries from miRNA.diff:

    a. Get IDs of New Entries:

     cat miRNA.diff | grep -v ^# | grep NEW | cut -f 2 > list

    b. Extract by ID List File:

     cat hairpin.fa | seqkit grep -f list > new.fa

  4. Extract Sequences Starting with AGGCG:

     cat hairpin.fa | seqkit grep -s -r -i -p ^aggcg

  5. Extract Sequences with TTSAA (AgsI Digest Site) in Sequence:

     cat hairpin.fa | seqkit grep -s -d -i -p UUSAA

    Equivalent to:

     cat hairpin.fa | seqkit grep -s -r -i -p UU[CG]AA

Locate

Description

Locate subsequences or motifs within sequences.

Motif Types:

Note: When using the -d flag, regular expressions may not function as expected.

Usage Command

seqkit locate [flags]

Flags

Examples

  1. Locate ORFs:

     cat hairpin.fa | seqkit locate -i -r -p "A[TU]G(?:.{3})+?[TU](?:AG|AA|GA)"

  2. Locate Motif:

     cat hairpin.fa | seqkit locate -i -d -p AUGGACUN

    Note: seqkit grep only searches on the positive strand, but seqkit locate can recognize both strands.

Rmdup

Description

Remove duplicated sequences by ID, name, or sequence.

Usage Command

seqkit rmdup [flags]

Flags

Examples

  1. General Use:

     cat hairpin.fa | seqkit rmdup -s -o clean.fa.gz
 zcat test.fastq.gz | seqkit rmdup -s -o clean.fa.gz

  2. Save Duplicated Sequences to File:

     cat hairpin.fa | seqkit rmdup -s -i -o clean.fa.gz -d duplicated.fa.gz -D duplicated.detail.txt
 cat duplicated.detail.txt

Common

Description

Find common sequences across multiple files by ID, name, or sequence.

Usage Command

seqkit common [flags]

Flags

Examples

  1. Check Files:

     cat file1.fa
 cat file2.fa

  2. By ID (Default):

     seqkit common file*.fa -o common.fasta

  3. By Full Name:

     seqkit common file*.fa -n -o common.fasta

  4. By Sequence:

     seqkit common file*.fa -s -i -o common.fasta

Split

Description

Split sequences into separate files based on:

Region Definition: 1-based indexing with custom design.

Examples:

1-based index    1 2 3 4 5 6 7 8 9 10
negative index    0-9-8-7-6-5-4-3-2-1
           seq    A C G T N a c g t n
           1:1    A
           2:4      C G T
         -4:-2                c g t
         -4:-1                c g t n
         -1:-1                      n
          2:-2      C G T N a c g t
          1:-1    A C G T N a c g t n

Usage Command

seqkit split [flags]

Flags

Examples

  1. Split Sequences into Parts with At Most 10,000 Sequences:

     seqkit split hairpin.fa -s 10000

  2. Split Sequences into 4 Parts:

     seqkit split hairpin.fa -p 4

    To reduce memory usage when splitting large files, always use the --two-pass flag:

     seqkit split hairpin.fa -p 4 -2

  3. Split Sequences by Species (Custom IDs - First Three Letters):

     seqkit split hairpin.fa -i --id-regexp "^([\w]+)\-" --two-pass

  4. Split Sequences by Sequence Region (e.g., Sequence Barcode):

     seqkit split hairpin.fa -r 1:3 -2

    Example Output:

     [INFO] split by region: 1:3
 [INFO] read and write sequences to temporary file: hairpin.fa.gz.fa ...
 [INFO] read sequence IDs and sequence region from FASTA file ...
 [INFO] create and read FASTA index ...
 [INFO] write 463 sequences to file: hairpin.region_1:3_AUG.fa.gz
 [INFO] write 349 sequences to file: hairpin.region_1:3_ACU.fa.gz
 [INFO] write 311 sequences to file: hairpin.region_1:3_CGG.fa.gz

    Sequence suffix can be defined as -r -12:-1.

Sample

Description

Sample sequences by number or proportion.

Usage Command

seqkit sample [flags]

Flags

Examples

  1. Sample by Proportion:

     cat hairpin.fa | seqkit sample -p 0.1 -o sample.fa.gz

  2. Sample by Number:

     cat hairpin.fa | seqkit sample -n 1000 -o sample.fa.gz

    To reduce memory usage when sampling large files, use the --two-pass flag:

     cat hairpin.fa | seqkit sample -n 1000 --two-pass -o sample.fa.gz

  3. Sample and Then Use head to Limit to 1,000 Sequences:

     cat hairpin.fa | seqkit sample -p 0.1 | seqkit head -n 1000 -o sample.fa.gz

  4. Set Random Seed to Reproduce Results:

     cat hairpin.fa | seqkit sample -p 0.1 -s 11

  5. Shuffle After Sampling:

     cat hairpin.fa | seqkit sample -p 0.1 | seqkit shuffle -o sample.fa.gz

    Note: When sampling FASTQ files, ensure using the same random seed with the -s (--rand-seed) flag.

Description

Print the first N FASTA/Q records from a file.

Usage Command

seqkit head [flags]

Flags

Examples

  1. FASTA:

     seqkit head -n 1 hairpin.fa

  2. FASTQ:

     seqkit head -n 1 test.fastq.gz

Replace

Description

Replace names or sequences using regular expressions.

Note: Replacement supports capture variables (e.g., $1 for the first submatch). Use single quotes ' instead of double quotes " in Unix-like systems.

Usage Command

seqkit replace [flags]

Flags

Examples

  1. Remove Descriptions:

     echo -e ">seq1 abc-123\nACGT-ACGT" | seqkit replace -p " .+"

  2. Replace “-“ with “=”:

     echo -e ">seq1 abc-123\nACGT-ACGT" | seqkit replace -p "\-" -r '='

  3. Remove Gaps in Sequences:

     echo -e ">seq1 abc-123\nACGT-ACGT" | seqkit replace -p " |-" -s

  4. Add Space to Every Base:

     echo -e ">seq1 abc-123\nACGT-ACGT" | seqkit replace -p "(.)" -r '$1 ' -s

  5. Transpose Sequence with csvtk:

     echo -e ">seq1\nACTGACGT\n>seq2\nactgccgt" | \
 seqkit replace -p "(.)" -r "\$1 " -s | \
 seqkit seq -s -u | \
 csvtk space2tab | \
 csvtk -t transpose

  6. Rename with Record Number:

     echo -e ">abc\nACTG\n>123\nATTT" | seqkit replace -p .+ -r "seq_{NR}"

More Information

Refer to the SeqKit Replace Usage for detailed information.

Shuffle

Description

Randomly shuffle sequences within a file.

Note: By default, all records are loaded into memory. For large FASTA files, use the --two-pass flag to reduce memory usage. FASTQ files are not supported.

Usage Command

seqkit shuffle [flags]

Flags

Examples

  1. General Use:

     seqkit shuffle hairpin.fa > shuffled.fa

  2. Shuffle Large Genome with Two-Pass Mode:

     time seqkit shuffle -2 GRCh38_latest_genomic.fna.gz > shuffle.fa

Sort

Description

Sort sequences by ID, name, sequence, or length.

Note: By default, all records are loaded into memory. For large FASTA files, use the --two-pass flag to reduce memory usage. FASTQ files are not supported.

Usage Command

seqkit sort [flags]

Flags

Examples

Note: For FASTA format, always use the --two-pass flag to minimize memory usage.

  1. Sort by ID:

     echo -e ">seq1\nACGTNcccc\n>SEQ2\nacgtnAAAA" | seqkit sort --quiet

  2. Sort by ID, Ignoring Case:

     echo -e ">seq1\nACGTNcccc\n>SEQ2\nacgtnAAAA" | seqkit sort --quiet -i

  3. Sort by Sequence, Ignoring Case:

     echo -e ">seq1\nACGTNcccc\n>SEQ2\nacgtnAAAA" | seqkit sort --quiet -s -i

  4. Sort by Sequence Length:

     echo -e ">seq1\nACGTNcccc\n>SEQ2\nacgtnAAAAnnn\n>seq3\nacgt" | seqkit sort --quiet -l

Quick Glance

  1. Sequence Number:

     seqkit stat hairpin.fa

  2. First 10 Bases:

     cat hairpin.fa | seqkit subseq -r 1:10 | seqkit sort -s | seqkit seq -s | head -n 10

Repeated Hairpin Sequences

Check for identical hairpins across different species and remove duplicates:

seqkit rmdup -s -i hairpin.fa -o clean.fa.gz -d duplicated.fa.gz -D duplicated.detail.txt

Results:

Hairpins in Different Species

  1. Inspect Sequence Names:

     seqkit seq hairpin.fa.gz -n | head -n 3

    Example Output:

     cel-let-7 MI0000001 Caenorhabditis elegans let-7 stem-loop
 cel-lin-4 MI0000002 Caenorhabditis elegans lin-4 stem-loop
 cel-mir-1 MI0000003 Caenorhabditis elegans miR-1 stem-loop

    Observation: The first three letters (e.g., cel) abbreviate species names.

  2. Split Sequences by Species:

    Use a custom regular expression to split based on species abbreviation:

     seqkit split hairpin.fa -i --id-regexp "^([\w]+)\-" --two-pass

    Note: Always use the --two-pass flag to handle large files efficiently.

  3. Identify Species with Most miRNA Hairpins:

     cd hairpin.fa.split/
 seqkit stat hairpin.part_* | csvtk space2tab | csvtk -t sort -k num_seqs:nr | csvtk -t pretty | more

Additional Notes

Other Datasets

Datasets from The miRBase Sequence Database – Release 21

Human genome from ensembl (For seqkit subseq)

Only DNA and gtf/bed data of Chr1 were used:

seq

Usage

transform sequences (extract ID, filter by length, remove gaps, reverse complement...)

Usage:
  seqkit seq [flags] 

Flags:
  -k, --color                     colorize sequences - to be piped into "less -R"
  -p, --complement                complement sequence, flag '-v' is recommended to switch on
      --dna2rna                   DNA to RNA
  -G, --gap-letters string        gap letters to be removed with -g/--remove-gaps (default "- \t.")
  -h, --help                      help for seq
  -l, --lower-case                print sequences in lower case
  -M, --max-len int               only print sequences shorter than or equal to the maximum length (-1
                                  for no limit) (default -1)
  -R, --max-qual float            only print sequences with average quality less than this limit (-1 for
                                  no limit) (default -1)
  -m, --min-len int               only print sequences longer than or equal to the minimum length (-1
                                  for no limit) (default -1)
  -Q, --min-qual float            only print sequences with average quality qreater or equal than this
                                  limit (-1 for no limit) (default -1)
  -n, --name                      only print names/sequence headers
  -i, --only-id                   print IDs instead of full headers
  -q, --qual                      only print qualities
  -b, --qual-ascii-base int       ASCII BASE, 33 for Phred+33 (default 33)
  -g, --remove-gaps               remove gaps letters set by -G/--gap-letters, e.g., spaces, tabs, and
                                  dashes (gaps "-" in aligned sequences)
  -r, --reverse                   reverse sequence
      --rna2dna                   RNA to DNA
  -s, --seq                       only print sequences
  -u, --upper-case                print sequences in upper case
  -v, --validate-seq              validate bases according to the alphabet
  -V, --validate-seq-length int   length of sequence to validate (0 for whole seq) (default 10000)

Examples

  1. Read and print

    • From file:

        $ seqkit seq hairpin.fa.gz
  >cel-let-7 MI0000001 Caenorhabditis elegans let-7 stem-loop
  UACACUGUGGAUCCGGUGAGGUAGUAGGUUGUAUAGUUUGGAAUAUUACCACCGGUGAAC
  UAUGCAAUUUUCUACCUUACCGGAGACAGAACUCUUCGA

  $ seqkit seq reads_1.fq.gz
  @HWI-D00523:240:HF3WGBCXX:1:1101:2574:2226 1:N:0:CTGTAG
  TGAGGAATATTGGTCAATGGGCGCGAGCCTGAACCAGCCAAGTAGCGTGAAGGATGACTGCCCTACGGG
  +
  HIHIIIIIHIIHGHHIHHIIIIIIIIIIIIIIIHHIIIIIHHIHIIIIIGIHIIIIHHHHHHGHIHIII

    • From stdin:

        zcat hairpin.fa.gz | seqkit seq

  2. Sequence types

    • By default, seqkit seq automatically detect the sequence type

        $ echo -e ">seq\nacgtryswkmbdhvACGTRYSWKMBDHV" | seqkit stats
  file  format  type  num_seqs  sum_len  min_len  avg_len  max_len
  -     FASTA   DNA          1       28       28       28       28

  $ echo -e ">seq\nACGUN ACGUN" | seqkit stats
  file  format  type  num_seqs  sum_len  min_len  avg_len  max_len
  -     FASTA   RNA          1       11       11       11       11

  $ echo -e ">seq\nabcdefghijklmnpqrstvwyz" | seqkit stats
  file  format  type     num_seqs  sum_len  min_len  avg_len  max_len
  -     FASTA   Protein         1       23       23       23       23

  $ echo -e "@read\nACTGCN\n+\n@IICCG" | seqkit stats
  file  format  type  num_seqs  sum_len  min_len  avg_len  max_len
  -     FASTQ   DNA          1        6        6        6        6

    • You can also set sequence type by flag -t (--seq-type). But this only take effect on subcommands seq and locate.

        $ echo -e ">seq\nabcdefghijklmnpqrstvwyz" | seqkit seq -t dna
  [INFO] when flag -t (--seq-type) given, flag -v (--validate-seq) is automatically switched on
  [ERRO] error when parsing seq: seq (invalid DNAredundant letter: e)

  3. Only print names

    • Full name:

        $ seqkit seq hairpin.fa.gz -n
  cel-let-7 MI0000001 Caenorhabditis elegans let-7 stem-loop
  cel-lin-4 MI0000002 Caenorhabditis elegans lin-4 stem-loop
  cel-mir-1 MI0000003 Caenorhabditis elegans miR-1 stem-loop

    • Only ID:

        $ seqkit seq hairpin.fa.gz -n -i
  cel-let-7
  cel-lin-4
  cel-mir-1

    • Custom ID region by regular expression (this could be applied to all subcommands):

        $ seqkit seq hairpin.fa.gz -n -i --id-regexp "^[^\s]+\s([^\s]+)\s"
  MI0000001
  MI0000002
  MI0000003

  4. Only print seq (global flag -w defines the output line width, 0 for no wrap)

     $ seqkit seq hairpin.fa.gz -s -w 0
 UACACUGUGGAUCCGGUGAGGUAGUAGGUUGUAUAGUUUGGAAUAUUACCACCGGUGAACUAUGCAAUUUUCUACCUUACCGGAGACAGAACUCUUCGA
 AUGCUUCCGGCCUGUUCCCUGAGACCUCAAGUGUGAGUGUACUAUUGAUGCUUCACACCUGGGCUCUCCGGGUACCAGGACGGUUUGAGCAGAU
 AAAGUGACCGUACCGAGCUGCAUACUUCCUUACAUGCCCAUACUAUAUCAUAAAUGGAUAUGGAAUGUAAAGAAGUAUGUAGAACGGGGUGGUAGU

  5. Convert multi-line FASTQ to 4-line FASTQ

     $ seqkit seq reads_1.fq.gz -w 0

  6. Reverse comlement sequence

     $ seqkit seq hairpin.fa.gz -r -p
 >cel-let-7 MI0000001 Caenorhabditis elegans let-7 stem-loop
 UCGAAGAGUUCUGUCUCCGGUAAGGUAGAAAAUUGCAUAGUUCACCGGUGGUAAUAUUCC
 AAACUAUACAACCUACUACCUCACCGGAUCCACAGUGUA

  7. Remove gaps and to lower/upper case

     $ echo -e ">seq\nACGT-ACTGC-ACC" | seqkit seq -g -u
 >seq
 ACGTACTGCACC

  8. RNA to DNA

     $ echo -e ">seq\nUCAUAUGCUUGUCUCAAAGAUUA" | seqkit seq --rna2dna
 >seq
 TCATATGCTTGTCTCAAAGATTA

  9. Filter by sequence length

     $ cat hairpin.fa | seqkit seq | seqkit stats
 file  format  type  num_seqs    sum_len  min_len  avg_len  max_len
 -     FASTA   RNA     28,645  2,949,871       39      103    2,354

 $ cat hairpin.fa | seqkit seq -m 100 | seqkit stats
 file  format  type  num_seqs    sum_len  min_len  avg_len  max_len
 -     FASTA   RNA     10,975  1,565,486      100    142.6    2,354

 $ cat hairpin.fa | seqkit seq -m 100 -M 1000 | seqkit stats
 file  format  type  num_seqs    sum_len  min_len  avg_len  max_len
 -     FASTA   RNA     10,972  1,560,270      100    142.2      938

subseq

Usage

get subsequences by region/gtf/bed, including flanking sequences.

Attentions:
  1. Use "seqkit grep" for extract subsets of sequences.
     "seqtk subseq seqs.fasta id.txt" equals to
     "seqkit grep -f id.txt seqs.fasta"

Recommendation:
  1. Use plain FASTA file, so seqkit could utilize FASTA index.
  2. The flag -U/--update-faidx is recommended to ensure the .fai file matches the FASTA file.

The definition of region is 1-based and with some custom design.

Examples:

 1-based index    1 2 3 4 5 6 7 8 9 10
negative index    0-9-8-7-6-5-4-3-2-1
           seq    A C G T N a c g t n
           1:1    A
           2:4      C G T
         -4:-2                c g t
         -4:-1                c g t n
         -1:-1                      n
          2:-2      C G T N a c g t
          1:-1    A C G T N a c g t n
          1:12    A C G T N a c g t n
        -12:-1    A C G T N a c g t n

Usage:
  seqkit subseq [flags] 

Flags:
      --bed string        by tab-delimited BED file
      --chr strings       select limited sequence with sequence IDs when using --gtf or --bed (multiple
                          value supported, case ignored)
  -d, --down-stream int   down stream length
      --feature strings   select limited feature types (multiple value supported, case ignored, only
                          works with GTF)
      --gtf string        by GTF (version 2.2) file
      --gtf-tag string    output this tag as sequence comment (default "gene_id")
  -h, --help              help for subseq
  -f, --only-flank        only return up/down stream sequence
  -r, --region string     by region. e.g 1:12 for first 12 bases, -12:-1 for last 12 bases, 13:-1 for
                          cutting first 12 bases. type "seqkit subseq -h" for more examples
  -u, --up-stream int     up stream length
  -U, --update-faidx      update the fasta index file if it exists. Use this if you are not sure whether
                          the fasta file changed

Examples

Recommendation: use plain FASTA file, so seqkit could utilize FASTA index.

  1. First 12 bases

     $ zcat hairpin.fa.gz | seqkit subseq -r 1:12

  2. Last 12 bases

     $ zcat hairpin.fa.gz | seqkit subseq -r -12:-1

  3. Subsequences without first and last 12 bases

     $ zcat hairpin.fa.gz | seqkit subseq -r 13:-13

  4. Get subsequence by GTF file

     $ cat t.fa
 >seq
 actgACTGactgn

 $ cat t.gtf
 seq     test    CDS     5       8       .       .       .       gene_id "A"; transcript_id "";
 seq     test    CDS     5       8       .       -       .       gene_id "B"; transcript_id "";

 $ seqkit subseq --gtf t.gtf t.fa
 >seq_5:8:. A
 ACTG
 >seq_5:8:- B
 CAGT

    Human genome example:

    AVOID loading all data from Homo_sapiens.GRCh38.84.gtf.gz, the uncompressed data are so big and may exhaust your RAM.

    We could specify chromesomes and features.

     $ seqkit subseq --gtf Homo_sapiens.GRCh38.84.gtf.gz --chr 1 --feature cds  hsa.fa > chr1.gtf.cds.fa

 $ seqkit stats chr1.gtf.cds.fa
 file             format  type  num_seqs    sum_len  min_len  avg_len  max_len
 chr1.gtf.cds.fa  FASTA   DNA     65,012  9,842,274        1    151.4   12,045

  5. Get CDS and 3bp up-stream sequences

     $ seqkit subseq --gtf t.gtf t.fa -u 3
 >seq_5:8:._us:3 A
 ctgACTG
 >seq_5:8:-_us:3 B
 agtCAGT

  6. Get 3bp up-stream sequences of CDS, not including CDS

     $ seqkit subseq --gtf t.gtf t.fa -u 3 -f
 >seq_5:8:._usf:3 A
 ctg
 >seq_5:8:-_usf:3 B
 agt

  7. Get subsequences by BED file.

    AVOID loading all data from Homo_sapiens.GRCh38.84.gtf.gz, the uncompressed data are so big and may exhaust your RAM.

     $ seqkit subseq --bed Homo_sapiens.GRCh38.84.bed.gz --chr 1 hsa.fa \
     >  chr1.bed.gz.fa

    We may need to remove duplicated sequences

     $ seqkit subseq --bed Homo_sapiens.GRCh38.84.bed.gz --chr 1 hsa.fa \
     | seqkit rmdup > chr1.bed.rmdup.fa
 [INFO] 141060 duplicated records removed

    Summary:

     $ seqkit stats chr1.gz.*.gz
 file               seq_format   seq_type   num_seqs   min_len   avg_len     max_len
 chr1.gz.fa         FASTA        DNA         231,974         1   3,089.5   1,551,957
 chr1.gz.rmdup.fa   FASTA        DNA          90,914         1   6,455.8   1,551,957

sliding

Usage

extract subsequences in sliding windows

Usage:
  seqkit sliding [flags]

Flags:
  -c, --circular          circular genome (same to -C/--circular-genome)
  -C, --circular-genome   circular genome (same to -c/--circular)
  -g, --greedy            greedy mode, i.e., exporting last subsequences even shorter than the windows size
  -h, --help              help for sliding
  -s, --step int          step size
  -S, --sufix string      suffix added to the sequence ID
  -W, --window int        window size

Examples

  1. General use

     $ echo -e ">seq\nACGTacgtNN" | seqkit sliding -s 3 -W 6
 >seq_sliding:1-6
 ACGTac
 >seq_sliding:4-9
 TacgtN

  2. Greedy mode

     $ echo -e ">seq\nACGTacgtNN" | seqkit sliding -s 3 -W 6 -g
 >seq_sliding:1-6
 ACGTac
 >seq_sliding:4-9
 TacgtN
 >seq_sliding:7-12
 gtNN
 >seq_sliding:10-15
 N

  3. Circular genome

     $ echo -e ">seq\nACGTacgtNN" | seqkit sliding -s 3 -W 6 -C
 >seq_sliding:1-6
 ACGTac
 >seq_sliding:4-9
 TacgtN
 >seq_sliding:7-2
 gtNNAC
 >seq_sliding:10-5
 NACGTa

  4. Generate GC content for ploting

     $ zcat hairpin.fa.gz \
     | seqkit sliding -s 5 -W 30 \
     | seqkit fx2tab -n -g
 cel-let-7_sliding:1-30          50.00
 cel-let-7_sliding:6-35          46.67
 cel-let-7_sliding:11-40         43.33
 cel-let-7_sliding:16-45         36.67
 cel-let-7_sliding:21-50         33.33
 cel-let-7_sliding:26-55         40.00
 ...

stats

Usage

simple statistics of FASTA/Q files

Columns:

  1.  file      input file, "-" for STDIN
  2.  format    FASTA or FASTQ
  3.  type      DNA, RNA, Protein or Unlimit
  4.  num_seqs  number of sequences
  5.  sum_len   number of bases or residues       , with gaps or spaces counted
  6.  min_len   minimal sequence length           , with gaps or spaces counted
  7.  avg_len   average sequence length           , with gaps or spaces counted
  8.  max_len   miximal sequence length           , with gaps or spaces counted
  9.  Q1        first quartile of sequence length , with gaps or spaces counted
  10. Q2        median of sequence length         , with gaps or spaces counted
  11. Q3        third quartile of sequence length , with gaps or spaces counted
  12. sum_gap   number of gaps
  13. N50       N50. https://en.wikipedia.org/wiki/N50,_L50,_and_related_statistics#N50
  14. Q20(%)    percentage of bases with the quality score greater than 20
  15. Q30(%)    percentage of bases with the quality score greater than 30
  16. GC(%)     percentage of GC content
  
Attentions:
  1. Sequence length metrics (sum_len, min_len, avg_len, max_len, Q1, Q2, Q3)
     count the number of gaps or spaces. You can remove them with "seqkit seq -g":
         seqkit seq -g input.fasta | seqkit stats

Tips:
  1. For lots of small files (especially on SDD), use big value of '-j' to
     parallelize counting.
  2. Extract one metric with csvtk (https://github.com/shenwei356/csvtk):
         seqkit stats -Ta input.fastq.gz | csvtk cut -t -f "Q30(%)" | csvtk del-header

Usage:
  seqkit stats [flags]

Aliases:
  stats, stat

Flags:
  -N, --N strings            other N50-like stats. value range [0, 100], multiple values supported,
                             e.g., -N 50,90 or -N 50 -N 90
  -a, --all                  all statistics, including quartiles of seq length, sum_gap, N50
  -b, --basename             only output basename of files
  -E, --fq-encoding string   fastq quality encoding. available values: 'sanger', 'solexa',
                             'illumina-1.3+', 'illumina-1.5+', 'illumina-1.8+'. (default "sanger")
  -G, --gap-letters string   gap letters (default "- .")
  -h, --help                 help for stats
  -e, --skip-err             skip error, only show warning message
  -i, --stdin-label string   label for replacing default "-" for stdin (default "-")
  -T, --tabular              output in machine-friendly tabular format

Eexamples

  1. General use

     $ seqkit stats *.f{a,q}.gz
 file           format  type  num_seqs    sum_len  min_len  avg_len  max_len
 hairpin.fa.gz  FASTA   RNA     28,645  2,949,871       39      103    2,354
 mature.fa.gz   FASTA   RNA     35,828    781,222       15     21.8       34
 reads_1.fq.gz  FASTQ   DNA      2,500    567,516      226      227      229
 reads_2.fq.gz  FASTQ   DNA      2,500    560,002      223      224      225

  2. Machine-friendly tabular format

     $ seqkit stats *.f{a,q}.gz -T
 file    format  type    num_seqs        sum_len min_len avg_len max_len
 hairpin.fa.gz   FASTA   RNA     28645   2949871 39      103.0   2354
 mature.fa.gz    FASTA   RNA     35828   781222  15      21.8    34
 Illimina1.8.fq.gz       FASTQ   DNA     10000   1500000 150     150.0   150
 reads_1.fq.gz   FASTQ   DNA     2500    567516  226     227.0   229
 reads_2.fq.gz   FASTQ   DNA     2500    560002  223     224.0   225

 # So you can process the result with tools like csvtk: http://bioinf.shenwei.me/csvtk

 $ seqkit stats *.f{a,q}.gz -T | csvtk pretty -t
 file                format   type   num_seqs   sum_len   min_len   avg_len   max_len
 -----------------   ------   ----   --------   -------   -------   -------   -------
 hairpin.fa.gz       FASTA    RNA    28645      2949871   39        103.0     2354
 mature.fa.gz        FASTA    RNA    35828      781222    15        21.8      34
 Illimina1.8.fq.gz   FASTQ    DNA    10000      1500000   150       150.0     150
 reads_1.fq.gz       FASTQ    DNA    2500       567516    226       227.0     229
 reads_2.fq.gz       FASTQ    DNA    2500       560002    223       224.0     225


 # To markdown

 $ seqkit stats *.f{a,q}.gz -T | csvtk csv2md -t
 |file             |format|type|num_seqs|sum_len|min_len|avg_len|max_len|
 |:----------------|:-----|:---|:-------|:------|:------|:------|:------|
 |hairpin.fa.gz    |FASTA |RNA |28645   |2949871|39     |103.0  |2354   |
 |mature.fa.gz     |FASTA |RNA |35828   |781222 |15     |21.8   |34     |
 |Illimina1.8.fq.gz|FASTQ |DNA |10000   |1500000|150    |150.0  |150    |
 |reads_1.fq.gz    |FASTQ |DNA |2500    |567516 |226    |227.0  |229    |
 |reads_2.fq.gz    |FASTQ |DNA |2500    |560002 |223    |224.0  |225    |
    file format type num_seqs sum_len min_len avg_len max_len
    hairpin.fa.gz FASTA RNA 28645 2949871 39 103.0 2354
    mature.fa.gz FASTA RNA 35828 781222 15 21.8 34
    Illimina1.8.fq.gz FASTQ DNA 10000 1500000 150 150.0 150
    reads_1.fq.gz FASTQ DNA 2500 567516 226 227.0 229
    reads_2.fq.gz FASTQ DNA 2500 560002 223 224.0 225

  3. Extra information

     $ seqkit stats *.f{a,q}.gz -a
 file               format  type  num_seqs    sum_len  min_len  avg_len  max_len   Q1   Q2   Q3  sum_gap  N50  Q20(%)  Q30(%)  GC(%)
 hairpin.fa.gz      FASTA   RNA     28,645  2,949,871       39      103    2,354   76   91  111        0  101       0       0  45.77
 mature.fa.gz       FASTA   RNA     35,828    781,222       15     21.8       34   21   22   22        0   22       0       0   47.6
 Illimina1.8.fq.gz  FASTQ   DNA     10,000  1,500,000      150      150      150   75  150   75        0  150   96.16   89.71  49.91
 nanopore.fq.gz     FASTQ   DNA      4,000  1,798,723      153    449.7    6,006  271  318  391        0  395   40.79   12.63  46.66
 reads_1.fq.gz      FASTQ   DNA      2,500    567,516      226      227      229  227  227  227        0  227   91.24   86.62  53.63
 reads_2.fq.gz      FASTQ   DNA      2,500    560,002      223      224      225  224  224  224        0  224   91.06   87.66  54.77

  4. Parallelize counting files, it’s much faster for lots of small files, especially for files on SSD

     seqkit stats -j 10 refseq/virual/*.fna.gz

  5. Skip error

     $ seqkit stats tests/*
 [ERRO] tests/hairpin.fa.fai: fastx: invalid FASTA/Q format

 $ seqkit stats tests/* -e
 [WARN] tests/hairpin.fa.fai: fastx: invalid FASTA/Q format
 [WARN] tests/hairpin.fa.seqkit.fai: fastx: invalid FASTA/Q format
 [WARN] tests/miRNA.diff.gz: fastx: invalid FASTA/Q format
 [WARN] tests/test.sh: fastx: invalid FASTA/Q format
 file                     format  type  num_seqs    sum_len  min_len  avg_len  max_len
 tests/contigs.fa         FASTA   DNA          9         54        2        6       10
 tests/hairpin.fa         FASTA   RNA     28,645  2,949,871       39      103    2,354
 tests/Illimina1.5.fq     FASTQ   DNA          1        100      100      100      100
 tests/Illimina1.8.fq.gz  FASTQ   DNA     10,000  1,500,000      150      150      150
 tests/hairpin.fa.gz      FASTA   RNA     28,645  2,949,871       39      103    2,354
 tests/reads_1.fq.gz      FASTQ   DNA      2,500    567,516      226      227      229
 tests/mature.fa.gz       FASTA   RNA     35,828    781,222       15     21.8       34
 tests/reads_2.fq.gz      FASTQ   DNA      2,500    560,002      223      224      225

  6. Output basename instead of full path (-b/--basename)

sum

Usage

compute message digest for all sequences in FASTA/Q files

Attentions:
  1. Sequence headers and qualities are skipped, only sequences matter.
  2. The order of sequences records does not matter.
  3. Circular complete genomes are supported with the flag -c/--circular.
     - The same double-stranded genomes with different start positions or
       in reverse complement strand will not affect the result.
     - For single-stranded genomes like ssRNA genomes, use -s/--single-strand.
     - The message digest would change with different values of k-mer size.
  4. Multiple files are processed in parallel (-j/--threads).

Method:
  1. Converting the sequences to low cases, optionally removing gaps (-g).
  2. Computing the hash (xxhash) for all sequences or k-mers of a circular
     complete genome (-c/--circular).
  3. Sorting all hash values, for ignoring the order of sequences.
  4. Computing MD5 digest from the hash values, sequences length, and
     the number of sequences.

Following the seqhash in Poly (https://github.com/TimothyStiles/poly/),
We add meta information to the message digest, with the format of:

    seqkit.<version>_<seq type><seq structure><strand>_<kmer size>_<seq digest>

    <version>:       digest version
    <seq type>:      'D' for DNA, 'R' for RNA, 'P' for protein, 'N' for others
    <seq structure>: 'L' for linear sequence, 'C' for circular genome
    <strand>:        'D' for double-stranded, 'S' for single-stranded
    <kmer size>:     0 for linear sequence, other values for circular genome

Examples:

    seqkit.v0.1_DLS_k0_176250c8d1cde6c385397df525aa1a94    DNA.fq.gz
    seqkit.v0.1_PLS_k0_c244954e4960dd2a1409cd8ee53d92b9    Protein.fasta
    seqkit.v0.1_RLS_k0_0f1fb263f0c05a259ae179a61a80578d    single-stranded RNA.fasta

    seqkit.v0.1_DCD_k31_e59dad6d561f1f1f28ebf185c6f4c183   double-stranded-circular DNA.fasta
    seqkit.v0.1_DCS_k31_dd050490cd62ea5f94d73d4d636b7d60   single-stranded-circular DNA.fasta

Usage:
  seqkit sum [flags]

Flags:
  -a, --all                  show all information, including the sequences length and the number of sequences
  -b, --basename             only output basename of files
  -c, --circular             the file contains a single cicular genome sequence
  -G, --gap-letters string   gap letters to delete with the flag -g/--remove-gaps (default "- \t.*")
  -h, --help                 help for sum
  -k, --kmer-size int        k-mer size for processing circular genomes (default 1000)
  -g, --remove-gaps          remove gap characters set in the option -G/gap-letters
      --rna2dna              convert RNA to DNA
  -s, --single-strand        only consider the positive strand of a circular genome, e.g., ssRNA virus
                             genomes

Examples:

A, B, C, D are the same vircular genomes with different starting positions or strands:

$ cat virus-{A,B,C,D}.fasta
>seq
TGGTAGGGAGTTGAGTAGCATGGGTATAGTATAGTGTCATGATGCCAGATTTTAAAAAAA
>seq.revcom
TTTTTTTAAAATCTGGCATCATGACACTATACTATACCCATGCTACTCAACTCCCTACCA
>seq.new-start
GGTAGGGAGTTGAGTAGCATGGGTATAGTATAGTGTCATGATGCCAGATTTTAAAAAAAT
>seq.revcom.new-start
TTTTTTAAAATCTGGCATCATGACACTATACTATACCCATGCTACTCAACTCCCTACCAT

# cat to one file
$ cat virus-{A,B,C,D}.fasta > virues.fasta

# shuffle and rename
$ cat virus-{A,B,C,D}.fasta \
    | seqkit shuffle \
    | seqkit replace -p '.*' -r '{nr}' \
    | tee virues.shuffled.fasta
>1
TTTTTTAAAATCTGGCATCATGACACTATACTATACCCATGCTACTCAACTCCCTACCAT
>2
TGGTAGGGAGTTGAGTAGCATGGGTATAGTATAGTGTCATGATGCCAGATTTTAAAAAAA
>3
GGTAGGGAGTTGAGTAGCATGGGTATAGTATAGTGTCATGATGCCAGATTTTAAAAAAAT
>4
TTTTTTTAAAATCTGGCATCATGACACTATACTATACCCATGCTACTCAACTCCCTACCA

Sum of all files (the sequences order does not matter):

$ seqkit sum viru*.fasta
seqkit.v0.1_DLS_k0_9bbe0abefc26013dffdde952a6725b17    virues.fasta
seqkit.v0.1_DLS_k0_9bbe0abefc26013dffdde952a6725b17    virues.shuffled.fasta
seqkit.v0.1_DLS_k0_176250c8d1cde6c385397df525aa1a94    virus-A.fasta
seqkit.v0.1_DLS_k0_7a813339f9ae686b376b1df55cd596ca    virus-B.fasta
seqkit.v0.1_DLS_k0_0fd51028bfbfa85ddbdd2b86ef7bd1c1    virus-C.fasta
seqkit.v0.1_DLS_k0_88b1d20dd0fe0dbf41c00b075fee4e4e    virus-D.fasta

Circular genomes (the same genomes with different start positions or in reverse complement strand will not affect the result):

$ seqkit sum -c -k 21  virus-*.fasta
seqkit.v0.1_DCD_k21_7efd18ce33380268d3aa335ffd2dd1cc   virus-A.fasta
seqkit.v0.1_DCD_k21_7efd18ce33380268d3aa335ffd2dd1cc   virus-B.fasta
seqkit.v0.1_DCD_k21_7efd18ce33380268d3aa335ffd2dd1cc   virus-C.fasta
seqkit.v0.1_DCD_k21_7efd18ce33380268d3aa335ffd2dd1cc   virus-D.fasta

$ seqkit sum -c -k 51  virus-*.fasta
seqkit.v0.1_DCD_k51_39e267864fddeafd7a5cacd77e0a6a11   virus-A.fasta
seqkit.v0.1_DCD_k51_39e267864fddeafd7a5cacd77e0a6a11   virus-B.fasta
seqkit.v0.1_DCD_k51_39e267864fddeafd7a5cacd77e0a6a11   virus-C.fasta
seqkit.v0.1_DCD_k51_39e267864fddeafd7a5cacd77e0a6a11   virus-D.fasta
    
# collect files with the same genomes
$ seqkit sum -c -k 51  virus-*.fasta | csvtk fold -Ht -f 1 -v 2 
seqkit.v0.1_DCD_k51_39e267864fddeafd7a5cacd77e0a6a11   virus-A.fasta; virus-B.fasta; virus-C.fasta; virus-D.fasta

faidx

Usage

create FASTA index file and extract subsequence

This command is similar with "samtools faidx" but has some extra features:

  1. output full header line with the flag -f
  2. support regular expression as sequence ID with the flag -r
  3. if you have large number of IDs, you can use:
        seqkit faidx seqs.fasta -l IDs.txt

Attentions:
  1. The flag -U/--update-faidx is recommended to ensure the .fai file matches the FASTA file.

The definition of region is 1-based and with some custom design.

Examples:

 1-based index    1 2 3 4 5 6 7 8 9 10
negative index    0-9-8-7-6-5-4-3-2-1
           seq    A C G T N a c g t n
           1:1    A
           2:4      C G T
         -4:-2                c g t
         -4:-1                c g t n
         -1:-1                      n
          2:-2      C G T N a c g t
          1:-1    A C G T N a c g t n
          1:12    A C G T N a c g t n
        -12:-1    A C G T N a c g t n

Usage:
  seqkit faidx [flags] <fasta-file> [regions...]

Flags:
  -f, --full-head            print full header line instead of just ID. New fasta index file ending with
                             .seqkit.fai will be created
  -h, --help                 help for faidx
  -i, --ignore-case          ignore case
  -I, --immediate-output     print output immediately, do not use write buffer
  -l, --region-file string   file containing a list of regions
  -U, --update-faidx         update the fasta index file if it exists. Use this if you are not sure
                             whether the fasta file changed
  -r, --use-regexp           IDs are regular expression. But subseq region is not supported here.

Example

  1. common usage like samtools faidx

     $ seqkit faidx tests/hairpin.fa hsa-let-7a-1 hsa-let-7a-2
 >hsa-let-7a-1
 UGGGAUGAGGUAGUAGGUUGUAUAGUUUUAGGGUCACACCCACCACUGGGAGAUAACUAU
 ACAAUCUACUGUCUUUCCUA
 >hsa-let-7a-2
 AGGUUGAGGUAGUAGGUUGUAUAGUUUAGAAUUACAUCAAGGGAGAUAACUGUACAGCCU
 CCUAGCUUUCCU

  2. output full header, not supported by samtools faidx

     $ seqkit faidx tests/hairpin.fa hsa-let-7a-1 hsa-let-7a-2 -f
 >hsa-let-7a-1 MI0000060 Homo sapiens let-7a-1 stem-loop
 UGGGAUGAGGUAGUAGGUUGUAUAGUUUUAGGGUCACACCCACCACUGGGAGAUAACUAU
 ACAAUCUACUGUCUUUCCUA
 >hsa-let-7a-2 MI0000061 Homo sapiens let-7a-2 stem-loop
 AGGUUGAGGUAGUAGGUUGUAUAGUUUAGAAUUACAUCAAGGGAGAUAACUGUACAGCCU
 CCUAGCUUUCCU

  3. extract subsequence of specific region

     $ seqkit faidx tests/hairpin.fa hsa-let-7a-1:1-10
 >hsa-let-7a-1:1-10
 UGGGAUGAGG

 $ seqkit faidx tests/hairpin.fa hsa-let-7a-1:-10--1
 >hsa-let-7a-1:-10--1
 GUCUUUCCUA

 $ seqkit faidx tests/hairpin.fa hsa-let-7a-1:1
 >hsa-let-7a-1:1-1
 U

  4. supporting begin > start, i.e., returning reverse complement sequence, not supported by samtools faidx

     $ seqkit faidx tests/hairpin.fa hsa-let-7a-1:10-1
 >hsa-let-7a-1:10-1
 CCUCAUCCCA

  5. use regular expression

     $ seqkit faidx tests/hairpin.fa hsa -r | seqkit stats
 file  format  type  num_seqs  sum_len  min_len  avg_len  max_len
 -     FASTA   RNA      1,881  154,002       41     81.9      180

watch

Usage

monitoring and online histograms of sequence features

Usage:
  seqkit watch [flags]

Flags:
  -B, --bins int                  number of histogram bins (default -1)
  -W, --delay int                 sleep this many seconds after online plotting (default 1)
  -y, --dump                      print histogram data to stderr instead of plotting
  -f, --fields string             target fields, available values: ReadLen, MeanQual, GC, GCSkew
                                  (default "ReadLen")
  -h, --help                      help for watch
  -O, --img string                save histogram to this PDF/image file
  -H, --list-fields               print out a list of available fields
  -L, --log                       log10(x+1) transform numeric values
  -x, --pass                      pass through mode (write input to stdout)
  -p, --print-freq int            print/report after this many records (-1 for print after EOF) (default -1)
  -b, --qual-ascii-base int       ASCII BASE, 33 for Phred+33 (default 33)
  -Q, --quiet-mode                supress all plotting to stderr
  -R, --reset                     reset histogram after every report
  -v, --validate-seq              validate bases according to the alphabet
  -V, --validate-seq-length int   length of sequence to validate (0 for whole seq) (default 10000)

Examples

  1. Histogram of sequence length

     seqkit watch --fields ReadLen nanopore.fq.gz -O len.png

  2. Dynamic histogram of log sequence length

     seqkit watch --log --fields ReadLen nanopore.fq.gz

  3. Histogram of mean base qualities every 500 record, also saved as png

     seqkit watch -p 500 -O qhist.png -f MeanQual nanopore.fq.gz

sana

Usage

sanitize broken single line FASTQ files

Usage:
  seqkit sana [flags]

Flags:
  -A, --allow-gaps            allow gap character (-) in sequences
  -i, --format string         input and output format: fastq or fasta (default "fastq")
  -h, --help                  help for sana
  -I, --in-format string      input format: fastq or fasta
  -O, --out-format string     output format: fastq or fasta
  -b, --qual-ascii-base int   ASCII BASE, 33 for Phred+33 (default 33)

Examples

  1. Rescue usable reads from fastq file with malformed records.

     seqkit sana broken.fq.gz -o rescued.fq.gz

scat

Usage

real time recursive concatenation and streaming of fastx files

Usage:
  seqkit scat [flags]

Flags:
  -A, --allow-gaps            allow gap character (-) in sequences
  -d, --delta int             minimum size increase in kilobytes to trigger parsing (default 5)
  -D, --drop-time string      Notification drop interval (default "500ms")
  -f, --find-only             concatenate exisiting files and quit
  -i, --format string         input and output format: fastq or fasta (fastq) (default "fastq")
  -g, --gz-only               only look for gzipped files (.gz suffix)
  -h, --help                  help for scat
  -I, --in-format string      input format: fastq or fasta (fastq)
  -O, --out-format string     output format: fastq or fasta
  -b, --qual-ascii-base int   ASCII BASE, 33 for Phred+33 (default 33)
  -r, --regexp string         regexp for watched files, by default guessed from the input format
  -T, --time-limit string     quit after inactive for this time period
  -p, --wait-pid int          after process with this PID exited (default -1)

Examples

  1. Concatenate all fastq files recursively under a directory

     seqkit scat -j 4 -f fastq_dir > all_records.fq

  2. Watch a directory and stream fastq records in real time until interrupt is recieved and plot read lengths using seqkit watch:

     seqkit scat -j 4 fastq_dir | seqkit watch -f ReadLen -

  3. Watch a directory and stream fastq records in real time until there is no write activity under the directory for 5 seconds:

     seqkit scat -j 4 -T "5s" fastq_dir > all_records.fq

  4. Watch a directory and stream fastq records in real time until a process with a specified PID is alive:

     seqkit scat -j 4 -p $PID fastq_dir > all_records.fq

Notes: You might need to increase the ulimit allowance on open files if you intend to stream fastx records from a large number of files.

fq2fa

Usage

convert FASTQ to FASTA

Usage:
  seqkit fq2fa [flags]

Examples

seqkit fq2fa reads_1.fq.gz -o reads_1.fa.gz

fa2fq

Usage

retrieve corresponding FASTQ records by a FASTA file

Attention:
  1. We assume the FASTA file comes from the FASTQ file,
     so they share sequence IDs, and sequences in FASTA
     should be subseq of sequences in FASTQ file.

Usage:
  seqkit fa2fq [flags]

Flags:
  -f, --fasta-file string      FASTA file)
  -h, --help                   help for fa2fq
  -P, --only-positive-strand   only search on positive strand

fx2tab & tab2fx

Usage (fx2tab)

convert FASTA/Q to tabular format, and provide various information,
like sequence length, GC content/GC skew.

Attention:
  1. Fixed three columns (ID, sequence, quality) are outputted for either FASTA
     or FASTQ, except when flag -n/--name is on. This is for format compatibility.

Usage:
  seqkit fx2tab [flags]

Flags:
  -a, --alphabet               print alphabet letters
  -q, --avg-qual               print average quality of a read
  -B, --base-content strings   print base content. (case ignored, multiple values supported) e.g. -B AT -B N
  -C, --base-count strings     print base count. (case ignored, multiple values supported) e.g. -C AT -C N
  -I, --case-sensitive         calculate case sensitive base content/sequence hash
  -g, --gc                     print GC content
  -G, --gc-skew                print GC-Skew
  -H, --header-line            print header line
  -h, --help                   help for fx2tab
  -l, --length                 print sequence length
  -n, --name                   only print names (no sequences and qualities)
  -Q, --no-qual                only output two column even for FASTQ file
  -i, --only-id                print ID instead of full head
  -b, --qual-ascii-base int    ASCII BASE, 33 for Phred+33 (default 33)
  -s, --seq-hash               print hash (MD5) of sequence

Usage (tab2fx)

convert tabular format (first two/three columns) to FASTA/Q format

Usage:
  seqkit tab2fx [flags] 

Flags:
  -b, --buffer-size string            size of buffer, supported unit: K, M, G. You need increase the
                                      value when "bufio.Scanner: token too long" error reported (default
                                      "1G")
  -p, --comment-line-prefix strings   comment line prefix (default [#,//])
  -h, --help                          help for tab2fx

Examples

  1. Default output

     $ seqkit fx2tab hairpin.fa.gz | head -n 2
 cel-let-7 MI0000001 Caenorhabditis elegans let-7 stem-loop      UACACUGUGGAUCCGGUGAGGUAGUAGGUUGUAUAGUUUGGAAUAUUACCACCGGUGAACUAUGCAAUUUUCUACCUUACCGGAGACAGAACUCUUCGA
 cel-lin-4 MI0000002 Caenorhabditis elegans lin-4 stem-loop      AUGCUUCCGGCCUGUUCCCUGAGACCUCAAGUGUGAGUGUACUAUUGAUGCUUCACACCUGGGCUCUCCGGGUACCAGGACGGUUUGAGCAGAU

  2. Print sequence length, GC content, and only print names (no sequences), we could also print title line by flag -H.

     $ seqkit fx2tab hairpin.fa.gz -l -g -n -i -H | head -n 4 | csvtk -t -C '&' pretty
 #name       seq   qual   length   GC
 cel-let-7                99       43.43
 cel-lin-4                94       54.26
 cel-mir-1                96       40.62

  3. Use fx2tab and tab2fx in pipe

     $ zcat hairpin.fa.gz | seqkit fx2tab | seqkit tab2fx

 $ zcat reads_1.fq.gz | seqkit fx2tab | seqkit tab2fx

  4. Sort sequences by length (use seqkit sort -l)

     $ zcat hairpin.fa.gz \
     | seqkit fx2tab -l \
     | sort -t"`echo -e '\t'`" -n -k4,4 \
     | seqkit tab2fx
 >cin-mir-4129 MI0015684 Ciona intestinalis miR-4129 stem-loop
 UUCGUUAUUGGAAGACCUUAGUCCGUUAAUAAAGGCAUC
 >mmu-mir-7228 MI0023723 Mus musculus miR-7228 stem-loop
 UGGCGACCUGAACAGAUGUCGCAGUGUUCGGUCUCCAGU
 >cin-mir-4103 MI0015657 Ciona intestinalis miR-4103 stem-loop
 ACCACGGGUCUGUGACGUAGCAGCGCUGCGGGUCCGCUGU

 $ seqkit sort -l hairpin.fa.gz

    Sorting or filtering by GC (or other base by -flag -B) content could also achieved in similar way.

  5. Get first 1000 sequences (use seqkit head -n 1000)

     $ seqkit fx2tab hairpin.fa.gz | head -n 1000 | seqkit tab2fx

 $ seqkit fx2tab reads_1.fq.gz | head -n 1000 | seqkit tab2fx

Extension

After converting FASTA to tabular format with seqkit fx2tab, it could be handled with CSV/TSV tools, e.g. csvtk, a cross-platform, efficient and practical CSV/TSV toolkit

convert

Usage

convert FASTQ quality encoding between Sanger, Solexa and Illumina

Usage:
  seqkit convert [flags]

Flags:
  -d, --dry-run                         dry run
  -f, --force                           for Illumina-1.8+ -> Sanger, truncate scores > 40 to 40
      --from string                     source quality encoding. if not given, we'll guess it
  -h, --help                            help for convert
  -n, --nrecords int                    number of records for guessing quality encoding (default 1000)
  -N, --thresh-B-in-n-most-common int   threshold of 'B' in top N most common quality for guessing
                                        Illumina 1.5. (default 2)
  -F, --thresh-illumina1.5-frac float   threshold of faction of Illumina 1.5 in the leading N records
                                        (default 0.1)
      --to string                       target quality encoding (default "Sanger")

Examples:

Note that seqkit convert always output sequences.

The test dataset contains score 41 (J):

$ seqkit head -n 1 tests/Illimina1.8.fq.gz
@ST-E00493:56:H33MFALXX:4:1101:23439:1379 1:N:0:NACAACCA
NCGTGGAAAGACGCTAAGATTGTGATGTGCTTCCCTGACGATTACAACTGGCGTAAGGACGTTTTGCCTACCTATAAGGCTAACCGTAAGGGTTCTCGCAAGCCTGTAGGTTACAAGAGGTTCGTAGCCGAAGTGATGGCTGACTCACGG
+
#AAAFAAJFFFJJJ<JJJJJFFFJFJJJJJFJJAJJJFJJFJFJJJJFAFJ<JA<FFJ7FJJFJJAAJJJJ<JJJJJJJFJJJAJJJJJFJJ77<JJJJ-F7A-FJFFJJJJJJ<FFJ-<7FJJJFJJ)A7)7AA<7--)<-7F-A7FA<

By default, nothing changes when converting Illumina 1.8 to Sanger. A warning message show that source and target quality encoding match.

$ seqkit convert tests/Illimina1.8.fq.gz  | seqkit head -n 1
[INFO] possible quality encodings: [Illumina-1.8+]
[INFO] guessed quality encoding: Illumina-1.8+
[INFO] converting Illumina-1.8+ -> Sanger
[WARN] source and target quality encoding match.
@ST-E00493:56:H33MFALXX:4:1101:23439:1379 1:N:0:NACAACCA
NCGTGGAAAGACGCTAAGATTGTGATGTGCTTCCCTGACGATTACAACTGGCGTAAGGACGTTTTGCCTACCTATAAGGCTAACCGTAAGGGTTCTCGCAAGCCTGTAGGTTACAAGAGGTTCGTAGCCGAAGTGATGGCTGACTCACGG
+
#AAAFAAJFFFJJJ<JJJJJFFFJFJJJJJFJJAJJJFJJFJFJJJJFAFJ<JA<FFJ7FJJFJJAAJJJJ<JJJJJJJFJJJAJJJJJFJJ77<JJJJ-F7A-FJFFJJJJJJ<FFJ-<7FJJJFJJ)A7)7AA<7--)<-7F-A7FA<

When switching flag --force on, J (41) was converted to I (40).

$ seqkit convert tests/Illimina1.8.fq.gz -f | seqkit head -n 1
[INFO] possible quality encodings: [Illumina-1.8+]
[INFO] guessed quality encoding: Illumina-1.8+
[INFO] converting Illumina-1.8+ -> Sanger
@ST-E00493:56:H33MFALXX:4:1101:23439:1379 1:N:0:NACAACCA
NCGTGGAAAGACGCTAAGATTGTGATGTGCTTCCCTGACGATTACAACTGGCGTAAGGACGTTTTGCCTACCTATAAGGCTAACCGTAAGGGTTCTCGCAAGCCTGTAGGTTACAAGAGGTTCGTAGCCGAAGTGATGGCTGACTCACGG
+
#AAAFAAIFFFIII<IIIIIFFFIFIIIIIFIIAIIIFIIFIFIIIIFAFI<IA<FFI7FIIFIIAAIIII<IIIIIIIFIIIAIIIIIFII77<IIII-F7A-FIFFIIIIII<FFI-<7FIIIFII)A7)7AA<7--)<-7F-A7FA<

Other cases:

To Illumina-1.5.

$ seqkit convert tests/Illimina1.8.fq.gz --to Illumina-1.5+ | seqkit head -n 1
[INFO] possible quality encodings: [Illumina-1.8+]
[INFO] guessed quality encoding: Illumina-1.8+
[INFO] converting Illumina-1.8+ -> Illumina-1.5+
@ST-E00493:56:H33MFALXX:4:1101:23439:1379 1:N:0:NACAACCA
NCGTGGAAAGACGCTAAGATTGTGATGTGCTTCCCTGACGATTACAACTGGCGTAAGGACGTTTTGCCTACCTATAAGGCTAACCGTAAGGGTTCTCGCAAGCCTGTAGGTTACAAGAGGTTCGTAGCCGAAGTGATGGCTGACTCACGG
+
B```e``ieeeiii[iiiiieeeieiiiiieii`iiieiieieiiiie`ei[i`[eeiVeiieii``iiii[iiiiiiieiii`iiiiieiiVV[iiiiLeV`Leieeiiiiii[eeiL[VeiiieiiH`VHV``[VLLH[LVeL`Ve`[

To Illumina-1.5 and back to Sanger.

$ seqkit convert tests/Illimina1.8.fq.gz --to Illumina-1.5+ | seqkit convert | seqkit head -n 1
[INFO] possible quality encodings: [Illumina-1.8+]
[INFO] guessed quality encoding: Illumina-1.8+
[INFO] converting Illumina-1.8+ -> Illumina-1.5+
[INFO] possible quality encodings: [Illumina-1.5+]
[INFO] guessed quality encoding: Illumina-1.5+
[INFO] converting Illumina-1.5+ -> Sanger
@ST-E00493:56:H33MFALXX:4:1101:23439:1379 1:N:0:NACAACCA
NCGTGGAAAGACGCTAAGATTGTGATGTGCTTCCCTGACGATTACAACTGGCGTAAGGACGTTTTGCCTACCTATAAGGCTAACCGTAAGGGTTCTCGCAAGCCTGTAGGTTACAAGAGGTTCGTAGCCGAAGTGATGGCTGACTCACGG
+
!AAAFAAJFFFJJJ<JJJJJFFFJFJJJJJFJJAJJJFJJFJFJJJJFAFJ<JA<FFJ7FJJFJJAAJJJJ<JJJJJJJFJJJAJJJJJFJJ77<JJJJ-F7A-FJFFJJJJJJ<FFJ-<7FJJJFJJ)A7)7AA<7--)<-7F-A7FA<

Checking encoding

$ seqkit convert tests/Illimina1.8.fq.gz --from Solexa
[INFO] converting Solexa -> Sanger
[ERRO] seq: invalid Solexa quality

Real Illumina 1.5+ data

$ seqkit seq tests/Illimina1.5.fq
@HWI-EAS209_0006_FC706VJ:5:58:5894:21141#ATCACG/1
TTAATTGGTAAATAAATCTCCTAATAGCTTAGATNTTACCTTNNNNNNNNNNTAGTTTCTTGAGATTTGTTGGGGGAGACATTTTTGTGATTGCCTTGAT
+
efcfffffcfeefffcffffffddf`feed]`]_Ba_^__[YBBBBBBBBBBRTT\]][]dddd`ddd^dddadd^BBBBBBBBBBBBBBBBBBBBBBBB

$ seqkit convert tests/Illimina1.5.fq | seqkit head -n 1
[INFO] possible quality encodings: [Illumina-1.5+]
[INFO] guessed quality encoding: Illumina-1.5+
[INFO] converting Illumina-1.5+ -> Sanger
@HWI-EAS209_0006_FC706VJ:5:58:5894:21141#ATCACG/1
TTAATTGGTAAATAAATCTCCTAATAGCTTAGATNTTACCTTNNNNNNNNNNTAGTTTCTTGAGATTTGTTGGGGGAGACATTTTTGTGATTGCCTTGAT
+
FGDGGGGGDGFFGGGDGGGGGGEEGAGFFE>A>@!B@?@@<:!!!!!!!!!!355=>><>EEEEAEEE?EEEBEE?!!!!!!!!!!!!!!!!!!!!!!!!

translate

Usage

translate DNA/RNA to protein sequence (supporting ambiguous bases)

Note:

  1. This command supports codons containing any ambiguous base.
     Please switch on flag -L INT for details. e.g., for standard table:

        ACN -> T
        CCN -> P
        CGN -> R
        CTN -> L
        GCN -> A
        GGN -> G
        GTN -> V
        TCN -> S
        
        MGR -> R
        YTR -> L

Translate Tables/Genetic Codes:

    # https://www.ncbi.nlm.nih.gov/Taxonomy/taxonomyhome.html/index.cgi?chapter=tgencodes

     1: The Standard Code
     2: The Vertebrate Mitochondrial Code
     3: The Yeast Mitochondrial Code
     4: The Mold, Protozoan, and Coelenterate Mitochondrial Code and the Mycoplasma/Spiroplasma Code
     5: The Invertebrate Mitochondrial Code
     6: The Ciliate, Dasycladacean and Hexamita Nuclear Code
     9: The Echinoderm and Flatworm Mitochondrial Code
    10: The Euplotid Nuclear Code
    11: The Bacterial, Archaeal and Plant Plastid Code
    12: The Alternative Yeast Nuclear Code
    13: The Ascidian Mitochondrial Code
    14: The Alternative Flatworm Mitochondrial Code
    16: Chlorophycean Mitochondrial Code
    21: Trematode Mitochondrial Code
    22: Scenedesmus obliquus Mitochondrial Code
    23: Thraustochytrium Mitochondrial Code
    24: Pterobranchia Mitochondrial Code
    25: Candidate Division SR1 and Gracilibacteria Code
    26: Pachysolen tannophilus Nuclear Code
    27: Karyorelict Nuclear
    28: Condylostoma Nuclear
    29: Mesodinium Nuclear
    30: Peritrich Nuclear
    31: Blastocrithidia Nuclear

Usage:
  seqkit translate [flags]

Flags:
  -x, --allow-unknown-codon                     translate unknown code to 'X'. And you may not use flag
                                                --trim which removes 'X'
  -F, --append-frame                            append frame information to sequence ID
      --clean                                   change all STOP codon positions from the '*' character
                                                to 'X' (an unknown residue)
  -f, --frame strings                           frame(s) to translate, available value: 1, 2, 3, -1, -2,
                                                -3, and 6 for all six frames (default [1])
  -h, --help                                    help for translate
  -M, --init-codon-as-M                         translate initial codon at beginning to 'M'
  -l, --list-transl-table int                   show details of translate table N, 0 for all (default -1)
  -L, --list-transl-table-with-amb-codons int   show details of translate table N (including ambigugous
                                                codons), 0 for all.  (default -1)
  -T, --transl-table int                        translate table/genetic code, type 'seqkit translate
                                                --help' for more details (default 1)
      --trim                                    remove all 'X' and '*' characters from the right end of
                                                the translation

Examples

  1. common usage

     $ seqkit translate tests/mouse-p53-cds.fna
 >lcl|AB021961.1_cds_BAA82344.1_1 [gene=p53] [protein=P53] [protein_id=BAA82344.1] [location=101..1273] [gbkey=CDS]
 MTAMEESQSDISLELPLSQETFSGLWKLLPPEDILPSPHCMDDLLLPQDVEEFFEGPSEA
 LRVSGAPAAQDPVTETPGPVAPAPATPWPLSSFVPSQKTYQGNYGFHLGFLQSGTAKSVM
 CTYSPPLNKLFCQLAKTCPVQLWVSATPPAGSRVRAMAIYKKSQHMTEVVRRCPHHERCS
 DGDGLAPPQHRIRVEGNLYPEYLEDRQTFRHSVVVPYEPPEAGSEYTTIHYKYMCNSSCM
 GGMNRRPILTIITLEDSSGNLLGRDSFEVRVCACPGRDRRTEEENFRKKEVLCPELPPGS
 AKRALPTCTSASPPQKKKPLDGEYFTLKIRGRKRFEMFRELNEALELKDAHATEESGDSR
 AHSSYLKTKKGQSTSRHKKTMVKKVGPDSD*

  2. trim the *

     $ seqkit translate tests/mouse-p53-cds.fna --trim
 >lcl|AB021961.1_cds_BAA82344.1_1 [gene=p53] [protein=P53] [protein_id=BAA82344.1] [location=101..1273] [gbkey=CDS]
 MTAMEESQSDISLELPLSQETFSGLWKLLPPEDILPSPHCMDDLLLPQDVEEFFEGPSEA
 LRVSGAPAAQDPVTETPGPVAPAPATPWPLSSFVPSQKTYQGNYGFHLGFLQSGTAKSVM
 CTYSPPLNKLFCQLAKTCPVQLWVSATPPAGSRVRAMAIYKKSQHMTEVVRRCPHHERCS
 DGDGLAPPQHRIRVEGNLYPEYLEDRQTFRHSVVVPYEPPEAGSEYTTIHYKYMCNSSCM
 GGMNRRPILTIITLEDSSGNLLGRDSFEVRVCACPGRDRRTEEENFRKKEVLCPELPPGS
 AKRALPTCTSASPPQKKKPLDGEYFTLKIRGRKRFEMFRELNEALELKDAHATEESGDSR
 AHSSYLKTKKGQSTSRHKKTMVKKVGPDSD

  3. different translate table

     $ cat tests/Lactococcus-lactis-phage-BK5-T-ORF25.fasta \
     | seqkit translate -T 11 --trim
 >CAC80166.1 hypothetical protein [Lactococcus phage BK5-T]
 MEEQAWREVLERLARIETKLDNYETVRDKAERALLIAQSNAKLIEKMEANNKWAWGFMLT
 LAVTVIGYLFTKIRF

  4. different frame

     $ cat tests/Lactococcus-lactis-phage-BK5-T-ORF25.fasta \
     | seqkit translate -T 11 --frame -1
 >CAC80166.1 hypothetical protein [Lactococcus phage BK5-T]
 SESNFSE*ITNNSYGKSKHKAPSPLIISFHFFYKFRI*LSY*ERSFCFISNCFIVI*LCF
 NSS*TFEDFSPCLFLH

 $ cat tests/Lactococcus-lactis-phage-BK5-T-ORF25.fasta \
     | seqkit seq -r -p \
     | seqkit translate -T 11 --frame -1
 >CAC80166.1 hypothetical protein [Lactococcus phage BK5-T]
 MEEQAWREVLERLARIETKLDNYETVRDKAERALLIAQSNAKLIEKMEANNKWAWGFMLT
 LAVTVIGYLFTKIRF*

  5. show details of translate table 1

     $ seqkit translate -l 1
 The Standard Code (transl_table=1)
 Source: https://www.ncbi.nlm.nih.gov/Taxonomy/taxonomyhome.html/index.cgi?chapter=tgencodes#SG1

 Initiation Codons:
     ATG, CTG, TTG

 Stop Codons:
     TAA, TAG, TGA

 Stranslate Table:
     AAA: K, AAC: N, AAG: K, AAT: N
     ACA: T, ACC: T, ACG: T, ACT: T
     AGA: R, AGC: S, AGG: R, AGT: S
     ATA: I, ATC: I, ATG: M, ATT: I

     CAA: Q, CAC: H, CAG: Q, CAT: H
     CCA: P, CCC: P, CCG: P, CCT: P
     CGA: R, CGC: R, CGG: R, CGT: R
     CTA: L, CTC: L, CTG: L, CTT: L

     GAA: E, GAC: D, GAG: E, GAT: D
     GCA: A, GCC: A, GCG: A, GCT: A
     GGA: G, GGC: G, GGG: G, GGT: G
     GTA: V, GTC: V, GTG: V, GTT: V

     TAA: *, TAC: Y, TAG: *, TAT: Y
     TCA: S, TCC: S, TCG: S, TCT: S
     TGA: *, TGC: C, TGG: W, TGT: C
     TTA: L, TTC: F, TTG: L, TTT: F

  6. show details of translate table 1, including ambigugous codons

     $ seqkit translate -L 1
 The Standard Code (transl_table=1)
 Source: https://www.ncbi.nlm.nih.gov/Taxonomy/taxonomyhome.html/index.cgi?chapter=tgencodes#SG1

 Initiation Codons:
     ATG, CTG, TTG

 Stop Codons:
     TAA, TAG, TGA

 Stranslate Table:
     AAA: K, AAC: N, AAG: K, AAR: K, AAT: N, AAY: N
     ACA: T, ACC: T, ACM: T, ACG: T, ACR: T, ACS: T, ACV: T, ACT: T, ACW: T, ACY: T, ACH: T, ACK: T, ACD: T, ACB: T, ACN: T
     AGA: R, AGC: S, AGG: R, AGR: R, AGT: S, AGY: S
     ATA: I, ATC: I, ATM: I, ATG: M, ATT: I, ATW: I, ATY: I, ATH: I

     CAA: Q, CAC: H, CAG: Q, CAR: Q, CAT: H, CAY: H
     CCA: P, CCC: P, CCM: P, CCG: P, CCR: P, CCS: P, CCV: P, CCT: P, CCW: P, CCY: P, CCH: P, CCK: P, CCD: P, CCB: P, CCN: P
     CGA: R, CGC: R, CGM: R, CGG: R, CGR: R, CGS: R, CGV: R, CGT: R, CGW: R, CGY: R, CGH: R, CGK: R, CGD: R, CGB: R, CGN: R
     CTA: L, CTC: L, CTM: L, CTG: L, CTR: L, CTS: L, CTV: L, CTT: L, CTW: L, CTY: L, CTH: L, CTK: L, CTD: L, CTB: L, CTN: L

     MGA: R, MGG: R, MGR: R

     GAA: E, GAC: D, GAG: E, GAR: E, GAT: D, GAY: D
     GCA: A, GCC: A, GCM: A, GCG: A, GCR: A, GCS: A, GCV: A, GCT: A, GCW: A, GCY: A, GCH: A, GCK: A, GCD: A, GCB: A, GCN: A
     GGA: G, GGC: G, GGM: G, GGG: G, GGR: G, GGS: G, GGV: G, GGT: G, GGW: G, GGY: G, GGH: G, GGK: G, GGD: G, GGB: G, GGN: G
     GTA: V, GTC: V, GTM: V, GTG: V, GTR: V, GTS: V, GTV: V, GTT: V, GTW: V, GTY: V, GTH: V, GTK: V, GTD: V, GTB: V, GTN: V

     TAA: *, TAC: Y, TAG: *, TAR: *, TAT: Y, TAY: Y
     TCA: S, TCC: S, TCM: S, TCG: S, TCR: S, TCS: S, TCV: S, TCT: S, TCW: S, TCY: S, TCH: S, TCK: S, TCD: S, TCB: S, TCN: S
     TGA: *, TGC: C, TGG: W, TGT: C, TGY: C
     TRA: *
     TTA: L, TTC: F, TTG: L, TTR: L, TTT: F, TTY: F

     YTA: L, YTG: L, YTR: L

grep

Usage

search sequences by ID/name/sequence/sequence motifs, mismatch allowed

Attentions:

  0. By default, we match sequence ID with patterns, use "-n/--by-name"
     for matching full name instead of just ID.
  1. Unlike POSIX/GNU grep, we compare the pattern to the whole target
     (ID/full header) by default. Please switch "-r/--use-regexp" on
     for partly matching.
  2. When searching by sequences, it's partly matching, and both positive
     and negative strands are searched.
     Mismatch is allowed using flag "-m/--max-mismatch", you can increase
     the value of "-j/--threads" to accelerate processing.
  3. Degenerate bases/residues like "RYMM.." are also supported by flag -d.
     But do not use degenerate bases/residues in regular expression, you need
     convert them to regular expression, e.g., change "N" or "X"  to ".".
  4. When providing search patterns (motifs) via flag '-p',
     please use double quotation marks for patterns containing comma, 
     e.g., -p '"A{2,}"' or -p "\"A{2,}\"". Because the command line argument
     parser accepts comma-separated-values (CSV) for multiple values (motifs).
     Patterns in file do not follow this rule.
  5. The order of sequences in result is consistent with that in original
     file, not the order of the query patterns. 
     But for FASTA file, you can use:
        seqkit faidx seqs.fasta --infile-list IDs.txt
  6. For multiple patterns, you can either set "-p" multiple times, i.e.,
     -p pattern1 -p pattern2, or give a file of patterns via "-f/--pattern-file".

You can specify the sequence region for searching with the flag -R (--region).
The definition of region is 1-based and with some custom design.

Examples:

 1-based index    1 2 3 4 5 6 7 8 9 10
negative index    0-9-8-7-6-5-4-3-2-1
           seq    A C G T N a c g t n
           1:1    A
           2:4      C G T
         -4:-2                c g t
         -4:-1                c g t n
         -1:-1                      n
          2:-2      C G T N a c g t
          1:-1    A C G T N a c g t n
          1:12    A C G T N a c g t n
        -12:-1    A C G T N a c g t n

Usage:
  seqkit grep [flags]

Flags:
  -n, --by-name                match by full name instead of just ID
  -s, --by-seq                 search subseq on seq, both positive and negative strand are searched, and
                               mismatch allowed using flag -m/--max-mismatch
  -c, --circular               circular genome
  -C, --count                  just print a count of matching records. with the -v/--invert-match flag,
                               count non-matching records
  -d, --degenerate             pattern/motif contains degenerate base
      --delete-matched         delete a pattern right after being matched, this keeps the firstly
                               matched data and speedups when using regular expressions
  -h, --help                   help for grep
  -i, --ignore-case            ignore case
  -I, --immediate-output       print output immediately, do not use write buffer
  -v, --invert-match           invert the sense of matching, to select non-matching records
  -m, --max-mismatch int       max mismatch when matching by seq. For large genomes like human genome,
                               using mapping/alignment tools would be faster
  -P, --only-positive-strand   only search on positive strand
  -p, --pattern strings        search pattern (multiple values supported. Attention: use double
                               quotation marks for patterns containing comma, e.g., -p '"A{2,}"')
  -f, --pattern-file string    pattern file (one record per line)
  -R, --region string          specify sequence region for searching. e.g 1:12 for first 12 bases,
                               -12:-1 for last 12 bases
  -r, --use-regexp             patterns are regular expression

Examples

  1. Searching with list of sequence IDs (do not contain whitespace)

     $ seqkit grep -f id.txt seqs.fq.gz -o result.fq.gz
    
 # ignore case
 $ seqkit grep -i -f id.txt seqs.fq.gz -o result.fq.gz

  2. Just print the matched number, like GNU grep (grep -c)

     $ seqkit grep -f id.txt seqs.fq.gz -C

  3. Serching non-canonical sequence IDs, Using --id-regexp to capture IDs. Refer to section Sequence ID and seqkit seq for examples.

  4. Searching with list of sequence names (they may contain whitespace).

     $ seqkit grep -n -f name.txt seqs.fa.gz -o result.fa.gz

  5. Useq -r/--use-regexp for partly matching, but this may produce “false positive” matches. For example, seq_1 matches seq_10 with -nri.

  6. Extract human hairpins (i.e. sequences with name starting with hsa)

     $ zcat hairpin.fa.gz | seqkit grep -r -p ^hsa
 >hsa-let-7a-1 MI0000060 Homo sapiens let-7a-1 stem-loop
 UGGGAUGAGGUAGUAGGUUGUAUAGUUUUAGGGUCACACCCACCACUGGGAGAUAACUAU
 ACAAUCUACUGUCUUUCCUA
 >hsa-let-7a-2 MI0000061 Homo sapiens let-7a-2 stem-loop
 AGGUUGAGGUAGUAGGUUGUAUAGUUUAGAAUUACAUCAAGGGAGAUAACUGUACAGCCU
 CCUAGCUUUCCU

  7. Remove human and mice hairpins (invert match with -v)

     $ zcat hairpin.fa.gz | seqkit grep -r -p ^hsa -p ^mmu -v

  8. Extract new entries by information from miRNA.diff.gz

    1. Get IDs of new entries.

       $ zcat miRNA.diff.gz | grep ^# -v | grep NEW | cut -f 2 > list
 $ more list
 cfa-mir-486
 cfa-mir-339-1
 pmi-let-7

    2. Extract by ID list file

       $ zcat hairpin.fa.gz | seqkit grep -f list > new.fa

  9. Extract sequences containing AGGCG

     $ cat hairpin.fa.gz | seqkit grep -s -i -p aggcg

  10. Circular genome

     $ echo -e ">seq\nACGTTGCA" 
 >seq
 ACGTTGCA
    
 $ echo -e ">seq\nACGTTGCA"  | seqkit grep -s -i -p AA
    
 $ echo -e ">seq\nACGTTGCA"  | seqkit grep -s -i -p AA -c
 >seq
 ACGTTGCA

  11. Extract sequences containing AGGCG (allow mismatch)

     $ time cat hairpin.fa.gz | seqkit grep -s -i -p aggcg | seqkit stats
 file  format  type  num_seqs  sum_len  min_len  avg_len  max_len
 -     FASTA   RNA      1,181  145,789       49    123.4    2,354

 real    0m0.058s
 user    0m0.100s
 sys     0m0.017s

 $ time zcat hairpin.fa.gz | seqkit grep -s -i -p aggcg -m 1 | seqkit stats
 file  format  type  num_seqs    sum_len  min_len  avg_len  max_len
 -     FASTA   RNA     22,290  2,375,819       39    106.6    2,354

 real    0m1.081s
 user    0m1.305s
 sys     0m0.158s

  12. Extract sequences starting with AGGCG

     $ zcat hairpin.fa.gz | seqkit grep -s -r -i -p ^aggcg

  13. Extract sequences with TTSAA (AgsI digest site) in SEQUENCE. Base S stands for C or G.

     $ zcat hairpin.fa.gz | seqkit grep -s -d -i -p TTSAA

    It’s equal to but simpler than:

     $ zcat hairpin.fa.gz | seqkit grep -s -r -i -p TT[CG]AA

  14. Specify sequence regions for searching. e.g., leading 30 bases.

     $ seqkit grep -s -R 1:30 -i -r -p GCTGG

locate

Usage

locate subsequences/motifs, mismatch allowed

Attentions:

  1. Motifs could be EITHER plain sequence containing "ACTGN" OR regular
     expression like "A[TU]G(?:.{3})+?[TU](?:AG|AA|GA)" for ORFs.     
  2. Degenerate bases/residues like "RYMM.." are also supported by flag -d.
     But do not use degenerate bases/residues in regular expression, you need
     convert them to regular expression, e.g., change "N" or "X"  to ".".
  3. When providing search patterns (motifs) via flag '-p',
     please use double quotation marks for patterns containing comma, 
     e.g., -p '"A{2,}"' or -p "\"A{2,}\"". Because the command line argument
     parser accepts comma-separated-values (CSV) for multiple values (motifs).
     Patterns in file do not follow this rule.     
  4. Mismatch is allowed using flag "-m/--max-mismatch",
     you can increase the value of "-j/--threads" to accelerate processing.
  5. When using flag --circular, end position of matched subsequence that 
     crossing genome sequence end would be greater than sequence length.

Usage:
  seqkit locate [flags]

Flags:
      --bed                       output in BED6 format
  -c, --circular                  circular genome. type "seqkit locate -h" for details
  -d, --degenerate                pattern/motif contains degenerate base
      --gtf                       output in GTF format
  -h, --help                      help for locate
  -M, --hide-matched              do not show matched sequences
  -i, --ignore-case               ignore case
  -I, --immediate-output          print output immediately, do not use write buffer
  -m, --max-mismatch int          max mismatch when matching by seq. For large genomes like human
                                  genome, using mapping/alignment tools would be faster
  -G, --non-greedy                non-greedy mode, faster but may miss motifs overlapping with others
  -P, --only-positive-strand      only search on positive strand
  -p, --pattern strings           pattern/motif (multiple values supported. Attention: use double
                                  quotation marks for patterns containing comma, e.g., -p '"A{2,}"')
  -f, --pattern-file string       pattern/motif file (FASTA format)
  -F, --use-fmi                   use FM-index for much faster search of lots of sequence patterns
  -r, --use-regexp                patterns/motifs are regular expression
  -V, --validate-seq-length int   length of sequence to validate (0 for whole seq) (default 10000)

Examples

  1. Locating subsequences (mismatch allowed)

     $ cat t.fa
 >seq
 agctggagctacc

 $ cat t.fa \
   | seqkit locate -p agc \
   | csvtk pretty -t
 seqID   patternName   pattern   strand   start   end   matched
 seq     agc           agc       +        1       3     agc
 seq     agc           agc       +        7       9     agc
 seq     agc           agc       -        8       10    agc
 seq     agc           agc       -        2       4     agc

 # do not show matched sequences
 $ cat t.fa \
   | seqkit locate -p agc -M \
   | csvtk pretty -t
 seqID   patternName   pattern   strand   start   end
 seq     agc           agc       +        1       3
 seq     agc           agc       +        7       9
 seq     agc           agc       -        8       10
 seq     agc           agc       -        2       4

 # max mismatch: 1
 $ cat t.fa \
   | seqkit locate -p agc -m 1 \
   | csvtk pretty -t
 seqID   patternName     pattern strand  start   end    matched
 seq     agc           agc       +        1       3     agc
 seq     agc           agc       +        7       9     agc
 seq     agc           agc       +        11      13    acc
 seq     agc           agc       -        8       10    agc
 seq     agc           agc       -        2       4     agc

 # max mismatch: 2
 $ cat t.fa \
   | seqkit locate -p agc -m 2 \
   | csvtk pretty -t
 seqID   patternName   pattern   strand   start   end   matched
 seq     agc           agc       +        1       3     agc
 seq     agc           agc       +        4       6     tgg
 seq     agc           agc       +        5       7     gga
 seq     agc           agc       +        7       9     agc
 seq     agc           agc       +        10      12    tac
 seq     agc           agc       +        11      13    acc
 seq     agc           agc       -        11      13    ggt
 seq     agc           agc       -        8       10    agc
 seq     agc           agc       -        6       8     ctc
 seq     agc           agc       -        5       7     tcc
 seq     agc           agc       -        2       4     agc

  2. Locate ORFs.

     $ zcat hairpin.fa.gz \
     | seqkit locate -i -p "A[TU]G(?:.{3})+?[TU](?:AG|AA|GA)" -r \
     | head -n 4 \
     | csvtk pretty -t
 seqID       patternName                        pattern                            strand   start   end   matched
 cel-lin-4   A[TU]G(?:.{3})+?[TU](?:AG|AA|GA)   A[TU]G(?:.{3})+?[TU](?:AG|AA|GA)   +        1       36    AUGCUUCCGGCCUGUUCCCUGAGACCUCAAGUGUGA
 cel-mir-1   A[TU]G(?:.{3})+?[TU](?:AG|AA|GA)   A[TU]G(?:.{3})+?[TU](?:AG|AA|GA)   +        54      95    AUGGAUAUGGAAUGUAAAGAAGUAUGUAGAACGGGGUGGUAG
 cel-mir-1   A[TU]G(?:.{3})+?[TU](?:AG|AA|GA)   A[TU]G(?:.{3})+?[TU](?:AG|AA|GA)   -        43      51    AUGAUAUAG

  3. Locate Motif.

     $ zcat hairpin.fa.gz \
     | seqkit locate -i -d -p AUGGACUN \
     | head -n 4 \
     | csvtk pretty -t 
 seqID         patternName   pattern    strand   start   end   matched
 cel-mir-58a   AUGGACUN      AUGGACUN   +        81      88    AUGGACUG
 ath-MIR163    AUGGACUN      AUGGACUN   -        122     129   AUGGACUC
 cel-mir-270   AUGGACUN      AUGGACUN   +        84      91    AUGGACUG

  4. Output in GTF or BED6 format, which you can use in seqkit subseq

     $ zcat hairpin.fa.gz | seqkit locate -i -d -p AUGGACUN --bed
 cel-mir-58a     80      88      AUGGACUN        0       +
 ath-MIR163      121     129     AUGGACUN        0       -

 $ zcat hairpin.fa.gz | seqkit locate -i -d -p AUGGACUN --gtf
 cel-mir-58a     SeqKit  location        81      88      0       +       .       gene_id "AUGGACUN";
 ath-MIR163      SeqKit  location        122     129     0       -       .       gene_id "AUGGACUN";

  5. Greedy mode (default)

      $ echo -e '>seq\nACGACGACGA' | seqkit locate -p ACGA | csvtk -t pretty
  seqID   patternName   pattern   strand   start   end   matched
  seq     ACGA          ACGA      +        1       4     ACGA
  seq     ACGA          ACGA      +        4       7     ACGA
  seq     ACGA          ACGA      +        7       10    ACGA

  6. Non-greedy mode (-G)

     $ echo -e '>seq\nACGACGACGA' | seqkit locate -p ACGA -G | csvtk -t pretty
 seqID   patternName   pattern   strand   start   end   matched
 seq     ACGA          ACGA      +        1       4     ACGA
 seq     ACGA          ACGA      +        7       10    ACGA

  7. Circular genome. Note that end position of matched subsequence that crossing genome sequence end would be greater than sequence length.

     $ echo -e ">seq\nACGTTGCA"
 >seq
 ACGTTGCA

 $ echo -e ">seq\nACGTTGCA" \
     | seqkit locate -i -p aa
 seqID   patternName     pattern strand  start   end     matched
 seq     aa      aa      -       4       5       aa
    
 $ echo -e ">seq\nACGTTGCA" \
     | seqkit locate -i -p aa -c \
     | csvtk pretty -t
 seqID   patternName   pattern   strand   start   end   matched
 seq     aa            aa        +        8       9     aa
 seq     aa            aa        -        4       5     aa

fish

Usage

look for short sequences in larger sequences using local alignment

Attention:
  1. output coordinates are BED-like 0-based, left-close and right-open.
  2. alignment information are printed to STDERR.

Usage:
  seqkit fish [flags]

Flags:
  -a, --all                       search all
  -p, --aln-params string         alignment parameters in format
                                  "<match>,<mismatch>,<gap_open>,<gap_extend>" (default "4,-4,-2,-1")
  -h, --help                      help for fish
  -i, --invert                    print out references not matching with any query
  -q, --min-qual float            minimum mapping quality (default 5)
  -b, --out-bam string            save aligmnets to this BAM file (memory intensive)
  -x, --pass                      pass through mode (write input to stdout)
  -g, --print-aln                 print sequence alignments
  -D, --print-desc                print full sequence header
  -f, --query-fastx string        query fasta
  -F, --query-sequences string    query sequences
  -r, --ranges string             target ranges, for example: ":10,30:40,-20:"
  -s, --stranded                  search + strand only
  -v, --validate-seq              validate bases according to the alphabet
  -V, --validate-seq-length int   length of sequence to validate (0 for whole seq) (default 10000)

Examples

  1. Find best local alignment of a short sequence in reads in a fasta file, print results as tabular

     $ seqkit fish -q 4.7 -F GGCGGCTGTGACC -g mouse-p53-cds.fna

  2. Compare to seqkit locate:

     $ echo -e '>seq\nACGACGACGA' \
     | seqkit locate -p ACGA -G | csvtk -t pretty
 seqID   patternName   pattern   strand   start   end   matched
 seq     ACGA          ACGA      +        1       4     ACGA
 seq     ACGA          ACGA      +        7       10    ACGA
    
 $ echo -e '>seq\nACGACGACGA' \
     | seqkit fish -F ACGA -a 2>&1 | csvtk -t pretty 
 Ref   RefStart   RefEnd   Query   QueryStart   QueryEnd   Strand   MapQual   RawScore   Acc      ClipAcc   QueryCov
 seq   6          10       q0      0            4          +        60.00     16         100.00   100.00    100.00
 seq   0          4        q0      0            4          +        60.00     16         100.00   100.00    100.00

  3. Find all local alignment of a short sequences in reads in a fasta file, print results as tabular and save as BAM

     seqkit fish -a -q 4.67 -f query.fas -b alignments.bam -g mouse-p53-cds.fna

amplicon

Usage

extract amplicon (or specific region around it) via primer(s).

Attentions:
  1. Only one (the longest) matching location is returned for every primer pair.
  2. Mismatch is allowed, but the mismatch location (5' or 3') is not controled. 
     You can increase the value of "-j/--threads" to accelerate processing.
  3. Degenerate bases/residues like "RYMM.." are also supported.
     But do not use degenerate bases/residues in regular expression, you need
     convert them to regular expression, e.g., change "N" or "X"  to ".".

Examples:
  0. no region given.
  
                    F
        -----===============-----
             F             R
        -----=====-----=====-----
             
             ===============         amplicon

  1. inner region (-r x:y).

                    F
        -----===============-----
             1 3 5                    x/y
                      -5-3-1          x/y
             F             R
        -----=====-----=====-----     x:y
        
             ===============          1:-1
             =======                  1:7
               =====                  3:7
                  =====               6:10
                  =====             -10:-6
                     =====           -7:-3
                                     -x:y (invalid)
                    
  2. flanking region (-r x:y -f)
        
                    F
        -----===============-----
         -3-1                        x/y
                            1 3 5    x/y
             F             R
        -----=====-----=====-----
        
        =====                        -5:-1
        ===                          -5:-3
                            =====     1:5
                              ===     3:5
            =================        -1:1
        =========================    -5:5
                                      x:-y (invalid)

Usage:
  seqkit amplicon [flags]

Flags:
      --bed                    output in BED6+1 format with amplicon as the 7th column
  -f, --flanking-region        region is flanking region
  -F, --forward string         forward primer (5'-primer-3'), degenerate bases allowed
  -h, --help                   help for amplicon
  -I, --immediate-output       print output immediately, do not use write buffer
  -m, --max-mismatch int       max mismatch when matching primers, no degenerate bases allowed
  -P, --only-positive-strand   only search on positive strand
  -M, --output-mismatches      append the total mismatches and mismatches of 5' end and 3' end
  -p, --primer-file string     3- or 2-column tabular primer file, with first column as primer name
  -r, --region string          specify region to return. type "seqkit amplicon -h" for detail
  -R, --reverse string         reverse primer (5'-primer-3'), degenerate bases allowed
  -u, --save-unmatched         also save records that do not match any primer
  -s, --strict-mode            strict mode, i.e., discarding seqs not fully matching (shorter) given
                               region range

Examples

  1. No region given.

     $ echo -ne ">seq\nacgcccactgaaatga\n" 
 >seq
 acgcccactgaaatga

 $ echo -ne ">seq\nacgcccactgaaatga\n" \
     | seqkit amplicon -F ccc -R ttt
 >seq
 cccactgaaa
    
 $ echo -ne ">seq\nacgcccactgaaatga\n" \
     | seqkit amplicon -F ccc -R ttt

 # BED6+1
 $ echo -ne ">seq\nacgcccactgaaatga\n" \
     | seqkit amplicon -F ccc -R ttt --bed
 seq     3       13      .       0       +       cccactgaaa
    
 # supporting degenerate bases.
 $ echo -ne ">seq\nacgcccactgaaatga\n" \
     | seqkit amplicon -F ccR -R ttt --bed
 seq     4       13      .       0       +       ccactgaaa

  2. Output mismatches:

     $ echo -ne ">seq\nacgcccactgaaatga\n" \
     | seqkit amplicon -F ccc -R ttt -M
 >seq mismatches=0(0+0)
 cccactgaaa
    
 $ echo -ne ">seq\nacgcccactgaaatga\n" \
     | seqkit amplicon -F ccc -R ttt -m 1 -M
 >seq mismatches=2(1+1)
 cgcccactgaaat
    
 $ echo -ne ">seq\nacgcccactgaaatga\n" \
     | seqkit amplicon -F ccc -R ttt -m 1 -M --bed
 seq     1       14      .       0       +       cgcccactgaaat   2       1       1

  3. Load primers from 3- or 2-column tabular primer file, with first column as primer name.

     $ cat seqs4amplicon.fa 
 >seq1
 ACGCCCACTGAAATGA
 >seq2
 ACGTACGGTCAGATCCA
    
 $ cat primers.tsv 
 p1      ccc     ttt
 p2      ttt     ccc
 p3      ttt

 p4      CG      TG
 P5      CG      GA

 # containing degenerate bases
 p6      TRC     WGG

    
 $ cat seqs4amplicon.fa | seqkit amplicon -p primers.tsv --bed
 seq1    3       13      p1      0       +       CCCACTGAAA
 seq1    1       7       p4      0       +       CGCCCA
 seq1    3       13      p2      0       -       TTTCAGTGGG
 seq1    10      13      p3      0       -       TTT
 seq2    1       17      p4      0       +       CGTACGGTCAGATCCA
 seq2    1       15      P5      0       +       CGTACGGTCAGATC
 seq2    3       17      p6      0       +       TACGGTCAGATCCA

  4. Inner region

     # region right behind forward primer
 $ echo -ne ">seq\nacgcccactgaaatga\n" \
     | seqkit amplicon -F ccc -R ttt -r 4:7
 >seq
 actg
    
 # BED
 $ echo -ne ">seq\nacgcccactgaaatga\n" \
     | seqkit amplicon -F ccc -R ttt -r 4:7 --bed
 seq     6       10      .       0       +       actg
    
 # more common case is triming primers
 $ echo -ne ">seq\nacgcccactgaaatga\n" \
     | seqkit amplicon -F ccc -R ttt -r 4:-4
 >seq
 actg

  5. flanking region

     # in one of my sequencing data, I only care about 
 # region downstream of forward primer
 $ echo -ne ">seq\nacgcccactgaaatga\n" \
     | seqkit amplicon -F ccc -f -r 3:6
 >seq
 tgaa
    
 $ echo -ne ">seq\nacgcccactgaaatga\n" \
     | seqkit amplicon -F ccc -f -r 3:6 --bed
 seq     8       12      .       0       +       tgaa

    
 # if given region if out scope of sequence. e.g,
 # 2-5bp downstream of aaa, we can get part of region (2-4) by default
 $ echo -ne ">seq\nacgcccactgaaatga\n" \
     | seqkit amplicon -F aaa -f -r 2:5
 >seq
 ga
    
 # you can also use strict mode to discard those cases
 $ echo -ne ">seq\nacgcccactgaaatga\n" \
     | seqkit amplicon -F aaa -f -r 2:5 -s

duplicate

Usage

duplicate sequences N times

You may need "seqkit rename" to make the the sequence IDs unique.

Usage:
  seqkit duplicate [flags]

Aliases:
  duplicate, dup

Flags:
  -h, --help        help for duplicate
  -n, --times int   duplication number (default 1)

Examples

  1. Data

     $ cat tests/hairpin.fa | seqkit head -n 1
 >cel-let-7 MI0000001 Caenorhabditis elegans let-7 stem-loop
 UACACUGUGGAUCCGGUGAGGUAGUAGGUUGUAUAGUUUGGAAUAUUACCACCGGUGAAC
 UAUGCAAUUUUCUACCUUACCGGAGACAGAACUCUUCGA

  2. Duplicate 2 times

     $ cat tests/hairpin.fa | seqkit head -n 1 \
     | seqkit duplicate -n 2
 >cel-let-7 MI0000001 Caenorhabditis elegans let-7 stem-loop
 UACACUGUGGAUCCGGUGAGGUAGUAGGUUGUAUAGUUUGGAAUAUUACCACCGGUGAAC
 UAUGCAAUUUUCUACCUUACCGGAGACAGAACUCUUCGA
 >cel-let-7 MI0000001 Caenorhabditis elegans let-7 stem-loop
 UACACUGUGGAUCCGGUGAGGUAGUAGGUUGUAUAGUUUGGAAUAUUACCACCGGUGAAC
 UAUGCAAUUUUCUACCUUACCGGAGACAGAACUCUUCGA

  3. use seqkit rename to make the the sequence IDs unique.

     $ cat tests/hairpin.fa | seqkit head -n 1 \
     | seqkit duplicate -n 2 | seqkit rename
 >cel-let-7 MI0000001 Caenorhabditis elegans let-7 stem-loop
 UACACUGUGGAUCCGGUGAGGUAGUAGGUUGUAUAGUUUGGAAUAUUACCACCGGUGAAC
 UAUGCAAUUUUCUACCUUACCGGAGACAGAACUCUUCGA
 >cel-let-7_2 MI0000001 Caenorhabditis elegans let-7 stem-loop
 UACACUGUGGAUCCGGUGAGGUAGUAGGUUGUAUAGUUUGGAAUAUUACCACCGGUGAAC
 UAUGCAAUUUUCUACCUUACCGGAGACAGAACUCUUCGA

rmdup

Usage

remove duplicated sequences by ID/name/sequence

Attentions:
  1. When comparing by sequences, both positive and negative strands are
     compared. Switch on -P/--only-positive-strand for considering the
     positive strand only.
  2. Only the first record is saved for duplicates.

Usage:
  seqkit rmdup [flags]

Flags:
  -n, --by-name                by full name instead of just id
  -s, --by-seq                 by seq
  -D, --dup-num-file string    file to save number and list of duplicated seqs
  -d, --dup-seqs-file string   file to save duplicated seqs
  -h, --help                   help for rmdup
  -i, --ignore-case            ignore case
  -P, --only-positive-strand   only considering positive strand when comparing by sequence

Examples

Similar to common.

  1. General use

     $ zcat hairpin.fa.gz | seqkit rmdup -s -o clean.fa.gz
 [INFO] 2226 duplicated records removed

 $ zcat reads_1.fq.gz | seqkit rmdup -s -o clean.fa.gz
 [INFO] 1086 duplicated records removed

  2. Save duplicated sequences to file

     $ zcat hairpin.fa.gz \
     | seqkit rmdup -s -i -o clean.fa.gz -d duplicated.fa.gz -D duplicated.detail.txt

 $ cat duplicated.detail.txt   # here is not the entire list
 3	hsa-mir-424, mml-mir-424, ppy-mir-424
 3	hsa-mir-342, mml-mir-342, ppy-mir-342
 2	ngi-mir-932, nlo-mir-932
 2	ssc-mir-9784-1, ssc-mir-9784-2

common

Usage

find common sequences of multiple files by id/name/sequence

Note:
  1. 'seqkit common' is designed to support 2 and MORE files.
  2. When comparing by sequences, both positive and negative strands are
     compared. Switch on -P/--only-positive-strand for considering the
     positive strand only.
  3. For 2 files, 'seqkit grep' is much faster and consumes lesser memory:
     seqkit grep -f <(seqkit seq -n -i small.fq.gz) big.fq.gz # by seq ID
     seqkit grep -s -f <(seqkit seq -s small.fq.gz) big.fq.gz # by seq
  4. Some records in one file may have same sequences/IDs. They will ALL be
     retrieved if the sequence/ID was shared in multiple files.
     So the records number may be larger than that of the smallest file.

Usage:
  seqkit common [flags]

Flags:
  -n, --by-name                match by full name instead of just id
  -s, --by-seq                 match by sequence
  -h, --help                   help for common
  -i, --ignore-case            ignore case
  -P, --only-positive-strand   only considering positive strand when comparing by sequence

Examples

  1. By ID (default)

     seqkit common file*.fa -o common.fasta

  2. By full name

     seqkit common file*.fa -n -o common.fasta

  3. By sequence

     seqkit common file*.fa -s -i -o common.fasta

split

Usage

split sequences into files by name ID, subsequence of given region,
part size or number of parts.

If you just want to split by parts or sizes, please use "seqkit split2",
which can apply to paired- and single-end FASTQ.

If you want to cut a sequence into multiple segments.
  1. For cutting into even chunks, please use 'kmcp utils split-genomes'
     (https://bioinf.shenwei.me/kmcp/usage/#split-genomes).
     E.g., cutting into 4 segments of equal size, with no overlap between adjacent segments:
        kmcp utils split-genomes -m 1 -k 1 --split-number 4 --split-overlap 0 input.fasta -O out
  2. For cutting into multiple chunks of fixed size, please using 'seqkit sliding'.
     E.g., cutting into segments of 40 bp and keeping the last segment which can be shorter than 40 bp.
        seqkit sliding -g -s 40 -W 40 input.fasta -o out.fasta

Attentions:
  1. For the two-pass mode (-2/--two-pass), The flag -U/--update-faidx is recommended to
     ensure the .fai file matches the FASTA file.

The definition of region is 1-based and with some custom design.

Examples:

 1-based index    1 2 3 4 5 6 7 8 9 10
negative index    0-9-8-7-6-5-4-3-2-1
           seq    A C G T N a c g t n
           1:1    A
           2:4      C G T
         -4:-2                c g t
         -4:-1                c g t n
         -1:-1                      n
          2:-2      C G T N a c g t
          1:-1    A C G T N a c g t n
          1:12    A C G T N a c g t n
        -12:-1    A C G T N a c g t n

Usage:
  seqkit split [flags]

Flags:
  -i, --by-id                     split squences according to sequence ID
      --by-id-prefix string       file prefix for --by-id
  -p, --by-part int               split sequences into N parts
      --by-part-prefix string     file prefix for --by-part
  -r, --by-region string          split squences according to subsequence of given region. e.g 1:12 for
                                  first 12 bases, -12:-1 for last 12 bases. type "seqkit split -h" for
                                  more examples
      --by-region-prefix string   file prefix for --by-region
  -s, --by-size int               split sequences into multi parts with N sequences
      --by-size-prefix string     file prefix for --by-size
  -d, --dry-run                   dry run, just print message and no files will be created.
  -e, --extension string          set output file extension, e.g., ".gz", ".xz", or ".zst"
  -f, --force                     overwrite output directory
  -h, --help                      help for split
  -k, --keep-temp                 keep temporary FASTA and .fai file when using 2-pass mode
  -O, --out-dir string            output directory (default value is $infile.split)
  -2, --two-pass                  two-pass mode read files twice to lower memory usage. (only for FASTA
                                  format)
  -U, --update-faidx              update the fasta index file if it exists. Use this if you are not sure
                                  whether the fasta file changed

Examples

  1. Split sequences into parts with at most 10000 sequences

     $ seqkit split hairpin.fa.gz -s 10000
 [INFO] split into 10000 seqs per file
 [INFO] write 10000 sequences to file: hairpin.fa.part_001.gz
 [INFO] write 10000 sequences to file: hairpin.fa.part_002.gz
 [INFO] write 8645 sequences to file: hairpin.fa.part_003.gz

  2. Split sequences into 4 parts

     $ seqkit split hairpin.fa.gz -p 4
 [INFO] split into 4 parts
 [INFO] read sequences ...
 [INFO] read 28645 sequences
 [INFO] write 7162 sequences to file: hairpin.fa.part_001.gz
 [INFO] write 7162 sequences to file: hairpin.fa.part_002.gz
 [INFO] write 7162 sequences to file: hairpin.fa.part_003.gz
 [INFO] write 7159 sequences to file: hairpin.fa.part_004.gz

    To reduce memory usage when spliting big file, we should alwasy use flag --two-pass

     $ seqkit split hairpin.fa.gz -p 4 -2
 [INFO] split into 4 parts
 [INFO] read and write sequences to temporary file: hairpin.fa.gz.fastx ...
 [INFO] 28645 sequences saved
 [INFO] create or read FASTA index ...
 [INFO] create FASTA index for hairpin.fa.gz.fastx
 [INFO]   28645 records loaded from hairpin.fa.gz.fastx.seqkit.fai
 [INFO] write 7162 sequences to file: hairpin.part_001.fa.gz
 [INFO] write 7162 sequences to file: hairpin.part_002.fa.gz
 [INFO] write 7162 sequences to file: hairpin.part_003.fa.gz
 [INFO] write 7159 sequences to file: hairpin.part_004.fa.gz

  3. Split sequences by species. i.e. by custom IDs (first three letters)

     $ seqkit split hairpin.fa.gz -i --id-regexp "^([\w]+)\-" -2
 [INFO] split by ID. idRegexp: ^([\w]+)\-
 [INFO] read and write sequences to temporary file: hairpin.fa.gz.fastx ...
 [INFO] 28645 sequences saved
 [INFO] create or read FASTA index ...
 [INFO] create FASTA index for hairpin.fa.gz.fastx
 [INFO]   28645 records loaded from hairpin.fa.gz.fastx.seqkit.fai
 [INFO] write 48 sequences to file: hairpin.id_cca.fa.gz
 [INFO] write 3 sequences to file: hairpin.id_hci.fa.gz
 [INFO] write 106 sequences to file: hairpin.id_str.fa.gz
 [INFO] write 1 sequences to file: hairpin.id_bkv.fa.gz
 ...

  4. Split sequences by sequence region (for example, sequence barcode)

     $ seqkit split hairpin.fa.gz -r 1:3 -2
 [INFO] split by region: 1:3
 [INFO] read and write sequences to temporary file: hairpin.fa.gz.fastx ...
 [INFO] 28645 sequences saved
 [INFO] read sequence IDs and sequence region from FASTA file ...
 [INFO] create or read FASTA index ...
 [INFO] create FASTA index for hairpin.fa.gz.fastx
 [INFO]   28645 records loaded from hairpin.fa.gz.fastx.seqkit.fai
 [INFO] write 463 sequences to file: hairpin.region_1:3_AUG.fa.gz
 [INFO] write 349 sequences to file: hairpin.region_1:3_ACU.fa.gz
 [INFO] write 311 sequences to file: hairpin.region_1:3_CGG.fa.gz

    Sequence suffix could be defined as -r -12:-1

split2

Usage

split sequences into files by part size or number of parts

This command supports FASTA and paired- or single-end FASTQ with low memory
occupation and fast speed.

The prefix of output files:
  1. For stdin: stdin
  2. Others: same to the input file
  3. Set via the options: --by-length-prefix, --by-part-prefix, or --by-size-prefix

The extension of output files:
  1. For stdin: .fast[aq]
  2. Others: same to the input file
  3. Additional extension via the option -e/--extension, e.g.， outputting
     gzipped files for plain text input:
         seqkit split2 -p 2 -O test tests/hairpin.fa -e .gz


If you want to cut a sequence into multiple segments.
  1. For cutting into even chunks, please use 'kmcp utils split-genomes'
     (https://bioinf.shenwei.me/kmcp/usage/#split-genomes).
     E.g., cutting into 4 segments of equal size, with no overlap between adjacent segments:
        kmcp utils split-genomes -m 1 -k 1 --split-number 4 --split-overlap 0 input.fasta -O out
  2. For cutting into multiple chunks of fixed size, please using 'seqkit sliding'.
     E.g., cutting into segments of 40 bp and keeping the last segment which can be shorter than 40 bp.
        seqkit sliding -g -s 40 -W 40 input.fasta -o out.fasta

Usage:
  seqkit split2 [flags]

Flags:
  -l, --by-length string          split sequences into chunks of >=N bases, supports K/M/G suffix
      --by-length-prefix string   file prefix for --by-length
  -p, --by-part int               split sequences into N parts
      --by-part-prefix string     file prefix for --by-part
  -s, --by-size int               split sequences into multi parts with N sequences
      --by-size-prefix string     file prefix for --by-size
  -e, --extension string          set output file extension, e.g., ".gz", ".xz", or ".zst"
  -f, --force                     overwrite output directory
  -h, --help                      help for split2
  -O, --out-dir string            output directory (default value is $infile.split)
  -1, --read1 string              (gzipped) read1 file
  -2, --read2 string              (gzipped) read2 file

Examples

  1. Split sequences into parts with at most 10000 sequences:

     $ seqkit split2 hairpin.fa -s 10000
 [INFO] split seqs from hairpin.fa
 [INFO] split into 10000 seqs per file
 [INFO] write 10000 sequences to file: hairpin.fa.split/hairpin.part_001.fa
 [INFO] write 10000 sequences to file: hairpin.fa.split/hairpin.part_002.fa
 [INFO] write 8645 sequences to file: hairpin.fa.split/hairpin.part_003.fa

  2. Force compression for plain text input by adding an extra extension:

     # gzip
 $ seqkit split2 hairpin.fa -O test -f -s 10000 -e .gz
 [INFO] split seqs from hairpin.fa
 [INFO] split into 10000 seqs per file
 [INFO] write 10000 sequences to file: test/hairpin.part_001.fa.gz
 [INFO] write 10000 sequences to file: test/hairpin.part_002.fa.gz
 [INFO] write 8645 sequences to file: test/hairpin.part_003.fa.gz
    
 # xz
 $ seqkit split2 hairpin.fa -O test -f -s 10000 -e .xz
 [INFO] split seqs from hairpin.fa
 [INFO] split into 10000 seqs per file
 [INFO] write 10000 sequences to file: test/hairpin.part_001.fa.xz
 [INFO] write 10000 sequences to file: test/hairpin.part_002.fa.xz
 [INFO] write 8645 sequences to file: test/hairpin.part_003.fa.xz
    
 # zstd
 $ seqkit split2 hairpin.fa -O test -f -s 10000 -e .zst
 [INFO] split seqs from hairpin.fa
 [INFO] split into 10000 seqs per file
 [INFO] write 10000 sequences to file: test/hairpin.part_001.fa.zst
 [INFO] write 10000 sequences to file: test/hairpin.part_002.fa.zst
 [INFO] write 8645 sequences to file: test/hairpin.part_003.fa.zst

  3. Change the prefix of output files:

     $ seqkit split2 hairpin.fa -O test -f -s 10000 -e .gz -o xxx
 [INFO] split seqs from hairpin.fa
 [INFO] split into 10000 seqs per file
 [INFO] write 10000 sequences to file: test/xxx.part_001.fa.gz
 [INFO] write 10000 sequences to file: test/xxx.part_002.fa.gz
 [INFO] write 8645 sequences to file: test/xxx.part_003.fa.gz
    
 # here, we also change the compression format from xz to zstd
 $ cat hairpin.fa.xz | seqkit split2 -O test -f -s 10000 -e .zst
 [INFO] split seqs from stdin
 [INFO] split into 10000 seqs per file
 [INFO] write 10000 sequences to file: test/stdin.part_001.fasta.zst
 [INFO] write 10000 sequences to file: test/stdin.part_002.fasta.zst
 [INFO] write 8645 sequences to file: test/stdin.part_003.fasta.zst

  4. Split sequences into 4 parts

     $ seqkit split hairpin.fa.gz -p 4 -f
 [INFO] split into 4 parts
 [INFO] read sequences ...
 [INFO] read 28645 sequences
 [INFO] write 7162 sequences to file: hairpin.fa.gz.split/hairpin.part_001.fa.gz
 [INFO] write 7162 sequences to file: hairpin.fa.gz.split/hairpin.part_002.fa.gz
 [INFO] write 7162 sequences to file: hairpin.fa.gz.split/hairpin.part_003.fa.gz
 [INFO] write 7159 sequences to file: hairpin.fa.gz.split/hairpin.part_004.fa.gz

  5. For FASTQ files (paired-end)

     $ seqkit split2 -1 reads_1.fq.gz -2 reads_2.fq.gz -p 2 -O out -f
 [INFO] split seqs from reads_1.fq.gz and reads_2.fq.gz
 [INFO] split into 2 parts
 [INFO] write 1250 sequences to file: out/reads_2.part_001.fq.gz
 [INFO] write 1250 sequences to file: out/reads_2.part_002.fq.gz
 [INFO] write 1250 sequences to file: out/reads_1.part_001.fq.gz
 [INFO] write 1250 sequences to file: out/reads_1.part_002.fq.gz

  6. For FASTA files (single-end)

     $ seqkit split2 -1 reads_1.fq.gz reads_2.fq.gz -p 2 -O out -f
 [INFO] flag -1/--read1 given, ignore: reads_2.fq.gz
 [INFO] split seqs from reads_1.fq.gz
 [INFO] split into 2 parts
 [INFO] write 1250 sequences to file: out/reads_1.part_001.fq.gz
 [INFO] write 1250 sequences to file: out/reads_1.part_002.fq.gz

 $ seqkit split2 reads_1.fq.gz -p 2 -O out -f
 [INFO] split seqs from reads_1.fq.gz
 [INFO] split into 2 parts
 [INFO] write 1250 sequences to file: out/reads_1.part_001.fq.gz
 [INFO] write 1250 sequences to file: out/reads_1.part_002.fq.gz

pair

Usage

match up paired-end reads from two fastq files

Attentions:
1. Orders of headers in the two files better be the same (not shuffled),
   otherwise, it consumes a huge number of memory for buffering reads in memory.
2. Unpaired reads are optional outputted with the flag -u/--save-unpaired.
3. If the flag -O/--out-dir is not given, the output will be saved in the same directory
   of input, with the suffix "paired", e.g., read_1.paired.fq.gz.
   Otherwise, names are kept untouched in the given output directory.
4. Paired gzipped files may be slightly larger than original files, because
   of using a different gzip package/library, don't worry.

Tips:
1. Support for '/1 'and '/2' tags for paired read files generated by platforms like MGI.
   You can simply specify the regular expression for extracting sequence IDs:
     --id-regexp '^(\S+)\/[12]'

Usage:
  seqkit pair [flags]

Flags:
  -f, --force            overwrite output directory
  -h, --help             help for pair
  -O, --out-dir string   output directory
  -1, --read1 string     (gzipped) read1 file
  -2, --read2 string     (gzipped) read2 file
  -u, --save-unpaired    save unpaired reads if there are

Examples

  1. Simple one

     $ seqkit pair -1 reads_1.fq.gz -2 reads_2.fq.gz
    
 # output
 reads_1.paired.fq.gz
 reads_2.paired.fq.gz

  2. Set output directory, file names are kept untouched.

     $ seqkit pair -1 reads_1.fq.gz -2 reads_2.fq.gz -O result
    
 $ tree result        
 result/
 ├── reads_1.fq.gz
 └── reads_2.fq.gz

  3. Save unpaired reads if there are.

     $ seqkit pair -1 reads_1.fq.gz -2 reads_2.fq.gz -O result -u
    
 $ tree result
 result
 ├── reads_1.fq.gz
 ├── reads_1.unpaired.fq.gz
 ├── reads_2.fq.gz
 └── reads_2.unpaired.fq.gz

  4. Support for ‘/1 ‘and ‘/2’ tags for paired read files generated by platforms like MGI. You can simply specify the regular expression for extracting sequence IDs:

     seqkit pair --id-regexp '^(\S+)\/[12]' -1 reads_1.fq.gz -2 reads_2.fq.gz

sample

Usage

sample sequences by number or proportion.

Attention:
1. Do not use '-n' on large FASTQ files, it loads all seqs into memory!
   use 'seqkit sample -p 0.1 seqs.fq.gz | seqkit head -n N' instead!

Usage:
  seqkit sample [flags]

Flags:
  -h, --help               help for sample
  -n, --number int         sample by number (result may not exactly match), DO NOT use on large FASTQ files.
  -p, --proportion float   sample by proportion
  -s, --rand-seed int      rand seed (default 11)
  -2, --two-pass           2-pass mode read files twice to lower memory usage. Not allowed when reading
                           from stdin

Examples

  1. Sample by proportion

     $ zcat hairpin.fa.gz | seqkit sample -p 0.1 -o sample.fa.gz
 [INFO] sample by proportion
 [INFO] 2814 sequences outputed

  2. Sample by number

     $ zcat hairpin.fa.gz | seqkit sample -n 1000 -o sample.fa.gz
 [INFO] sample by number
 [INFO] 949 sequences outputed

    949 != 1000 ??? see Effect of random seed on results of seqkit sample

    To reduce memory usage when spliting big file, we could use flag --two-pass

    We can also use seqkit sample -p followed with seqkit head -n:

     $ zcat hairpin.fa.gz \
     | seqkit sample -p 0.1 \
     | seqkit head -n 1000 -o sample.fa.gz

  3. Set rand seed to reproduce the result

     $ zcat hairpin.fa.gz \
     | seqkit sample -p 0.1 -s 11

  4. Most of the time, we could shuffle after sampling

     $ zcat hairpin.fa.gz \
     | seqkit sample -p 0.1 \
     | seqkit shuffle -o sample.fa.gz

Note that when sampling on FASTQ files, make sure using same random seed by flag -s (--rand-seed)

head

Usage

print first N FASTA/Q records

For returning the last N records, use:
    seqkit range -N:-1 seqs.fasta

Usage:
  seqkit head [flags]

Flags:
  -n, --number int   print first N FASTA/Q records (default 10)

Examples

  1. FASTA

     $ seqkit head -n 1 hairpin.fa.gz
 >cel-let-7 MI0000001 Caenorhabditis elegans let-7 stem-loop
 UACACUGUGGAUCCGGUGAGGUAGUAGGUUGUAUAGUUUGGAAUAUUACCACCGGUGAAC
 UAUGCAAUUUUCUACCUUACCGGAGACAGAACUCUUCGA

  2. FASTQ

     $ seqkit head -n 1 reads_1.fq.gz
 @HWI-D00523:240:HF3WGBCXX:1:1101:2574:2226 1:N:0:CTGTAG
 TGAGGAATATTGGTCAATGGGCGCGAGCCTGAACCAGCCAAGTAGCGTGAAGGATGACTGCCCTACGGGTTGTAA
 +
 HIHIIIIIHIIHGHHIHHIIIIIIIIIIIIIIIHHIIIIIHHIHIIIIIGIHIIIIHHHHHHGHIHIIIIIIIII

head-genome

Usage


print sequences of the first genome with common prefixes in name

For a FASTA file containing multiple contigs of strains (see example below),
these's no list of IDs available for retrieving sequences of a certain strain,
while descriptions of each strain share the same prefix.

This command is used to restrieve sequences of the first strain,
i.e., "Vibrio cholerae strain M29".

>NZ_JFGR01000001.1 Vibrio cholerae strain M29 Contig_1, whole genome shotgun sequence
>NZ_JFGR01000002.1 Vibrio cholerae strain M29 Contig_2, whole genome shotgun sequence
>NZ_JFGR01000003.1 Vibrio cholerae strain M29 Contig_3, whole genome shotgun sequence
>NZ_JSTP01000001.1 Vibrio cholerae strain 2012HC-12 NODE_79, whole genome shotgun sequence
>NZ_JSTP01000002.1 Vibrio cholerae strain 2012HC-12 NODE_78, whole genome shotgun sequence

Attention:

  1. Sequences in file should be well organized.

Usage:
  seqkit head-genome [flags]

Flags:
  -h, --help                    help for head-genome
  -m, --mini-common-words int   minimal shared prefix words (default 4)

range

Usage

print FASTA/Q records in a range (start:end)

Examples:
  1. leading 100 records (head -n 100)
      seqkit range -r 1:100
  2. last 100 records (tail -n 100)
      seqkit range -r -100:-1
  3. remove leading 100 records (tail -n +101)
      seqkit range -r 101:-1
  4. other ranges:
      seqkit range -r 10:100
      seqkit range -r -100:-10

Usage:
  seqkit range [flags]

Flags:
  -h, --help           help for range
  -r, --range string   range. e.g., 1:12 for first 12 records (head -n 12), -12:-1 for last 12 records (tail -n 12)

Examples

  1. leading N records (head)

     $ cat tests/hairpin.fa | seqkit head -n 100 | md5sum
 f65116af7d9298d93ba4b3d19077bbf1  -
 $ cat tests/hairpin.fa | seqkit range -r 1:100 | md5sum
 f65116af7d9298d93ba4b3d19077bbf1  -

  2. last N records (tail)

     $ cat tests/hairpin.fa | seqkit range -r -100:-1 | seqkit stats
 file  format  type  num_seqs  sum_len  min_len  avg_len  max_len
 -     FASTA   RNA        100    8,656       58     86.6      172

  3. remove leading 100 records (tail -n +101)

     $ seqkit range -r 101:-1 tests/hairpin.fa | seqkit sum -a
 seqkit.v0.1_RLS_k0_e1feced9bb0be653afa8205dec4987db     -       28545   2940929
    
 $ seqkit fx2tab tests/hairpin.fa | tail -n +101 | seqkit tab2fx | seqkit sum -a
 seqkit.v0.1_RLS_k0_e1feced9bb0be653afa8205dec4987db     -       28545   2940929

  4. Other ranges

     $ cat tests/hairpin.fa | seqkit range -r 101:150 | seqkit stats
 file  format  type  num_seqs  sum_len  min_len  avg_len  max_len
 -     FASTA   RNA         50    3,777       63     75.5       96

 $ cat tests/hairpin.fa | seqkit range -r -100:-2 | seqkit stats
 file  format  type  num_seqs  sum_len  min_len  avg_len  max_len
 -     FASTA   RNA         99    8,484       58     85.7      146

replace

Usage

replace name/sequence by regular expression.

Note that the replacement supports capture variables.
e.g. $1 represents the text of the first submatch.
ATTENTION: use SINGLE quote NOT double quotes in *nix OS.

Examples: Adding space to all bases.

    seqkit replace -p "(.)" -r '$1 ' -s

Or use the \ escape character.

    seqkit replace -p "(.)" -r "\$1 " -s

more on: http://bioinf.shenwei.me/seqkit/usage/#replace

Special replacement symbols (only for replacing name not sequence):

    {nr}    Record number, starting from 1
    {kv}    Corresponding value of the key (captured variable $n) by key-value file,
            n can be specified by flag -I (--key-capt-idx) (default: 1)
            
Special cases:
  1. If replacements contain '$', 
    a). If using '{kv}', you need use '$$$$' instead of a single '$':
            -r '{kv}' -k <(sed 's/\$/$$$$/' kv.txt)
    b). If not, use '$$':
            -r 'xxx$$xx'

Filtering records to edit:
  You can use flags similar to those in "seqkit grep" to choose partly records to edit.

Usage:
  seqkit replace [flags]

Flags:
  -s, --by-seq                   replace seq (only FASTA)
      --f-by-name                [target filter] match by full name instead of just ID
      --f-by-seq                 [target filter] search subseq on seq, both positive and negative strand
                                 are searched, and mismatch allowed using flag -m/--max-mismatch
      --f-ignore-case            [target filter] ignore case
      --f-invert-match           [target filter] invert the sense of matching, to select non-matching records
      --f-only-positive-strand   [target filter] only search on positive strand
      --f-pattern strings        [target filter] search pattern (multiple values supported. Attention:
                                 use double quotation marks for patterns containing comma, e.g., -p
                                 '"A{2,}"')
      --f-pattern-file string    [target filter] pattern file (one record per line)
      --f-use-regexp             [target filter] patterns are regular expression
  -h, --help                     help for replace
  -i, --ignore-case              ignore case
  -K, --keep-key                 keep the key as value when no value found for the key (only for
                                 sequence name)
  -U, --keep-untouch             do not change anything when no value found for the key (only for
                                 sequence name)
  -I, --key-capt-idx int         capture variable index of key (1-based) (default 1)
  -m, --key-miss-repl string     replacement for key with no corresponding value
  -k, --kv-file string           tab-delimited key-value file for replacing key with value when using
                                 "{kv}" in -r (--replacement) (only for sequence name)
      --nr-width int             minimum width for {nr} in flag -r/--replacement. e.g., formatting "1"
                                 to "001" by --nr-width 3 (default 1)
  -p, --pattern string           search regular expression
  -r, --replacement string       replacement. supporting capture variables.  e.g. $1 represents the text
                                 of the first submatch. ATTENTION: for *nix OS, use SINGLE quote NOT
                                 double quotes or use the \ escape character. Record number is also
                                 supported by "{nr}".use ${1} instead of $1 when {kv} given!

Examples

  1. Prepend or append to the header

     $ echo -e ">seq1 abc-123\nACGT-ACGT"
 >seq1 abc-123
 ACGT-ACGT
    
 # prefix
 $ echo -e ">seq1 abc-123\nACGT-ACGT" \
     | seqkit replace -p ^ -r _prefix_
 >_prefix_seq1 abc-123
 ACGT-ACGT
    
 # suffix
 $ echo -e ">seq1 abc-123\nACGT-ACGT" \
     | seqkit replace -p $ -r _suffix_
 >seq1 abc-123_suffix_
 ACGT-ACGT
 $ echo -e ">seq1 abc-123\nACGT-ACGT" \
     | seqkit seq -i \
     | seqkit replace -p $ -r _suffix_
 >seq1_suffix_
 ACGT-ACGT

  2. Remove descriptions

     $ echo -e ">seq1 abc-123\nACGT-ACGT"
 >seq1 abc-123
 ACGT-ACGT
    
 $ echo -e ">seq1 abc-123\nACGT-ACGT" \
     | seqkit replace -p "\s.+"
 >seq1
 ACGT-ACGT

  3. Replace “-“ with “=”

     $ echo -e ">seq1 abc-123\nACGT-ACGT" \
     | seqkit replace -p "\-" -r '='
 >seq1 abc=123
 ACGT-ACGT

  4. Remove gaps in sequences.

     $ echo -e ">seq1 abc-123\nACGT-ACGT" \
     | seqkit replace -p " |-" -s
 >seq1 abc-123
 ACGTACGT

  5. Add space to every base. ATTENTION: use SINGLE quote NOT double quotes in *nix OS

     $ echo -e ">seq1 abc-123\nACGT-ACGT" \
     | seqkit replace -p "(.)" -r '$1 ' -s
 >seq1 abc-123
 A C G T - A C G T
    
 $ echo -e ">seq1 abc-123\nACGT-ACGT" \
     | seqkit replace -p "(.)" -r "\$1 " -s
 >seq1 abc-123
 A C G T - A C G T

  6. Transpose sequence with csvtk

     $ echo -e ">seq1\nACTGACGT\n>seq2\nactgccgt" \
     | seqkit replace -p "(.)" -r     "\$1 " -s \
     | seqkit seq -s -u \
     | csvtk space2tab \
     | csvtk -t transpose
 A       A
 C       C
 T       T
 G       G
 A       C
 C       C
 G       G
 T       T

  7. Rename with number of record

     $ echo -e ">abc\nACTG\n>123\nATTT" \
     |  seqkit replace -p .+ -r "seq_{nr}"
 >seq_1
 ACTG
 >seq_2
 ATTT

 $ echo -e ">abc\nACTG\n>123\nATTT" \
     |  seqkit replace -p .+ -r "seq_{nr}" --nr-width 5
 >seq_00001
 ACTG
 >seq_00002
 ATTT

  8. Replace key with value by key-value file

     $ more test.fa
 >seq1 name1
 CCCCAAAACCCCATGATCATGGATC
 >seq2 name2
 CCCCAAAACCCCATGGCATCATTCA
 >seq3 name3
 CCCCAAAACCCCATGTTGCTACTAG

 $ more alias.txt
 name0   ABC
 name1   123
 name3   Hello
 name4   World

 $ seqkit replace -p ' (.+)$' -r ' {kv}' -k alias.txt test.fa
 [INFO] read key-value file: alias.txt
 [INFO] 4 pairs of key-value loaded
 >seq1 123
 CCCCAAAACCCCATGATCATGGATC
 >seq2
 CCCCAAAACCCCATGGCATCATTCA
 >seq3 Hello
 CCCCAAAACCCCATGTTGCTACTAG

 $ seqkit replace -p ' (.+)$' -r ' {kv}' -k alias.txt test.fa --keep-key
 [INFO] read key-value file: alias.txt
 [INFO] 4 pairs of key-value loaded
 >seq1 123
 CCCCAAAACCCCATGATCATGGATC
 >seq2 name2
 CCCCAAAACCCCATGGCATCATTCA
 >seq3 Hello
 CCCCAAAACCCCATGTTGCTACTAG

    If the target to replace is the sequence ID, meanwhile you would like to keep the description. Please use this:

     # ^(\S+)   # ID
 # (.+)$    # description
 seqkit replace -p '^(\S+)(.+?)$' -r '{kv}$2' -k alias.txt test.fa

  9. convert fasta to genbank style

     $ cat seq.fa
 >seq1
 TTTAAAGAGACCGGCGATTCTAGTGAAATCGAACGGGCAGGTCAATTTCCAACCAGCGAT
 GACGTAATAGATAGATACAAGGAAGTCATTTTTCTTTTAAAGGATAGAAACGGTTAATGC
 TCTTGGGACGGCGCTTTTCTGTGCATAACT
 >seq2
 AAGGATAGAAACGGTTAATGCTCTTGGGACGGCGCTTTTCTGTGCATAACTCGATGAAGC
 CCAGCAATTGCGTGTTTCTCCGGCAGGCAAAAGGTTGTCGAGAACCGGTGTCGAGGCTGT
 TTCCTTCCTGAGCGAAGCCTGGGGATGAACG

 $ cat seq.fa \
     | seqkit replace -s -p '(\w{10})' -r '$1 ' -w 66 \
     | perl -ne 'if (/^>/) {print; $n=1} else {s/ \r?\n$/\n/; printf "%9d %s", $n, $_; $n+=60;}'
 >seq1
         1 TTTAAAGAGA CCGGCGATTC TAGTGAAATC GAACGGGCAG GTCAATTTCC AACCAGCGAT
        61 GACGTAATAG ATAGATACAA GGAAGTCATT TTTCTTTTAA AGGATAGAAA CGGTTAATGC
       121 TCTTGGGACG GCGCTTTTCT GTGCATAACT
 >seq2
         1 AAGGATAGAA ACGGTTAATG CTCTTGGGAC GGCGCTTTTC TGTGCATAAC TCGATGAAGC
        61 CCAGCAATTG CGTGTTTCTC CGGCAGGCAA AAGGTTGTCG AGAACCGGTG TCGAGGCTGT
       121 TTCCTTCCTG AGCGAAGCCT GGGGATGAAC G

  10. only edit some of the records via patterns. Eight flags starting with --f- (with the same usage to these in seqkit grep) can be used to filter records to edit.

     # prepend some text to header
 $ echo -e ">abc\nACTG\n>123\nATTT" \
     | seqkit replace -p ^ -r _
 >_abc
 ACTG
 >_123
 ATTT

 # only edit some records
 $ echo -e ">abc\nACTG\n>123\nATTT" \
     | seqkit replace -p ^ -r _ --f-use-regexp --f-pattern abc
 [INFO] 1 records matched by the filter
 >_abc
 ACTG
 >123
 ATTT
    
 $ echo -e ">abc\nACTG\n>123\nATTT" \
     | seqkit replace -p ^ -r _ --f-use-regexp --f-invert-match --f-pattern abc 
 [INFO] 1 records matched by the filter
 >abc
 ACTG
 >_123
 ATTT

rename

Usage

rename duplicated IDs

Attention:
  1. This command only appends "_N" to duplicated sequence IDs to make them unique.
  2. Use "seqkit replace" for editing sequence IDs/headers using regular expression.

Example:

    $ seqkit seq seqs.fasta 
    >id comment
    actg
    >id description
    ACTG

    $ seqkit rename seqs.fasta
    >id comment
    actg
    >id_2 description
    ACTG

Usage:
  seqkit rename [flags]

Flags:
  -n, --by-name             check duplication by full name instead of just id
  -f, --force               overwrite output directory
  -h, --help                help for rename
  -m, --multiple-outfiles   write results into separated files for multiple input files
  -O, --out-dir string      output directory (default "renamed")
  -1, --rename-1st-rec      rename the first record as well
  -s, --separator string    separator between original ID/name and the counter (default "_")
  -N, --start-num int       starting count number for *duplicated* IDs/names, should be greater than
                            zero (default 2

Examples

$ echo -e ">a comment\nacgt\n>b comment of b\nACTG\n>a comment\naaaa"
>a comment
acgt
>b comment of b
ACTG
>a comment
aaaa

$ echo -e ">a comment\nacgt\n>b comment of b\nACTG\n>a comment\naaaa" \
    | seqkit rename
>a comment
acgt
>b comment of b
ACTG
>a_2 comment
aaaa

$ echo -e ">a comment\nacgt\n>b comment of b\nACTG\n>a comment\naaaa"  \
    | seqkit rename -s '|' -N 2 -1
>a|1 comment
acgt
>b|1 comment of b
ACTG
>a|2 comment
aaaa

restart

Usage

reset start position for circular genome

Examples

    $ echo -e ">seq\nacgtnACGTN"
    >seq
    acgtnACGTN

    $ echo -e ">seq\nacgtnACGTN" | seqkit restart -i 2
    >seq
    cgtnACGTNa

    $ echo -e ">seq\nacgtnACGTN" | seqkit restart -i -2
    >seq
    TNacgtnACG

Usage:
  seqkit restart [flags]

Flags:
  -i, --new-start int   new start position (1-base, supporting negative value counting from the end) (default 1)

concat

Usage

concatenate sequences with same ID from multiple files

Attentions:
   1. By default, only sequences with IDs that appear in all files are outputted.
      use -f/--full to output all sequences.
   2. If there are more than one sequences of the same ID, we output the Cartesian
      product of sequences.
   3. Description are also concatenated with a separator (-s/--separator).
   4. Order of sequences with different IDs are random.

Usage:
  seqkit concat [flags]

Aliases:
  concat, concate

Flags:
  -f, --full               keep all sequences, like full/outer join
  -h, --help               help for concat
  -s, --separator string   separator for descriptions of records with the same ID (default "|")

Examples

$ cat a.fa 
>A 1
a1-
>A 2
a2-
>B 1
b1-

$ cat b.fa 
>A x
ax-
>C 1
c1-

$ seqkit concat a.fa b.fa 
>A 1|x
a1-ax-
>A 2|x
a2-ax-

$ seqkit concat a.fa b.fa --full
>C 1
c1-
>A 1|x
a1-ax-
>A 2|x
a2-ax-
>B 1
b1-

mutate

Usage

edit sequence (point mutation, insertion, deletion)

Attentions:

  1. Mutiple point mutations (-p/--point) are allowed, but only single 
     insertion (-i/--insertion) OR single deletion (-d/--deletion) is allowed.
  2. Point mutation takes place before insertion/deletion.

Notes:

  1. You can choose certain sequences to edit using similar flags in
     'seqkit grep'.

The definition of position is 1-based and with some custom design.

Examples:

 1-based index    1 2 3 4 5 6 7 8 9 10
negative index    0-9-8-7-6-5-4-3-2-1
           seq    A C G T N a c g t n
           1:1    A
           2:4      C G T
         -4:-2                c g t
         -4:-1                c g t n
         -1:-1                      n
          2:-2      C G T N a c g t
          1:-1    A C G T N a c g t n
          1:12    A C G T N a c g t n
        -12:-1    A C G T N a c g t n

Usage:
  seqkit mutate [flags]

Flags:
  -n, --by-name               [match seqs to mutate] match by full name instead of just id
  -d, --deletion string       deletion mutation: deleting subsequence in a range. e.g., -d 1:2 for
                              deleting leading two bases, -d -3:-1 for removing last 3 bases
  -h, --help                  help for mutate
  -I, --ignore-case           [match seqs to mutate] ignore case of search pattern
  -i, --insertion string      insertion mutation: inserting bases behind of given position, e.g., -i
                              0:ACGT for inserting ACGT at the beginning, -1:* for add * to the end
  -v, --invert-match          [match seqs to mutate] invert the sense of matching, to select
                              non-matching records
  -s, --pattern strings       [match seqs to mutate] search pattern (multiple values supported.
                              Attention: use double quotation marks for patterns containing comma, e.g.,
                              -p '"A{2,}"'))
  -f, --pattern-file string   [match seqs to mutate] pattern file (one record per line)
  -p, --point strings         point mutation: changing base at given position. e.g., -p 2:C for setting
                              2nd base as C, -p -1:A for change last base as A
  -r, --use-regexp            [match seqs to mutate] search patterns are regular expression

Examples:

  1. Point mutation:

     $ echo -ne ">1\nACTGNactgn\n>2\nactgnACTGN\n"
 >1
 ACTGNactgn
 >2
 actgnACTGN

 # first base
 $ echo -ne ">1\nACTGNactgn\n>2\nactgnACTGN\n" \
     | seqkit mutate -p 1:x
 [INFO] edit seq: 1
 [INFO] edit seq: 2
 >1
 xCTGNactgn
 >2
 xctgnACTGN

 # 5th base
 $ echo -ne ">1\nACTGNactgn\n>2\nactgnACTGN\n" \
     | seqkit mutate -p 5:x --quiet
 >1
 ACTGxactgn
 >2
 actgxACTGN

 # last base
 $ echo -ne ">1\nACTGNactgn\n>2\nactgnACTGN\n" \
     | seqkit mutate -p -1:x --quiet
 >1
 ACTGNactgx
 >2
 actgnACTGx

 # mutiple locations:

 $ echo -ne ">1\nACTGNactgn\n>2\nactgnACTGN\n" \
     | seqkit mutate -p 1:x -p -1:x --quiet
 >1
 xCTGNactgx
 >2
 xctgnACTGx

  2. Deletion

     # first base
 $ echo -ne ">1\nACTGNactgn\n>2\nactgnACTGN\n" \
     | seqkit mutate -d 1:1 --quiet
 >1
 CTGNactgn
 >2
 ctgnACTGN

 # last 3 bases
 $ echo -ne ">1\nACTGNactgn\n>2\nactgnACTGN\n" \
     | seqkit mutate -d -3:-1 --quiet
 >1
 ACTGNac
 >2
 actgnAC

  3. Insertion: inserting bases behind of given position

     # at the beginning
 $ echo -ne ">1\nACTGNactgn\n>2\nactgnACTGN\n" \
     | seqkit mutate -i 0:xx --quiet
 >1
 xxACTGNactgn
 >2
 xxactgnACTGN

 # at the end
 $ echo -ne ">1\nACTGNactgn\n>2\nactgnACTGN\n" \
     | seqkit mutate -i -1:xx --quiet
 >1
 ACTGNactgnxx
 >2
 actgnACTGNxx

 # behind of 5th base
 $ echo -ne ">1\nACTGNactgn\n>2\nactgnACTGN\n" \
     | seqkit mutate -i 5:x --quiet
 >1
 ACTGNxactgn
 >2
 actgnxACTGN

  4. Choosing which sequences to edit, using similar flags in seqkit grep.

     $ cat tests/hsa.fa
 >chr1 1th seq
 ACTGNactgn
 >chr2 2nd seq
 actgnACTGN
 >chr11 11th seq
 ACTGNACTGN
 >MT mitochondrial seq
 actgnactgn

 # only edit chr1 and chr2
 # or cat tests/hsa.fa | seqkit mutate -p -1:X -s chr1 -s chr2
 $ cat tests/hsa.fa \
     | seqkit mutate -p -1:X -s chr1,chr2
 [INFO] edit seq: chr1 1th seq
 [INFO] edit seq: chr2 2nd seq
 >chr1 1th seq
 ACTGNactgX
 >chr2 2nd seq
 actgnACTGX
 >chr11 11th seq
 ACTGNACTGN
 >MT mitochondrial seq
 actgnactgn

 # using regular expression to match.
 # e,g., editing all chrosomes:
 $ cat tests/hsa.fa \
     | seqkit mutate -p -1:X -r -s chr
 [INFO] edit seq: chr1 1th seq
 [INFO] edit seq: chr2 2nd seq
 [INFO] edit seq: chr11 11th seq
 >chr1 1th seq
 ACTGNactgX
 >chr2 2nd seq
 actgnACTGX
 >chr11 11th seq
 ACTGNACTGX
 >MT mitochondrial seq
 actgnactgn

 # excluding seqs
 $ cat tests/hsa.fa \
     | seqkit mutate -p -1:X -s chr1 -s chr2 -v 
 [INFO] edit seq: chr11 11th seq
 [INFO] edit seq: MT mitochondrial seq
 >chr1 1th seq
 ACTGNactgn
 >chr2 2nd seq
 actgnACTGN
 >chr11 11th seq
 ACTGNACTGX
 >MT mitochondrial seq
 actgnactgX

shuffle

Usage

shuffle sequences.

By default, all records will be readed into memory.
For FASTA format, use flag -2 (--two-pass) to reduce memory usage. FASTQ not
supported.

Firstly, seqkit reads the sequence IDs. If the file is not plain FASTA file,
seqkit will write the sequences to temporary files, and create FASTA index.

Secondly, seqkit shuffles sequence IDs and extract sequences by FASTA index.

Attentions:
  1. For the two-pass mode (-2/--two-pass), The flag -U/--update-faidx is recommended to
     ensure the .fai file matches the FASTA file.

Usage:
  seqkit shuffle [flags]

Flags:
  -h, --help            help for shuffle
  -k, --keep-temp       keep temporary FASTA and .fai file when using 2-pass mode
  -s, --rand-seed int   rand seed for shuffle (default 23)
  -2, --two-pass        two-pass mode read files twice to lower memory usage. (only for FASTA format)
  -U, --update-faidx    update the fasta index file if it exists. Use this if you are not sure whether
                        the fasta file changed

Examples

  1. General use.

     $ seqkit shuffle hairpin.fa.gz > shuffled.fa
 [INFO] read sequences ...
 [INFO] 28645 sequences loaded
 [INFO] shuffle ...
 [INFO] output ...

  2. For big genome, you’d better use two-pass mode so seqkit could use FASTA index to reduce memory usage

     $ time seqkit shuffle -2 hsa.fa > shuffle.fa
 [INFO] create and read FASTA index ...
 [INFO] create FASTA index for hsa.fa
 [INFO] read sequence IDs from FASTA index ...
 [INFO] 194 sequences loaded
 [INFO] shuffle ...
 [INFO] output ...

 real    0m35.080s
 user    0m45.521s
 sys     0m3.411s

Note that when sampling on FASTQ files, make sure using same random seed by flag -s (--rand-seed) for read 1 and 2 files.

sort

Usage

sort sequences by id/name/sequence/length.

By default, all records will be readed into memory.
For FASTA format, use flag -2 (--two-pass) to reduce memory usage. FASTQ not
supported.

Firstly, seqkit reads the sequence head and length information.
If the file is not plain FASTA file,
seqkit will write the sequences to temporary files, and create FASTA index.

Secondly, seqkit sorts sequence by head and length information
and extracts sequences by FASTA index.

Attentions:
  1. For the two-pass mode (-2/--two-pass), The flag -U/--update-faidx is recommended to
     ensure the .fai file matches the FASTA file.

Usage:
  seqkit sort [flags]

Flags:
  -b, --by-bases                by non-gap bases
  -l, --by-length               by sequence length
  -n, --by-name                 by full name instead of just id
  -s, --by-seq                  by sequence
  -G, --gap-letters string      gap letters (default "- \t.")
  -h, --help                    help for sort
  -i, --ignore-case             ignore case
  -k, --keep-temp               keep temporary FASTA and .fai file when using 2-pass mode
  -N, --natural-order           sort in natural order, when sorting by IDs/full name
  -r, --reverse                 reverse the result
  -L, --seq-prefix-length int   length of sequence prefix on which seqkit sorts by sequences (0 for
                                whole sequence) (default 10000)
  -2, --two-pass                two-pass mode read files twice to lower memory usage. (only for FASTA format)
  -U, --update-faidx            update the fasta index file if it exists. Use this if you are not sure
                                whether the fasta file changed

Examples

For FASTA format, use flag -2 (–two-pass) to reduce memory usage

  1. sort by ID

     $ echo -e ">seq1\nACGTNcccc\n>SEQ2\nacgtnAAAA" \
     | seqkit sort --quiet
 >SEQ2
 acgtnAAAA
 >seq1
 ACGTNcccc

  2. sort by ID and in natural order

     $ echo -e ">3\na\n>1\na\n>Y\na\n>x\na\n>Mt\na\n>11\na\n>2\na\n" \
     | seqkit seq -n -i
 3
 1
 Y
 x
 Mt
 11
 2

 $ echo -e ">3\na\n>1\na\n>Y\na\n>x\na\n>Mt\na\n>11\na\n>2\na\n" \
     | seqkit sort -N -i -2 \
     | seqkit seq -n -i
 1
 2
 3
 11
 Mt
 x
 Y

  3. sort by ID, ignoring case.

     $ echo -e ">seq1\nACGTNcccc\n>SEQ2\nacgtnAAAA" \
     | seqkit sort --quiet -i
 >seq1
 ACGTNcccc
 >SEQ2
 acgtnAAAA

  4. sort by seq, ignoring case.

     $ echo -e ">seq1\nACGTNcccc\n>SEQ2\nacgtnAAAA" \
     | seqkit sort --quiet -s -i
 >SEQ2
 acgtnAAAA
 >seq1
 ACGTNcccc

  5. sort by sequence length

     $ echo -e ">seq1\nACGTNcccc\n>SEQ2\nacgtnAAAAnnn\n>seq3\nacgt" \
     | seqkit sort --quiet -l
 >seq3
 acgt
 >seq1
 ACGTNcccc
 >SEQ2
 acgtnAAAAnnn

bam

monitoring and online histograms of BAM record features

Usage:
  seqkit bam [flags]

Flags:
  -B, --bins int             number of histogram bins (default -1)
  -N, --bundle int           partition BAM file into loci (-1) or bundles with this minimum size
  -c, --count string         count reads per reference and save to this file
  -W, --delay int            sleep this many seconds after plotting (default 1)
  -y, --dump                 print histogram data to stderr instead of plotting
  -G, --exclude-ids string   exclude records with IDs contained in this file
  -e, --exec-after string    execute command after reporting
  -E, --exec-before string   execute command before reporting
  -f, --field string         target fields
  -g, --grep-ids string      only keep records with IDs contained in this file
  -h, --help                 help for bam
  -C, --idx-count            fast read per reference counting based on the BAM index
  -i, --idx-stat             fast statistics based on the BAM index
  -O, --img string           save histogram to this PDF/image file
  -H, --list-fields          list all available BAM record features
  -L, --log                  log10(x+1) transform numeric values
  -q, --map-qual int         minimum mapping quality
  -x, --pass                 passthrough mode (forward filtered BAM to output)
  -k, --pretty               pretty print certain TSV outputs
  -F, --prim-only            filter out non-primary alignment records
  -p, --print-freq int       print/report after this many records (-1 for print after EOF) (default -1)
  -Q, --quiet-mode           supress all plotting to stderr
  -M, --range-max float      discard record with field (-f) value greater than this flag (default NaN)
  -m, --range-min float      discard record with field (-f) value less than this flag (default NaN)
  -R, --reset                reset histogram after every report
  -Z, --silent-mode          supress TSV output to stderr
  -s, --stat                 print BAM satistics of the input files
  -T, --tool string          invoke toolbox in YAML format (see documentation)
  -@, --top-bam string       save the top -? records to this bam file
  -?, --top-size int         size of the top-mode buffer (default 100)

Examples

  1. Get detailed statistics from multiple BAM files.

     seqkit bam -s *.bam

  2. Get rough statistics from multiple indexed BAM files.

     seqkit bam -i *.bam

  3. Count reads mapped to references from a BAM stream.

     cat sample.bam | seqkit bam -c counts.tsv  -

  4. Count reads mapped to references using the BAM index.

     seqkit bam -C sorted_indexed.bam

  5. Monitor alignment accuracy from a bam stream and report after every 1000 records, use 20 bins.

     cat sample.bam | seqkit bam -B -f Acc -p 1000 -

  6. Dump selected fields to TSV.

     seqkit bam -f Ref,Acc,RefCov,Strand sample.bam

  7. Save the best 100 records in terms of alignment accuracy to a BAM file.

     seqkit bam -f Acc -@ top_acc_100.bam -? 100 -Q sample.bam

  8. Inkvoke the BAM toolbox.

The BAM toolbox is a collection of filters acting on a stream of BAM records, configured via YAML. The currently available tools can be listed by seqkit bam -T help:

Tool        Description
----        -----------
AccStats    calculates mean accuracy weighted by aligment lengths
AlnContext  filter records by the sequence context at start and end
Dump        dump various record properties in TSV format
help        list all tools with description

Example YAML configs:

Invoking the AccStats tool directly from the command line or YAML config:

seqkit bam -T '{AccStats: {Tsv: "-"}, Sink: True}' input.bam
seqkit bam -T '{Yaml: "tests/examples/bam_tool_acc_stats.yml"}' input.bam

Where the contents of bam_tool_acc_stats.yml are:

AccStats:
  Tsv: "-"
Sink: True

Invoking the AlnContext tool using YAML:

AlnContext:
  Tsv: "-"
  Ref: "../SIRV_150601a.fasta"
  LeftShift: -10
  RightShift: 10
  RegexStart: "T{4,}"
  RegexEnd: "A{4,}"
  Stranded: True
  Invert: True
Sink: True

Invoking the Dump tool using YAML:

Dump:
  Tsv: "-"
  Fields: ["Read", "Ref", "Pos", "EndPos", "MapQual", "Acc", "Match", "Mismatch", "Ins", "Del", "AlnLen", "  ReadLen", "RefLen", "RefAln", "RefCov", "ReadAln", "ReadCov", "Strand", "MeanQual", "LeftClip", "RightClip", "Flags", "IsSec", "  IsSup", "ReadSeq", "ReadAlnSeq", "LeftSoftClipSeq", "RightSoftClip", "LeftHardClip", "RightHardClip"]
Sink: True

The tools can be chained together, for example the YAML using all three tools look like:

AlnContext:
  Tsv: "context.tsv"
  Ref: "../SIRV_150601a.fasta"
  LeftShift: -10
  RightShift: 10
  RegexStart: "T{4,}"
  RegexEnd: "A{4,}"
  Stranded: True
  Invert: True
Dump:
  Tsv: "dump.tsv"
  Fields: ["Read", "Ref", "Pos", "EndPos", "MapQual", "Acc", "Match", "Mismatch", "Ins", "Del", "AlnLen", "  ReadLen", "RefLen", "RefAln", "RefCov", "ReadAln", "ReadCov", "Strand", "MeanQual", "LeftClip", "RightClip", "Flags", "IsSec", "  IsSup", "ReadSeq", "ReadAlnSeq", "LeftSoftClipSeq", "RightSoftClip", "LeftHardClip", "RightHardClip"]
AccStats:
  Tsv: "-"

If the “Sink” parameter is not specified in the last pipeline step, the output BAM records are streamed to the standard output and can be piped into standard tools, for example:

seqkit bam -T '{Yaml: "bam_tool_pipeline.yml"}' ../pcs109_5k_spliced.bam | samtools flagstat -

merge-slides

Usage

merge sliding windows generated from seqkit sliding

For example,

    ref.contig00001_sliding:454531-454680
    ref.contig00001_sliding:454561-454710
    ref.contig00001_sliding:454591-454740
    ref.contig00002_sliding:362281-362430
    ref.contig00002_sliding:362311-362460
    ref.contig00002_sliding:362341-362490
    ref.contig00002_sliding:495991-496140
    ref.contig00044_sliding:1-150
    ref.contig00044_sliding:31-180
    ref.contig00044_sliding:61-210
    ref.contig00044_sliding:91-240

could be merged into

    ref.contig00001 454530  454740
    ref.contig00002 362280  362490
    ref.contig00002 495990  496140
    ref.contig00044 0       240

Output (BED3 format):
    1. chrom      - chromosome name
    2. chromStart - starting position (0-based)
    3. chromEnd   - ending position (0-based)

Usage:
  seqkit merge-slides [flags]

Flags:
  -b, --buffer-size string            size of buffer, supported unit: K, M, G. You need increase the
                                      value when "bufio.Scanner: token too long" error reported (default
                                      "1G")
  -p, --comment-line-prefix strings   comment line prefix (default [#,//])
  -h, --help                          help for merge-slides
  -g, --max-gap int                   maximum distance of starting positions of two adjacent regions, 0
                                      for no limitation, 1 for no merging.
  -l, --min-overlap int               minimum overlap of two adjacent regions, recommend
                                      $sliding_step_size - 1. (default 1)
  -r, --regexp string                 regular expression for extract the reference name and window
                                      position. (default "^(.+)_sliding:(\\d+)\\-(\\d+)")

Example:

seqkit merge-slides sliding_windows.tsv -l 50 -o sliding_windows.merged.tsv

genautocomplete

Usage

generate shell autocompletion script

Supported shell: bash|zsh|fish|powershell

Bash:

    # generate completion shell
    seqkit genautocomplete --shell bash

    # configure if never did.
    # install bash-completion if the "complete" command is not found.
    echo "for bcfile in ~/.bash_completion.d/* ; do source \$bcfile; done" >> ~/.bash_completion
    echo "source ~/.bash_completion" >> ~/.bashrc

Zsh:

    # generate completion shell
    seqkit genautocomplete --shell zsh --file ~/.zfunc/_seqkit

    # configure if never did
    echo 'fpath=( ~/.zfunc "${fpath[@]}" )' >> ~/.zshrc
    echo "autoload -U compinit; compinit" >> ~/.zshrc

fish:

    seqkit genautocomplete --shell fish --file ~/.config/fish/completions/seqkit.fish

Usage:
  seqkit genautocomplete [flags]

Flags:
      --file string   autocompletion file (default "/home/shenwei/.bash_completion.d/seqkit.sh")
  -h, --help          help for genautocomplete
      --type string   autocompletion type (currently only bash supported) (default "bash")