Supported file formats

seqmagick works with Biostrings objects and provides readers/writers for a set of widely used sequence file formats. Compressed files (gzip, bzip2 and xz) are supported transparently — just use them like plain text files.

Format	Read	Write	Related functions
FASTA (.fas/.fa/.fasta)	`fa_read()`	`fa_write()`	`fa_to_interleaved()`, `fa_to_sequential()`, `fa_combine()`, `fas2phy()`, `phy2fas()`
PHYLIP (.phy)	`phy_read()`	`phy_write()`	`fas2phy()`, `phy2fas()`
CLUSTAL (.clw)	`clw_read()`	-	-
STOCKHOLM (.sth)	`sth_read()`	-	-
MEGA (.meg, .mega)	`mega_read()`	-	-
GenBank (.gb)	`gb_read()`	-	`download_genbank()`
NCBI (online)	`ncbi_fa_read()`	-	`download_genbank()`
BAM (.bam)	`bam2DNAStringSet()`	-	-

FASTA and PHYLIP support both sequential and interleaved layouts; specify via the type parameter where applicable.

Sequence I/O

FASTA

fa_file <- system.file("extdata/HA.fas", package="seqmagick")
x <- fa_read(fa_file)
x

## DNAStringSet object of length 100:
##       width sequence                                        names               
##   [1]  1737 GGATAATTCTATTAACCATGAA...TGAGTGCATTAATTAAAAACAC CY168567 A/Boston...
##   [2]  1734 AAAAGCAACAAAAATGAAGGCA...AACATTAGGATTTCAGAAGCAT CY053309 A/Browns...
##   [3]  1762 AGCAAAAGCAGGGGATAATTCT...TAAAAACACCCTTGTTTCTACT LC111587 A/Fukuok...
##   [4]  1723 GGGAAAATAAAAACAACCGAAA...CAGTGCAGAATATGTATTTAAA CY084969 A/swine/...
##   [5]  1701 ATGAAGACTATCATTGCTTTGA...TAGGTGCAACATTTGCATTTGA EU885538 A/New Yo...
##   ...   ... ...
##  [96]  1698 ATGAAAGTAAAACTACTGACCC...GCAATGTAGAATATGCATCTGA KC013577 A/swine/...
##  [97]  1701 ATGAAGACTATCATTGCTTTGA...TAGGTGCAACATTTGCATTTGA GQ385825 A/New Ha...
##  [98]  1710 ATGAAGACTATCATTGCTTTGA...CATTTGCATTTGAGTGCATTAA CY043768 A/Hong K...
##  [99]  1701 ATGAAGGCAATACTAGTAGTTC...ACAGTGTAGAATATGTATTTAA GU371272 A/Hunan ...
## [100]  1701 ATGAAGACTATCATTGCTTTGA...TAGGTGCAACATTTGCATTTGA KM208509 A/Hangzh...

Compressed files can be written and read directly:

tmpgz <- tempfile(fileext=".fas.gz")
fa_write(x[1:5], tmpgz)          ## gzip is used automatically
y <- fa_read(tmpgz)
identical(width(y), width(x[1:5]))

## [1] TRUE

PHYLIP

phy_file <- system.file("extdata/HA.phy", package="seqmagick")
p <- phy_read(phy_file)

CLUSTAL and STOCKHOLM

x <- clw_read("alignment.clw")
x <- sth_read("alignment.sth")

MEGA

meg <- system.file("extdata/mega/Crab_rRNA.meg", package="seqmagick")
m <- mega_read(meg)
m

## DNAStringSet object of length 13:
##      width sequence                                         names               
##  [1]   421 TAATTAAAGGGCCGTGGTATA-C...AAAGAGTTTGCGAGCCTCGATG Artemia_salina
##  [2]   421 GTGTTGAAGGGCCGCGGTATTTT...AAAAAGATTGCGA-CCTCGATG Clibanarius_vittatus
##  [3]   421 ATATTGAAGGGCCGCGGTATTTC...GAAAAGTTTGCGA-CCTCGATG Coenobita_sp.
##  [4]   421 AAATTAAAGAGCCGCAGTATTT-...AAAGAGTTTGCGAGCCTCGATG Lithodes_aequispina
##  [5]   421 AAATTAAAGAGCCGCAGTATTT-...AAAGAGTTTGCGAGCCTCGATG Paralithodes_camt...
##  ...   ... ...
##  [9]   421 GTATTAAAGAGCCGCAGTATTC-...AAAGAGTTTGCGAGCCTCGATG Pagurus_acadianus
## [10]   421 TGATTAAAGAGCCGCAGTATTT-...AAAGAGTTTGCGAGCCTCGATG Pagurus_pollicari...
## [11]   421 TGATTAAAGAGCCGCAGTATTT-...AAAGAGTTTGCGAGCCTCGATG Pagurus_pollicari...
## [12]   421 TAATTAAAGAGCCGCAGTATTT-...AAAGAGTTTGCGAGCCTCGATG Pagurus_longicarp...
## [13]   421 TAATTAAAGAGCCGCAGTATTT-...AAAGAGTTTGCGAGCCTCGATG Pagurus_longicarp...

GenBank

gb <- system.file("extdata/AB115403.gb", package="seqmagick")
x <- gb_read(gb)

Download sequences from NCBI

Both download_genbank() (writes to files) and ncbi_fa_read() (returns Biostrings objects directly) fetch records by accession number.

## save to local files
download_genbank(acc='AB115403', format='genbank', outfile='out.gb')
download_genbank(acc='AB115403', format='fasta', outfile='out.fa')

## or read into R directly
x <- ncbi_fa_read(acc=c('AB115403', 'CY084969'))

Format conversion

fasta and phylip conversion

Note that PHYLIP is only meaningful for aligned sequences. We first subset the FASTA file, align it with bs_aln(), then convert it to PHYLIP:

## use a small subset to keep this demo fast
fa2 <- tempfile(fileext = '.fa')
fa_read(fa_file) |> bs_filter('ATGAAAGTAAAA', by='sequence') |> fa_write(fa2, type='interleaved')

alnfas <- tempfile(fileext = ".fas")
fa_read(fa2) |> bs_aln(quiet=TRUE) |> fa_write(alnfas)

tmpphy <- tempfile(fileext = ".phy")
fas2phy(alnfas, tmpphy, type = 'sequential')

Converting back from PHYLIP to FASTA is symmetric:

phy2fas(tmpphy, alnfas, type = 'interleaved')

interleaved and sequential format conversion

Use the type parameter in fa_write()/phy_write(), or the shortcuts fa_to_interleaved()/fa_to_sequential():

## inter-convert by read + write with different 'type'
fa_read(fa2) |> fa_write("out.fas", type="sequential")
phy_read(tmpphy) |> phy_write("out.phy", type="interleaved")

## or one-liner conversions for FASTA files
fa_to_interleaved("in.fas", "interleaved.fas")
fa_to_sequential("in.fas", "sequential.fas")

Sequence manipulation

Filtering, alignment and consensus

bs <- fa_read(fa_file)

## keep only sequences containing the pattern
f <- bs_filter(bs, 'ATGAAAGTAAAA', by='sequence')

## multiple sequence alignment requires the 'muscle' package
aln <- f |> bs_aln(quiet=TRUE)

## consensus sequence of the alignment
bs_consensus(aln)

Renaming sequences

Rename sequences by supplying a two-column mapping (old name -> new name). Unmatched names are kept unchanged.

mapping <- data.frame(old = c(names(x)[1], names(x)[2]),
                      new = c("HA_HK", "HA_BR"))
z <- bs_rename(x, mapping)
names(z)[1:4]

## [1] "HA_HK"                                                
## [2] "HA_BR"                                                
## [3] "LC111587 A/Fukuoka/DS-261v/2014 2014/01/21 4 (HA)"    
## [4] "CY084969 A/swine/Hong Kong/158/1993 1993/07/19 4 (HA)"

With sep and position, only a specific token of each name (defined by splitting the name with sep) is replaced by the mapped value:

m2 <- data.frame(old = c("seq1", "seq2"),
                 new = c("HA_human", "HA_swine"))
demo <- Biostrings::BStringSet(c("seq1 HK01"="AAAA",
                                 "seq2 TW02"="CCCC",
                                 "seq3 JP03"="GGGG"))
z2 <- bs_rename(demo, m2, sep=" ", position=1)
names(z2)

## [1] "HA_human HK01" "HA_swine TW02" "seq3 JP03"

For files, fa_rename() reads a FASTA file plus a two-column table and writes renamed output in one step:

fa_rename("seqs.fas", "map.txt", outfile="renamed.fas")

Summary statistics

The new fa_summary() gives a quick overview of a FASTA file (or an existing Biostrings object): number of sequences, length range, GC content and proportion of ambiguous characters.

fa_summary(fa_file)

##   n_seq min_len max_len mean_len GC_pct ambiguous_pct
## 1   100    1689    1777  1716.15  41.36          0.01

Works on in-memory objects as well:

fa_summary(x[1:10])

##   n_seq min_len max_len mean_len GC_pct ambiguous_pct
## 1    10    1698    1762   1716.5   41.3          0.01

For per-sequence lengths (gap characters excluded), use the classic seqlen():

head(seqlen(fa_file))

##      [,1]
## [1,] 1737
## [2,] 1734
## [3,] 1762
## [4,] 1723
## [5,] 1701
## [6,] 1701

Bugs/Feature requests

If you have any, let me know. Thx!

Session info

Here is the output of sessionInfo() on the system on which this document was compiled:

## R version 4.6.1 (2026-06-24)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 26.04.1 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so;  LAPACK version 3.12.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Etc/UTC
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] seqmagick_0.1.9      Biostrings_2.81.9    Seqinfo_1.3.2       
##  [4] XVector_0.53.0       IRanges_2.47.5       S4Vectors_0.51.10   
##  [7] BiocGenerics_0.59.12 generics_0.1.4       yulab.utils_0.2.5   
## [10] prettydoc_0.4.1     
## 
## loaded via a namespace (and not attached):
##  [1] crayon_1.5.3     cli_3.6.6        knitr_1.52       rlang_1.3.0     
##  [5] xfun_0.61        otel_0.2.0       jsonlite_2.0.0   buildtools_1.0.0
##  [9] htmltools_0.5.9  maketools_1.3.2  sys_3.4.3        sass_0.4.10     
## [13] rmarkdown_2.32   rappdirs_0.3.4   evaluate_1.0.5   jquerylib_0.1.4 
## [17] fastmap_1.2.0    yaml_2.3.12      lifecycle_1.0.5  compiler_4.6.1  
## [21] fs_2.1.0         digest_0.6.39    R6_2.6.1         bslib_0.12.0    
## [25] tools_4.6.1      cachem_1.1.0

seqmagick: sequence manipulation