## ----style, echo=FALSE, results="asis", message=FALSE-------------------------
knitr::opts_chunk$set(tidy = FALSE,
		   message = FALSE)

## ----echo=FALSE, results="hide", message=FALSE--------------------------------
library(yulab.utils)
library(Biostrings)
library("seqmagick")

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

## -----------------------------------------------------------------------------
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]))

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

## ----eval=FALSE---------------------------------------------------------------
# x <- clw_read("alignment.clw")
# x <- sth_read("alignment.sth")

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

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

## ----eval=FALSE---------------------------------------------------------------
# ## 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'))

## ----eval=FALSE---------------------------------------------------------------
# ## 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')

## ----eval=FALSE---------------------------------------------------------------
# phy2fas(tmpphy, alnfas, type = 'interleaved')

## ----eval=FALSE---------------------------------------------------------------
# ## 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")

## ----eval=FALSE---------------------------------------------------------------
# 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)

## -----------------------------------------------------------------------------
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]

## -----------------------------------------------------------------------------
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)

## ----eval=FALSE---------------------------------------------------------------
# fa_rename("seqs.fas", "map.txt", outfile="renamed.fas")

## -----------------------------------------------------------------------------
fa_summary(fa_file)

## -----------------------------------------------------------------------------
fa_summary(x[1:10])

## -----------------------------------------------------------------------------
head(seqlen(fa_file))

## ----echo=FALSE---------------------------------------------------------------
sessionInfo()