seqkit.md

SEQKIT

Summary Statistics

# reports per fasta entry the name, the sequence, the length and %GC
seqkit fx2tab -lg [fasta]

# reports per fasta entry the name, length and %GC, but not the sequence
seqkit fx2tab -nlg [fasta]

# get number of contigs, total length, min/max length, N50
seqkit stats -a [fasta]

Visual aid

# color FASTA files
seqkit seq --color [fasta] | less -R

Selecting sequences

# extract a single entry
seqkit grep -p [regexp] [fasta]

# delete a single entry
seqkit grep -v -p [regexp] [fasta]

# extract using a list of headers
seqkit grep -f headers.list [fasta]
# NOTE: this does not print the out FASTA in the same order as in headers.list
# its in the order as they appear in the FASTA

# extract using a list of regexps
seqkit grep -rf patterns.list [fasta]

# extract using a list of headers - in order
# -rp if patterns are regexp
# -p  if patterns are not regexp
for header in $(cat headers.list); do seqkit grep -rp $header [fasta]; done

# select sequences that are at least x long
seqkit seq -m 30 [fasta]

# select sequences that are in between these lengths
seqkit seq -m 70 -M 150 [fasta]

Removing duplicate sequences

# remove all duplicates from a fasta
seqkit rmdup -s [fasta]

# remove all duplicates but only if their header is also identical
seqkit fx2tab [fasta] | sort -u | seqkit tab2fx > [deduped.fasta]

Sorting FASTA files

# sort from largest to smallest sequences
seqkit sort -l -r [fasta]

Select parts of sequences

# report sequence for 'ergo_tig00000012' but only for region 16000-17000
# coordinates are 1-indexed


# select the last 20000 bases of a contig
seqkit subseq --chr ergo_tig00000012 --region -20000:-1 consensus.fasta

# extract sequences of particular regions using a BED file
seqkit subseq --bed [bed] [fasta]

Splitting FASTA files

Split a single multiFASTA file into many singleFASTA files

seqkit split -i [fasta]

# split it up into 10 roughly equal parts
seqkit split -p 10 -o output_dir input.fasta

Search for motifs

By default searches the given strand and its reverse complement

seqkit locate -p AFLEADRTGQA [fasta]

# using regexp
seqkit locate -rp AFLEADRTGQA.*VAPARS [fasta]

# if you use a comma in your pattern to denote a range of counts,
# {25,45} - meaning 25 or 26 or 27 ... or 45 times
# you have to wrap it around double quotes and wrap that around single quotes
seqkit locate -rp '"C.{3}C.{25,45}C.{3}C"' [fasta]

Formatting FASTA files

# make FASTA single line
seqkit seq -w 0 [fasta]

# make FASTA multi line
## good if you've edited a multiline fasta file and
## the linewidths are not equal anymore for some reason
seqkit seq -w 60 [fasta]

# only print sequence, not headers
seqkit seq -s [fasta]

# force upper case
seqkit seq -u [fasta]

# reverse complement sequences in a FASTA file
seqkit seq -rp [fasta]

Index the FASTA file

seqkit faidx [fasta]

Translate DNA sequences

# translate over all six frames
seqkit translate -f 6 [fasta]
cat <fasta> | seqkit translate -f 6

# append frame information in sequence identifier
seqkit translate -F -f6 [fasta]
# >CU914152.1_frame=1
# >CU914152.1_frame=2
# >CU914152.1_frame=3
# >CU914152.1_frame=-1
# >CU914152.1_frame=-2
# >CU914152.1_frame=-3