readmmseq <- function( mmseq_files = grep("gene|identical", dir( pattern="\\.mmseq$" ), value=TRUE, invert=TRUE),
		 sample_names = sub("\\.mmseq","", mmseq_files), normalize=TRUE, counts=FALSE) {
	
	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()
	lengths = c()
  n_reads <- 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 )
    if(any(colnames(dat)=="log_mu_gibbs") && any(colnames(dat)=="mcse")) {
      mu <- cbind(mu, dat$log_mu_gibbs)
      mcse <- cbind(mcse, dat$mcse)
    }
    if(i==1) {
      rownames(mu) = rownames(mcse)=dat[,grepl("id", colnames(dat))]
    } else {
      if(! all(rownames(mu) == dat[,1])) {
        stop("Attempting to join incompatible mmseq files (different features).\n")
      }
    }

    if(length(lengths)==0) {
      lengths <- dat[,grepl("length", colnames(dat))]
    }

    ta <- readLines(mmseq_files[i], n=100) # header should be no more than 100 lines
    ta <- strsplit(ta[grep("# Mapped reads", 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)=sample_names
  names(n_reads) <- sample_names
	
  if(!counts) {
    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")
      }
    }
    return(list(log_mu_gibbs=mu, mcse=mcse))
  } else {

    # 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]*lengths)/(1e3*1e6) } )
    counts[is.na(counts)] = 0 # the ones that had no expression in MMSEQ
    
    if(normalize) {
      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(counts)
  }
}
