R: Read/write an XStringSet object from/to a file

XStringSet-io {Biostrings}

R Documentation

Read/write an XStringSet object from/to a file

Description

Functions to read/write an XStringSet object from/to a file.

Usage

## Read FASTA (or FASTQ) files in an XStringSet object:
readBStringSet(filepath, format="fasta",
               nrec=-1L, skip=0L, seek.first.rec=FALSE, use.names=TRUE)
readDNAStringSet(filepath, format="fasta",
               nrec=-1L, skip=0L, seek.first.rec=FALSE, use.names=TRUE)
readRNAStringSet(filepath, format="fasta",
               nrec=-1L, skip=0L, seek.first.rec=FALSE, use.names=TRUE)
readAAStringSet(filepath, format="fasta",
               nrec=-1L, skip=0L, seek.first.rec=FALSE, use.names=TRUE)

## Extract basic information about FASTA (or FASTQ) files
## without actually loading the sequence data:
fasta.seqlengths(filepath,
                 nrec=-1L, skip=0L, seek.first.rec=FALSE, seqtype="B",
                 use.names=TRUE)
fasta.index(filepath,
            nrec=-1L, skip=0L, seek.first.rec=FALSE, seqtype="B")

fastq.geometry(filepath, nrec=-1L, skip=0L, seek.first.rec=FALSE)

## Write an XStringSet object to a FASTA (or FASTQ) file:
writeXStringSet(x, filepath, append=FALSE,
                compress=FALSE, compression_level=NA, format="fasta", ...)

## Serialize an XStringSet object:
saveXStringSet(x, objname, dirpath=".", save.dups=FALSE, verbose=TRUE)

Arguments

`filepath`	A character vector (of arbitrary length when reading, of length 1 when writing) containing the path(s) to the file(s) to read or write. Reading files in gzip format (which usually have the '.gz' extension) is supported. Note that special values like `""` or `"\|cmd"` (typically supported by other I/O functions in R) are not supported here. Also `filepath` cannot be a connection.
`format`	Either `"fasta"` (the default) or `"fastq"`.
`nrec`	Single integer. The maximum of number of records to read in. Negative values are ignored.
`skip`	Single non-negative integer. The number of records of the data file(s) to skip before beginning to read in records.
`seek.first.rec`	`TRUE` or `FALSE` (the default). If `TRUE`, then the reading function starts by setting the file position indicator at the beginning of the first line in the file that looks like the beginning of a FASTA (if `format` is `"fasta"`) or FASTQ (if `format` is `"fastq"`) record. More precisely this is the first line in the file that starts with a '>' (for FASTA) or a '@' (for FASTQ). An error is raised if no such line is found. Normal parsing then starts from there, and everything happens like if the file actually started there. In particular it will be an error if this first record is not a valid FASTA or FASTQ record. Using `seek.first.rec=TRUE` is useful for example to parse GFF3 files with embedded FASTA data.
`use.names`	Should the returned vector be named? For FASTA the names are taken from the record description lines. For FASTQ they are taken from the record sequence ids. Dropping the names can help reducing memory footprint e.g. for a FASTQ file containing millions of reads.
`seqtype`	A single string specifying the type of sequences contained in the FASTA file(s). Supported sequence types: `"B"` for anything i.e. any letter is a valid one-letter sequence code. `"DNA"` for DNA sequences i.e. only letters in `DNA_ALPHABET` (case ignored) are valid one-letter sequence codes. `"RNA"` for RNA sequences i.e. only letters in `RNA_ALPHABET` (case ignored) are valid one-letter sequence codes. `"AA"` for Amino Acid sequences. Currently treated as `"B"` but this will change in the near future i.e. only letters in `AA_ALPHABET` (case ignored) will be valid one-letter sequence codes. Invalid one-letter sequence codes are ignored with a warning.
`x`	For `writeXStringSet`, the object to write to `file`. For `saveXStringSet`, the object to serialize.
`append`	`TRUE` or `FALSE`. If `TRUE` output will be appended to `file`; otherwise, it will overwrite the contents of `file`. See `?cat` for the details.
`compress`	Like for the `save` function in base R, must be `TRUE` or `FALSE` (the default), or a single string specifying whether writing to the file is to use compression. The only type of compression supported at the moment is `"gzip"`. Passing `TRUE` is equivalent to passing `"gzip"`.
`compression_level`	Not implemented yet.
`...`	Further format-specific arguments. If `format="fasta"`, the `width` argument (single integer) can be used to specify the maximum number of letters per line of sequence. If `format="fastq"`, the `qualities` argument (BStringSet object) can be used to specify the qualities. If the qualities are omitted, then the fake quality ';' is assigned to each letter in `x` and written to the file.
`objname`	The name of the serialized object.
`dirpath`	The path to the directory where to save the serialized object.
`save.dups`	`TRUE` or `FALSE`. If `TRUE` then the `Dups` object describing how duplicated elements in `x` are related to each other is saved too. For advanced users only.
`verbose`	`TRUE` or `FALSE`.

Details

gzip compression is supported by reading and writing functions on all platforms.

readDNAStringSet and family (i.e. readBStringSet, readDNAStringSet, readRNAStringSet and readAAStringSet) load sequences from an input file (or multiple input files) into an XStringSet object. When multiple input files are specified, all must have the same format (i.e. FASTA or FASTQ) and files with different compression types can be mixed with non-compressed files. The files are read in the order they were specified and the sequences are stored in the returned object in the order they were read.

Only FASTA and FASTQ files are supported for now. The read qualities stored in FASTQ files are ignored by readDNAStringSet and family. When multiple input FASTQ files are specified, all must have the same "width" (i.e. all their sequences must have the same length).

The fasta.seqlengths utility returns an integer vector with one element per FASTA record in the input files. Each element is the length of the sequence found in the corresponding record, that is, the number of valid one-letter sequence codes in the record. See description of the seqtype argument above for how to control the set of valid one-letter sequence codes.

The fasta.index utility returns a data frame with 1 row per FASTA record in the input files and the following columns:

recno: The rank of the record in the (virtually) concatenated input files.
fileno: The rank of the file where the record is located.
offset: The offset of the record relative to the start of the file where it's located. Measured in bytes.
desc: The description line (a.k.a. header) of the record.
seqlength: The length of the sequence in the record (not counting invalid letters).
filepath: The path to the file where the record is located. Always a local file, so if the user specified a remote file, this column will contain the path to the downloaded file.

A subset of this data frame can be passed to readDNAStringSet and family for direct access to an arbitrary subset of sequences. More precisely, if fai is a FASTA index that was obtained with fasta.index(filepath, ..., seqtype="DNA"), then readDNAStringSet(fai[i, ]) is equivalent to readDNAStringSet(filepath, ...)[i] for any valid subscript i, except that the former only loads the requested sequences in memory and thus will be more memory efficient if only a small subset of sequences is requested.

The fastq.geometry utility returns an integer vector describing the "geometry" of the FASTQ files i.e. a vector of length 2 where the first element is the total number of FASTQ records in the files and the second element the common "width" of these files (this width is NA if the files contain no FASTQ records or records with different widths).

writeXStringSet writes an XStringSet object to a file. Like with readDNAStringSet and family, only FASTA and FASTQ files are supported for now. WARNING: Please be aware that using writeXStringSet on a BStringSet object that contains the '\n' (LF) or '\r' (CR) characters or the FASTA markup characters '>' or ';' is almost guaranteed to produce a broken FASTA file!

Serializing an XStringSet object with saveXStringSet is equivalent to using the standard save mechanism. But it will try to reduce the size of x in memory first before calling save. Most of the times this leads to a much reduced size on disk.

References

http://en.wikipedia.org/wiki/FASTA_format

Examples

## ---------------------------------------------------------------------
## A. READ/WRITE FASTA FILES
## ---------------------------------------------------------------------

## Read a non-compressed FASTA files:
filepath1 <- system.file("extdata", "someORF.fa", package="Biostrings")
fasta.seqlengths(filepath1, seqtype="DNA")
x1 <- readDNAStringSet(filepath1)
x1

## Read a gzip-compressed FASTA file:
filepath2 <- system.file("extdata", "someORF.fa.gz", package="Biostrings")
fasta.seqlengths(filepath2, seqtype="DNA")
x2 <- readDNAStringSet(filepath2)
x2

## Sanity check:
stopifnot(identical(as.character(x1), as.character(x2)))

## Read 2 FASTA files at once:
filepath3 <- system.file("extdata", "fastaEx.fa", package="Biostrings")
fasta.seqlengths(c(filepath2, filepath3), seqtype="DNA")
x23 <- readDNAStringSet(c(filepath2, filepath3))
x23

## Sanity check:
x3 <- readDNAStringSet(filepath3)
stopifnot(identical(as.character(x23), as.character(c(x2, x3))))

## Use a FASTA index to load only an arbitrary subset of sequences:
filepath4 <- system.file("extdata", "dm3_upstream2000.fa.gz",
                         package="Biostrings")
fai <- fasta.index(filepath4, seqtype="DNA")
head(fai)
head(fai$desc)
i <- sample(nrow(fai), 10)  # randomly pick up 10 sequences
x4 <- readDNAStringSet(fai[i, ])

## Sanity check:
stopifnot(identical(as.character(readDNAStringSet(filepath4)[i]),
                    as.character(x4)))

## Write FASTA files:
out23a <- tempfile()
writeXStringSet(x23, out23a)
out23b <- tempfile()
writeXStringSet(x23, out23b, compress=TRUE)
file.info(c(out23a, out23b))$size

## Sanity checks:
stopifnot(identical(as.character(readDNAStringSet(out23a)),
                    as.character(x23)))
stopifnot(identical(readLines(out23a), readLines(out23b)))

## ---------------------------------------------------------------------
## B. READ/WRITE FASTQ FILES
## ---------------------------------------------------------------------

filepath <- system.file("extdata", "s_1_sequence.txt",
                        package="Biostrings")
fastq.geometry(filepath)
readDNAStringSet(filepath, format="fastq")

library(BSgenome.Celegans.UCSC.ce2)
## Create a "sliding window" on chr I:
sw_start <- seq.int(1, length(Celegans$chrI)-50, by=50)
sw <- Views(Celegans$chrI, start=sw_start, width=10)
my_fake_shortreads <- as(sw, "XStringSet")
my_fake_ids <- sprintf("ID%06d",  seq_len(length(my_fake_shortreads)))
names(my_fake_shortreads) <- my_fake_ids
my_fake_shortreads

## Fake quality ';' will be assigned to each base in 'my_fake_shortreads':
out2 <- tempfile()
writeXStringSet(my_fake_shortreads, out2, format="fastq")

## Passing qualities thru the 'qualities' argument:
my_fake_quals <- rep.int(BStringSet("DCBA@?>=<;"),
                         length(my_fake_shortreads))
my_fake_quals
out3 <- tempfile()
writeXStringSet(my_fake_shortreads, out3, format="fastq",
                qualities=my_fake_quals)

## ---------------------------------------------------------------------
## C. SERIALIZATION
## ---------------------------------------------------------------------
saveXStringSet(my_fake_shortreads, "my_fake_shortreads", dirpath=tempdir())

[Package Biostrings version 2.46.0 Index]