R/readSFF.R

In msettles/rSFFreader: rSFFreader reads in sff files generated by Roche 454 and Life Sciences Ion Torrent sequencers

### Sample Data
load454SampleData <- function() {
  readSff(system.file("extdata","Small454Test.sff",package="rSFFreader"))
}
loadIonSampleData <- function() {
  readSff(system.file("extdata","SmallTorrentTest.sff",package="rSFFreader"))
}

#lkup_seq <- get_seqtype_conversion_lookup("B", "DNA")
#lkup_seq <- get_seqtype_conversion_lookup("B", "DNA")
#ans <- .Call("read_sff","inst/extdata/SmallTest.sff",TRUE, TRUE, lkup_seq, NULL,TRUE,"rSFFreader")
#new("SffReadsQ",sread=ans[[2]],quality=FastqQuality(ans[[3]]),
#  qualityClip=ans[[4]],adapterClip=ans[[5]],clipMode="Full",header=ans[[1]])

### Some clip points (specifically adapter clip points) will have 0 values, reset to full length)
.fixSFFclipPoints2Iranges <- function(widths,IR){
  clipL <- pmax(1, start(IR)) 
  clipR <- ifelse(end(IR) == 0 , widths,end(IR))
  clipR <- pmin(widths,clipR)
  return(solveUserSEW(widths,clipL,clipR))  
}

## functions 
## returns the contents of the SFF file into either a SffReads or SffReadsQ class, which acts and behaves similar to
## the ShortRead and ShortReadQ classes from package ShortRead
readSff <- function(filenames, use.qualities=TRUE, use.names=TRUE, use.flowgrams=FALSE, clipMode=c("full","adapter","quality","raw"), verbose=TRUE){
  if (!use.names) warning ("Currently use.names is not used, by default names will always be returned.")
  stopifnot(file.exists(filenames))
  clipMode <- match.arg(clipMode)
	if (!isTRUEorFALSE(use.names))
		stop("'use.names' must be TRUE or FALSE")
  if (!isTRUEorFALSE(use.flowgrams))
    stop("'use.flowgrams' must be TRUE or FALSE")
	if (!isTRUEorFALSE(verbose))
		stop("'verbose' must be TRUE or FALSE")
	lkup_seq <- get_seqtype_conversion_lookup("B", "DNA")
	ans <- .Call("read_sff",filenames,use.names, use.flowgrams, lkup_seq, NULL,verbose)
  widths <- width(ans[["sread"]])
  if (use.flowgrams){
    SffReadsF(sread=ans[["sread"]],quality=ans[["quality"]],
              qualityIR=.fixSFFclipPoints2Iranges(widths,ans[["qualityClip"]]),
              adapterIR=.fixSFFclipPoints2Iranges(widths,ans[["adapterClip"]]),
              clipMode=clipMode,header=ans[["header"]],
              flowgram=matrix(ans[["flowgrams"]], byrow=TRUE, nrow=length(ans[["sread"]])),
              flowindices=split(ans[["flowIndices"]], rep(seq.int(1,length(widths)), times=widths)))
  } else if (use.qualities){
    	SffReadsQ(sread=ans[["sread"]],quality=ans[["quality"]],
                qualityIR=.fixSFFclipPoints2Iranges(widths,ans[["qualityClip"]]),
                adapterIR=.fixSFFclipPoints2Iranges(widths,ans[["adapterClip"]]),
                clipMode=clipMode,header=ans[["header"]])
    } else {
    	SffReads(ans[["sread"]],
               qualityIR=.fixSFFclipPoints2Iranges(widths,ans[["qualityClip"]]),
               adapterIR=.fixSFFclipPoints2Iranges(widths,ans[["adapterClip"]]),
               clipMode=clipMode,header=ans[["header"]])
    }
}

## Returns a list of size 2
readSffGeometry <- function(filenames) {
	stopifnot(file.exists(filenames))
 	sffgeometry <- .Call("sff_geometry", filenames)
    names(sffgeometry) <- c("nReads","Read_Widths")
	return(sffgeometry)
}


readSffHeader <- function(filenames,verbose=TRUE) {
    stopifnot(file.exists(filenames))
    if (!isTRUEorFALSE(verbose))
		stop("'verbose' must be TRUE or FALSE")
	ans <- .Call("read_sff_header", filenames,verbose)
	new("SffHeader", header=ans)
}