readmmseq <- function( mmseq_files = grep("gene|identical", dir( pattern="\\.mmseq$" ), value=TRUE, invert=TRUE),
		 sample_names = sub("\\.mmseq","", mmseq_files), normalize=TRUE) {
	
	cat( "Using mmseq files:\n", paste("\t", mmseq_files, sep="", collapse="\n"), "\n") 
	cat( "Using sample names:\n", paste("\t", sample_names, sep="", collapse="\n"), "\n") 
	if( length(mmseq_files) == 0 || length(sample_names)  != length(mmseq_files) ) {
		cat( "ERROR:\n\tplease define mmseq files and sample names for each file.\n" )
		return(NULL)
	}
	
	# read in mmseq data
	mu = c()
	mcse = c()
	effective_lengths = c()
	true_lengths = c()
  n_reads <- c()
  unique_hits <- c()
	cat( "Reading mmseq output files\n" )
	for( i in 1:length(mmseq_files) ) {
		cat( "\t", mmseq_files[i], "\t" )
		dat = read.table( mmseq_files[i], header=TRUE, stringsAsFactors=FALSE, sep="\t" )
    if(! any(colnames(dat)=="log_mu_gibbs") ||
       ! any(colnames(dat)=="mcse") || 
       ! any(colnames(dat)=="unique_hits") ||
       ! any(colnames(dat)=="effective_length") ||
       ! any(colnames(dat)=="true_length")) {
      stop("Colnames must contain log_mu_gibbs, mcse, effective_length, true_length and unique_hits.")
    }
    mu <- cbind(mu, dat$log_mu_gibbs)
    mcse <- cbind(mcse, dat$mcse)
    unique_hits <- cbind(unique_hits, dat$unique_hits)
    effective_lengths <- cbind(effective_lengths, dat$effective_length)
    true_lengths <- cbind(true_lengths, dat$true_length)
    if(i==1) {
      rownames(mu) = rownames(mcse)= rownames(unique_hits)=rownames(true_lengths)=rownames(effective_lengths)=dat[,grepl("id", colnames(dat))]
    } else {
      if(! all(rownames(mu) == dat[,1])) {
        stop("Attempting to join incompatible mmseq files (different features).\n")
      }
    }

    ta <- readLines(mmseq_files[i], n=100) # header should be no more than 100 lines
    ta <- strsplit(ta[grep("# Mapped fragments", ta)], " ")[[1]]
    n_reads <- c(n_reads, as.numeric(ta[length(ta)]))
		cat( "found", n_reads[length(n_reads)], "aligned reads (or read pairs)\n" )
	}
  colnames(mu)=colnames(mcse)=colnames(unique_hits)=colnames(effective_lengths)=colnames(true_lengths)=sample_names
  names(n_reads) <- sample_names

  if(!all(is.na(true_lengths)) && any(apply(true_lengths,1, function(x) { !all(x==x[1])}))) {
    stop("True transcript lengths differ across files.\n")
  }

  # un-log the estimates 
  cat( "Unlogging and unstandardising the estimates\n" )
  counts = exp(1)^(mu)
  
  # multiply by library size
  counts= sapply( colnames(counts), function(x) { counts[,x]*n_reads[x] })
  rownames(counts)=rownames(mu)

  # multiply by transcript length and divide by 10^9
  counts = sapply( colnames(mu), function(x) { counts[,x] = (counts[,x]*effective_lengths[,x])/(1e3*1e6) } )
  counts[is.na(counts)] = 0 # the ones that had no expression in MMSEQ
	
  if(normalize) {
    cat("Normalizing expression estimates with log factors:\n")
    sf <- apply(mu - rowMeans(mu) ,2, function(x){median(x[is.finite(x)])})
    for(i in 1:length(sf)) {
      mu[,i] <- mu[,i] - sf[i]
      cat("\t", sample_names[i],"\t", signif(sf[i],4),"\n")
    }
    cat("Normalizing counts with factors:\n")
    sf <- apply(log(counts) - rowMeans(log(counts)), 2, function(x){median(x[is.finite(x)])})
    for(i in 1:length(sf)) {
      counts[,i] <- counts[,i]/exp(sf)[i]
      cat("\t", sample_names[i],"\t", signif(exp(sf[i]),4),"\n")
    }
  }
  
  return(list(log_mu_gibbs=mu, mcse=mcse, effective_length=effective_lengths, true_length=true_lengths[,1], unique_hits=unique_hits, counts=counts))
}


