configure_Cap3.R

opt <- list(samplesFile="samples.txt",capFolder="04-Cap3",newblerFolder="03-Assemblies")

cap3_header <- function(output){
  config_output <- c(
    "#!/bin/sh",
    "## Script autogenerated by configure_Cap3.R")
  writeLines(config_output,output)
}

write_out_cap3_assembly <- function(reads,folder,sample,script){
  if (length(reads) != 1 ) stop("Cap3 only accepts a single file")
  call <- "cap3"
  info <- paste(">",file.path(paste(sample,"cap","out",sep=".")),sep=" ")
  fullCall <- paste("nohup",call,reads,info,"&",sep=" ")
  writeLines(paste("cd", file.path(getwd(),folder,sample),sep=" "),script)
  writeLines(fullCall,script)
}

dir.create(opt$capFolder,showWarnings=FALSE,recursive=TRUE)

if(!file.exists(dir(pattern=opt$samplesFile,full.names=TRUE))) stop("Samples file does not exist")
samples <- read.table(opt$samplesFile,sep="",header=T,as.is=T)

zz <- file(file.path(opt$capFolder,"runCap3.sh"), "w")  # open an output file connection
cap3_header(zz)

library(ShortRead)
library(rSFFreader)
for (i in unique(samples$SAMPLE_ID)){
  dir.create(file.path(opt$capFolder,i),showWarnings=FALSE,recursive=TRUE)
  scaffolds <- readFastaQual(file.path(getwd(),opt$newblerFolder,i),fastaPattern="454Scaffolds.fna",qualPattern="454Scaffolds.qual")
  allContigs <- readFastaQual(file.path(getwd(),opt$newblerFolder,i),fastaPattern="454LargeContigs.fna",qualPattern="454LargeContigs.qual")
  joinReads <- c(sread(scaffolds),sread(allContigs))
  joinQual <- c(quality(quality(scaffolds)),quality(quality(allContigs)))
  joinID <- c(as.character(id(scaffolds)),as.character(id(allContigs)))
  joinSRQ <- ShortReadQ(joinReads,joinQual,BStringSet(joinID))
  dups <-  duplicated(joinReads)
  joinSRQ <- joinSRQ[!dups]
  
  writeFastaQual(joinSRQ,file.path(opt$capFolder,i,paste(i,"NewblerContigs",sep="_")))
  write_out_cap3_assembly(reads=paste(i,"NewblerContigs.fasta",sep="_"),folder=opt$capFolder,sample=i,script=zz)         
}

writeLines(paste("cd", file.path(getwd(),opt$capFolder),sep=" "),zz)

close(zz)

Sys.chmod(file.path(opt$capFolder,"runCap3.sh"),"775")