micropan (version 1.2)

hmmerScan: Scanning a profile Hidden Markov Model database

Description

Scanning FASTA formatted protein files against a database of pHMMs using the HMMER3 software.

Usage

hmmerScan(in.files, db, out.folder, threads = 0, verbose = TRUE)

Arguments

in.files

A character vector of file names.

db

The full name of the database to scan.

out.folder

The name of the folder to put the result files.

threads

Number of CPU's to use.

verbose

Logical indicating if textual output should be given to monitor the progress.

Value

This function produces files in the folder specified by out.folder. Existing files are never overwritten by hmmerScan, if you want to re-compute something, delete the corresponding result files first.

Details

The HMMER3 software is purpose-made for handling profile Hidden Markov Models (pHMM) describing patterns in biological sequences (Eddy, 2008). This function will make calls to the HMMER3 software to scan FASTA files of proteins against a pHMM database.

The files named in in.files must contain FASTA formatted protein sequences. These files should be prepared by panPrep to make certain each sequence, as well as the file name, has a GID-tag identifying their genome. The database named in db must be a HMMER3 formatted database. It is typically the Pfam-A database, but you can also make your own HMMER3 databases, see the HMMER3 documentation for help.

hmmerScan will query every input file against the named database. The database contains profile Hidden Markov Models describing position specific sequence patterns. Each sequence in every input file is scanned to see if some of the patterns can be matched to some degree. Each input file results in an output file with the same GID-tag in the name. The result files give tabular output, and are plain text files. See readHmmer for how to read the results into R.

Scanning large databases like Pfam-A takes time, usually several minutes per genome. The scan is set up to use only 1 cpu per scan by default. By increasing threads you can utilize multiple CPUs, typically on a computing cluster. Our experience is that from a multi-core laptop it is better to start this function in default mode from mutliple R-sessions. This function will not overwrite an existing result file, and multiple parallel sessions can write results to the same folder.

References

Eddy, S.R. (2008). A Probabilistic Model of Local Sequence Alignment That Simplifies Statistical Significance Estimation. PLoS Computational Biology, 4(5).

See Also

panPrep, readHmmer.

Examples

# NOT RUN {
# This example requires the external HMMER software
# Using two files in the micropan package
xpth <- file.path(path.package("micropan"),"extdata")
prot.file <- file.path(xpth,"Example_proteins_GID1.fasta.xz")
db <- "microfam.hmm"
db.files <- file.path(xpth,paste(db,c(".h3f.xz",".h3i.xz",".h3m.xz",".h3p.xz"),sep=""))

# We need to uncompress them first...
prot.tf <- tempfile(pattern="GID1.fasta",fileext=".xz")
s <- file.copy(prot.file,prot.tf)
prot.tf <- xzuncompress(prot.tf)
db.tf <- paste(tempfile(),c(".h3f.xz",".h3i.xz",".h3m.xz",".h3p.xz"),sep="")
s <- file.copy(db.files,db.tf)
db.tf <- unlist(lapply(db.tf,xzuncompress))
db.name <- gsub("\\",.Platform$file.sep,sub(".h3f$","",db.tf[1]),fixed=T)

# Scanning the FASTA-file against microfam0...
tmp.dir <- tempdir()
hmmerScan(in.files=prot.tf,db=db.name,out.folder=tmp.dir)

# Reading results
db.nm <- rev(unlist(strsplit(db.name,split=.Platform$file.sep)))[1]
hmm.file <- file.path(tmp.dir,paste("GID1_vs_",db.nm,".txt",sep=""))
hmm.tab <- readHmmer(hmm.file)

# ...and cleaning...
s <- file.remove(prot.tf)
s <- file.remove(sub(".xz","",db.tf))
s <- file.remove(hmm.file)
# }
# NOT RUN {
# }