---
title: "seqmagick: sequence manipulation"
author: 
- name: Guangchuang Yu
  email: guangchuangyu@gmail.com
  affiliation: Department of Bioinformatics, School of Basic Medical Sciences, Southern Medical University
date: "`r Sys.Date()`"
bibliography: seqmagick.bib
biblio-style: apalike
output:
  prettydoc::html_pretty:
    toc: true
    highlight: github
    theme: cayman
  pdf_document:
    toc: true
vignette: >
  %\VignetteIndexEntry{seqmagick introduction}
  %\VignetteEngine{knitr::rmarkdown}
  %\usepackage[utf8]{inputenc}
---

```{r style, echo=FALSE, results="asis", message=FALSE}
knitr::opts_chunk$set(tidy = FALSE,
		   message = FALSE)
```


```{r echo=FALSE, results="hide", message=FALSE}
library(yulab.utils)
library(Biostrings)
library("seqmagick")
```

# Supported file formats

`seqmagick` works with `r Biocpkg("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

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

Compressed files can be written and read directly:

```{r}
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]))
```

## PHYLIP

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

## CLUSTAL and STOCKHOLM

```{r eval=FALSE}
x <- clw_read("alignment.clw")
x <- sth_read("alignment.sth")
```

## MEGA

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

## GenBank

```{r eval=FALSE}
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 `r Biocpkg("Biostrings")` objects directly) fetch records by accession number.

```{r 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'))
```

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

```{r 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')
```

Converting back from PHYLIP to FASTA is symmetric:

```{r eval=FALSE}
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()`:

```{r 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")
```

# Sequence manipulation

## Filtering, alignment and consensus

```{r 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)
```

## Renaming sequences

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

```{r}
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]
```

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

```{r}
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)
```

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

```{r eval=FALSE}
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 `r Biocpkg("Biostrings")` object): number of sequences, length range, GC content and proportion of ambiguous characters.

```{r}
fa_summary(fa_file)
```

Works on in-memory objects as well:

```{r}
fa_summary(x[1:10])
```

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

```{r}
head(seqlen(fa_file))
```

# Bugs/Feature requests

If you have any, [let me know](https://github.com/YuLab-SMU/seqmagick/issues). Thx!

# Session info

Here is the output of `sessionInfo()` on the system on which this document was compiled:
```{r echo=FALSE}
sessionInfo()
```


# References