I’ve been working on a lot of AFLP data this winter. I’d really like to be able to do all the analysis in R, for a few reasons. First, it would mean no more fighting with GeneMapper, which is incredibly frustrating: it’s Windows-only, expensive, closed-source and painfully underpowered for the job. Second, presumably if I can figure out how to code this myself I will develop a deeper understanding of the system. And third, if I can get the code working in R, I will be able to automate most of the process.

There are two R projects already in progress for working with AFLP data. RawGeno is one option. It doesn’t yet allow for importing fsa files directly, but the example scripts provide some clues about how to do this. I couldn’t get the code to work as written, but I was able to steal some ideas from it.

The other R package is AFLP. This package includes a read.fsa() function, but it doesn’t seem to work yet. I understand they’ve only recently switched to ABI sequencers, and haven’t yet updated their code. AFLP also combines reading the fsa files, calibrating the sizing, and defining the bins into one step. That’s a sensible thing to do, but I’m not that clever. I need to break things into small pieces if I hope to get anywhere.

Since one of my goals is self-education, I’m not concerned about duplicating some of the effort of these other projects. In fact, I’m going to try and steal as much as I can from them. That’s one of the benefits of Free Software, we get to learn from each other.

Step one, reading the raw data

Lucky for me, most of the work involved in actually getting the contents of an .fsa file into R has already been done, via the package seqinr. All that I need to do is extract the useful bits and reformat it into a data.frame.

library(seqinr)

read.fsa <- function(files = NULL, path = "./", sig.channel = 1:3, lad.channel = 105, pretrim = FALSE,
                     posttrim = ".fsa", thresh = -100, verbose = TRUE){

  if(is.null(files))
    files <- list.files(path, pattern = "\\.fsa$", full.names = TRUE)
  else
    files <- paste(path, files, sep = "")

  res <- do.call(rbind, lapply(files, function(file) {
    if (verbose) message(file)
    abif <- read.abif(file)
    tag <- tag.trimmer(basename(file), pretrim, posttrim)
    
    lad.dat <- abif$Data[[paste('DATA.', lad.channel, sep='')]]
    
    res1 <- data.frame(tag = as.character(rep(tag, length(lad.dat))),
                       chan = as.character(rep("standard", length(lad.dat))),
                       time = as.numeric(1:length(lad.dat)),
                       peak = as.numeric(lad.dat))
    
    for (i in sig.channel) {
      chan.dat <- abif$Data[[paste('DATA.', i, sep='')]]
      res1 <- rbind(res1, data.frame(tag = as.character(rep(tag, length(chan.dat))),
                                     chan = as.character(rep(i, length(chan.dat))),
                                     time = as.numeric(1:length(chan.dat)),
                                     peak = as.numeric(chan.dat)))
    }
    res1
  }))
    
  if (thresh > -10) res <- subset(res, peak > thresh)
  return(res)
}

tag.trimmer <- function(x, pretrim = FALSE, posttrim = FALSE) {
  if(! is.na(pretrim)) {
    x <- sub(paste("^", pretrim, sep = ""), "", x)
  }
  if(! is.na(posttrim)){
    x <- sub(paste(posttrim, "$", sep = ""), "", x)
  }
  x
}

sig.channel is a vector of the DATA channels to read from the fsa file. I’m using FAM dye, which gets recorded in channel 1. lad.channel is the DATA channel where the size standard is found. We use the orange dye for the ladder, which is in channel 105. pretrim and posttrim are conveniences, for removing leading and trailing strings from the filenames, via tag.trimmer.

fsa <- read.fsa(path = "./path/to/fsa/files/", sig.channel = 1,
                pretrim = "AFLP.*AFLP_", posttrim = "-5_Frag.*fsa")   

head(fsa)
      tag     chan time peak
1 QCWR-25 standard    1   -3
2 QCWR-25 standard    2    3
3 QCWR-25 standard    3    1
4 QCWR-25 standard    4   -3
5 QCWR-25 standard    5   -2
6 QCWR-25 standard    6   -1

The actual data, in my case, is composed of 8959 rows for each sample x dye combination. Each row is the reading from the laser at that point in the run (the time). In other words, the size of the fragments that are migrating past the window at that particular time. Since we have multiple readings for each time, the time column allows us to refer to information from different dyes and different samples that were detected at the same time. peak is the strength of the fluorescence associated with each sample x dye x time combination. The negative numbers are obviously noise. You can use the thresh argument to clear out all rows that are below a particular fluorescence value. This isn’t necessary unless you’ve got a really big data set. I ran this on nearly 200 samples and had no problems – I don’t think you’re likely to run into issues with less than 1000 samples.

This isn’t useful yet. First we need to convert the times into actual fragment sizes, in base pairs. In the meantime, we can at least plot our raw data in R now:

plot(peak ~ time, data = subset(fsa, tag == "QCWR-25" & chan == "standard"), 
       col = "orange", type = 'l', ylim = c(0, 4000), xlim = c(700, 4000))
points(peak ~ time, data = subset(fsa, tag == "QCWR-25" & chan == "1"), 
           col = "blue", type = 'l')

Raw fsa data produced by read.fsa

Next up is finding the peaks in each channel, matching up the size standard peaks to the known values, and using that to convert the rest of the peaks from time to base-pairs.