Rbowtie (version 1.12.0)

SpliceMap: Interface to SpliceMap

Description

The following function can be used to call the SpliceMap binaries.

We recommend to use the QuasR package instead of using SpliceMap directly. The QuasR package provides a simpler interface then Rbowtie and covers the whole analysis workflow of typical ultra-high throughput sequencing experiments, starting from the raw sequence reads, over pre-processing and alignment, up to quantification.

Usage

SpliceMap(cfg)

Arguments

cfg

A list containing named elements with SpliceMap parameters (see sQuoteDetails or SpliceMap documentation).

Value

An invisible character vector of length one with the file name of the generated output SAM file.An exception is thrown if a failure is detected at one of the many steps.

Details

The SpliceMap function performes the same steps as the runSpliceMap binary from the SpliceMap software package, but using R functions to improve compatibility on Windows.

While the original SpliceMap software package is able to call different tools to find sub-read alignments, the SpliceMap function from the Rbowtie package works only with bowtie contained in the package itself. Further modifications from the original version include the reporting of unmapped reads at the end of the output sam file, and the restriction to a single (pair) of input sequence file(s). The cfg argument is a list with SpliceMap configuration parameters that would normally be specified using the SpliceMap configuration file. Here is a list of supported parameters extracted from the sample config file that is distributed with SpliceMap:

• genome_dir (single character value) Directory of the chromosome files in FASTA format, or path to a single FASTA file containing all chromosomes. If a directory, each chromosome can be in a separate file or chromosomes can be concatenated, ie. chr1.fa, chr2.fa, ...

reads_list1 and reads_list2 (single character values) These are the two input sequence files. reads_list2 can be missing if reads are not paired-end. Note: reads_list1 must be the first pair, and pair-reads should be in the “forward-reverse” format.

read_format (single character value) Format of the sequencer reads, also make sure reads are not split over multiple lines. Choices are: FASTA, FASTQ, RAW

quality_format (single character value) Format of the quality string if FASTQ is used. Choices are:

• phred-33 Phred base 33 (same as Sanger format)
• phred-64 Phred base 64 (same as Illumina 1.3+)
• solexa Format used by solexa machines

outfile (single character value) Name of the output file (spliced alignments in SAM format). If a file with this name already exists, SpliceMap will stop with an exception. Note that unmapped reads will be appended to the end of the output file (for paired-end experiments, only pairs without alignments for any read are considered unmapped).

temp_path (single character value) Directory name of the directory that stores temporary files. All temporary files will be created in a subfolder of temp_path and will be removed when SpliceMap finishes successfully or failed at an intermediate step.

max_intron (single integer value) Maximum intron size, this is absolute 99th-percentile maximum. Introns beyond this size will be ignored. If you don't set this, we will assume a mamalian genome (400,000)

min_intron (single integer value) 25-th intron size, this is the lower 25th-percentile intron size. This is not the smallest size that SpliceMap will search. That is about ~25bp. If you don't set this, we will assume a mammalian genome (20,000)

max_multi_hit (single integer value) Maximum number of multi-hits. If a 25-mer seed has more than this many multi-hits, it will be discarded. Default is 10.

full_read_length (single integer value) Full read length. SpliceMap will only use the first “full_read_length” bp for mapping. If the read is shorter than “full_read_length”, the full read will be used before head clip. If you don't set this parameter, SpliceMap will use as many as possible. This is for the case where the reads might have N's at the end. It is always desireable to cut off the N's

head_clip_length (single integer value) Number of bases to clip off the head of the read. This clipping is applied after “full_read_length”

seed_mismatch (single integer value) Number of mismatches allowed in half-seeding. Choices are 0,1(default) or 2

read_mismatch (single integer value) Maximum number of mismatches allowed in entire read. No limit on value, however SpliceMap can only identify reads with a maximum of 2 mismatches per 25bp. Default is 2.

max_clip_allowed (single integer value) Maximum number of bases allowed to be soft clipped from the ends of reads during alignment. This is required as mismatches near junctions could cause parts of a a read to not map. Default is 40.

num_chromosome_together (single integer value) Number of chromosomes to process at once, to take advantage of multi-core systems. The child processes are created using the makeCluster function from the parallel package. This is not threading, so it will take extra memory. However, running 2 at a time should be fine on current hardware. Default = 2

bowtie_base_dir (single character value) Base of bowtie index, this should be the same genome as the chromosome files eg. if you bowtie files are “genome/hg18/genome.1.ewbt”, ... then your base dir is “genome/hg18/genome”

num_threads (single integer values) Number of threads to use for bowtie mapping. Default value is 2

try_hard (single character value) Try hard? Choices are “yes” or “no”. Default value is “yes” (about 15% slower).

selectSingleHit (single logical value) If TRUE and multiple alignments are found for a read (pair), only a single alignment is selected randomly and reported. This is a new paramater only available in the R version of SpliceMap, but not in the original implementation that could report several alignments per read (pair).

The following parameters are mandatory:

• genome_dir
• reads_list1 (reads_list2 for paired-end experiments)
• read_format
• quality_format
• bowtie_base_dir
• outfile
• temp_path
• num_threads
• quality_format
• selectSingleHit

References

Au KF, Jiang H, Lin L, Xing Y, Wong WH. Detection of splice junctions from paired-end RNA-seq data by SpliceMap. Nucleic Acids Research, 38(14):4570-8 (2010).

See Also

makeCluster from package parallel

Examples

## Building a bowtie index
refDir <- system.file(package="Rbowtie", "samples", "refs")
indexDir <- file.path(tempdir(), "refsIndex")

tmp <- bowtie_build(references=dir(refDir, full=TRUE), outdir=indexDir, prefix="index", force=TRUE)

## Alignments
readsFiles <- system.file(package="Rbowtie", "samples", "reads", "reads.fastq")
samFiles <- file.path(tempdir(), "splicedAlignments.sam")

cfg <- list(genome_dir=refDir,
            reads_list1=readsFiles,
            read_format="FASTQ",
            quality_format="phred-33",
            outfile=samFiles,
            temp_path=tempdir(),
            max_intron=400000,
            min_intron=20000,
            max_multi_hit=10,
            seed_mismatch=1,
            read_mismatch=2,
            num_chromosome_together=2,
            bowtie_base_dir=file.path(indexDir, "index"),
            num_threads=4,
            try_hard="yes",
            selectSingleHit=TRUE)
            
res <- SpliceMap(cfg)