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 seq                                             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 seq                                              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 24.04.4 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.12.0 
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0  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=C              
##  [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: Asia/Shanghai
## 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.80.1   Seqinfo_1.2.0      
## [4] XVector_0.52.0      IRanges_2.46.0      S4Vectors_0.50.2   
## [7] BiocGenerics_0.58.1 generics_0.1.4      yulab.utils_0.2.5  
## 
## loaded via a namespace (and not attached):
##  [1] crayon_1.5.3    cli_3.6.6       knitr_1.51      rlang_1.3.0    
##  [5] xfun_0.60       otel_0.2.0      jsonlite_2.0.0  htmltools_0.5.9
##  [9] sass_0.4.10     rmarkdown_2.31  rappdirs_0.3.4  evaluate_1.0.5 
## [13] jquerylib_0.1.4 prettydoc_0.4.1 fastmap_1.2.0   yaml_2.3.12    
## [17] lifecycle_1.0.5 compiler_4.6.1  fs_2.1.0        digest_0.6.39  
## [21] R6_2.6.1        bslib_0.12.0    tools_4.6.1     cachem_1.1.0

References