stringlist2patternandkinship.Rscript

#!/usr/bin/env Rscript
options(error = quote({
	dump.frames(to.file=TRUE, dumpto="Rcoredump")
	load("Rcoredump.rda")
	print(Rcoredump)
	q()
}))

# Start time
start.time = proc.time()[3]
args = commandArgs(trailingOnly=TRUE)

# Process arguments
help = c("stringlist2patternandkinship.Rscript define presence/absence of a kmer list among kmer count files and build kinship matrix",
			"Adapted from Daniel Wilson scripts (2018) kmerlist2pattern.Rscript, patternmerge.Rscript, pattern2kinship.Rscript, kinshipmerge.Rscript",
			"Usage: stringlist2patternandkinship.Rscript task_id p id_file fullkmerlistfile kmercountslistfile analysis_dir output_prefix kmertype software_file [kmer_length=31 mincount=5]")

if(length(args!=0)){
	if(args[1]=="-help" | args[1]=="-h"){
		cat(help,sep="\n")
		q("no")
	}
}

if(length(args)!=9 & length(args)!=11) {
	cat(help,sep="\n")
	cat("Received arguments: ",args,"\n")
	stop("\nIncorrect usage\n")
}

###################################################################################################
## Functions and software paths
###################################################################################################

create_pattern_batch = function(fullkmerlistfile = NULL, p = NULL, t = NULL, stringlist2patternpath = NULL, kmerlist2patternpath = NULL, kmercountslistfile = NULL, kmerlen = NULL, mincount = NULL, kmertype = NULL, output_prefix = NULL, output_dir = NULL){
	
	# Pattern batches directory - create if doesn't exist
	batches_dir = file.path(output_dir,paste0(kmertype,"kmer", kmerlen,"_patternbatches/"))
	if(!dir.exists(batches_dir)){
		cat("Creating temp directory for pattern batches:", batches_dir, "\n")
		dir.create(batches_dir)
	}
	
	# Read kmers
	# Total number of kmers
	n = as.integer(scan(pipe(paste0("zcat ", fullkmerlistfile," | wc -l")), quiet = TRUE))  # SGE added quiet = TRUE
	if(n<1) stop("No kmers found in", fullkmerlistfile)
	# Number of kmers (batch size) per process
	nkmersbatch = get_kmer_batch_numbers(n = n, p = p) # SGE added
	beg = nkmersbatch[t,1] # SGE updated
	end = nkmersbatch[t,2] # SGE updated
	
	output_prefix_batch = paste0(output_prefix,"_",kmertype,kmerlen,".",beg,"-",end) # SGE changed ":" to "-" in out_prefix

	kmersublistfile = paste0(batches_dir, output_prefix_batch,".",t,".temp_kmerlist.txt.gz")
	cat("Creating temp file:",kmersublistfile,"\n")
	cmd = paste0("zcat ", fullkmerlistfile," | head -n ",end," | tail -n ",end-beg+1," | gzip -c > ",kmersublistfile)
	stopifnot(system(cmd)==0)
	softwarepath = ifelse(kmertype=="protein", stringlist2patternpath, kmerlist2patternpath)
	cmd = paste(softwarepath, kmersublistfile, kmercountslistfile, paste0(batches_dir, output_prefix_batch), kmerlen, mincount)
	stopifnot(system(cmd)==0)

	cmd = paste("rm",kmersublistfile)
	stopifnot(system(cmd)==0)

	# Check that files have been created and are not empty
	outfiles = paste0(batches_dir, output_prefix_batch,c(".patternKey.txt.gz",".patternIndex.txt.gz"))
	outfiles_size = sapply(outfiles, function(x) as.numeric(system(paste0("ls -l ",x," | cut -d ' ' -f5"), intern = T)), USE.NAMES = F) # SGE added this - could set higher than zero
	if(any(outfiles_size==0)) stop("One or more task_id ",t, " ", output_prefix_batch, " patternKey, patternKeySize or patternIndex files are empty")
	
	
	# Write output file prefix to output files
	cat(paste0("Output file prefix: ", output_prefix_batch),"\n")
	outfile_completed = paste0(batches_dir, output_prefix_batch,".patternbatch.completed.txt")
	stopifnot(system2("/bin/bash",paste0("-c 'touch ",outfile_completed,"'"), wait = TRUE)==0)
	
	return(list("nkmersbatch" = nkmersbatch, "output_prefix_batch" = output_prefix_batch, "batches_dir" = batches_dir))

}

get_kmer_batch_numbers = function(n = NULL, p = NULL){
	
	# Number of kmers (batch size) per process
	b = as.integer(ceiling(n/p))

	# Define kmers to process # SGE updated to output for all batches
	# beg = (t-1)*b+1
	# end = min(t*b,n)
	beg = sapply(1:p, function(t) (t-1)*b+1, USE.NAMES = F)
	end = sapply(1:p, function(t) min(t*b,n), USE.NAMES = F)
	
	# output_prefix_batch = paste0(output_prefix,".",beg,"-",end) # SGE changed ":" to "-" in out_prefix
	return(cbind(beg, end))
	
}

# Function created from DJW patternmerge.Rscript
merge_pattern_batches_parameters = function(nkmersbatch = NULL, p = NULL, t = NULL, output_dir = NULL, output_prefix = NULL, kmertype = NULL, kmerlen = NULL){
	
	# SGE changed from patternmerge to use the number of total kmers to get the batch numbers used to create the patterns, rather than search for the files in the directory, incase has been run multiple times with the same prefix but different sample sizes
	output_prefix_batches = apply(nkmersbatch, 1, function(b, prefix) paste0(prefix, ".", b[1], "-", b[2]), prefix = paste0(output_prefix, "_", kmertype, kmerlen))
	batches_dir = file.path(output_dir,paste0(kmertype,"kmer", kmerlen,"_patternbatches/"))
	infiles = paste0(batches_dir, output_prefix_batches, ".patternKey.txt.gz")
	infiles_completed = paste0(batches_dir, output_prefix_batches,".patternbatch.completed.txt")
	nattempts = 0
	while(!all(file.exists(infiles_completed)) | !all(file.exists(infiles))) {
		nattempts=nattempts+1
		if(nattempts>100) stop("Could not find files",infiles_completed)
		Sys.sleep(60) # SGE changed from 1 to 60
	}
	
	n = as.integer(length(infiles))
	if(!all(1+nkmersbatch[,2][-n]==nkmersbatch[,1][-1])) stop("Input files do not have consecutive ranges")
	
	# Determine remaining parameters
	b = as.integer(ceiling(n/p))
	if(n==1) {
		imax = as.integer(1)
	} else {
		if(b==1) stop("Pattern merging - cannot have batchsize = 1. Try p < n/2")
		if(p != as.integer(ceiling(n/b))) {
			cat("Warning: pattern merging - adjusting number of processes to equal ceiling(n/b)\n")
			p = as.integer(ceiling(n/b))
		}
		imax = as.integer(ceiling(log(n)/log(b)))
	}

	return(list("n" = n, "b" = b, "imax" = imax, "infiles" = infiles))
	
}

# Function created from DJW patternmerge.Rscript
merge_patterns = function(t = NULL, n = NULL, b = NULL, p = NULL, imax = NULL, prefix = NULL, files = NULL, patternmergepath = NULL, output_dir = NULL, kmertype = NULL, kmerlen = NULL, batches_dir = NULL){

	# Merge
	i = 0
	while(TRUE) {
    	i=i+1
	    if((n==1 & i>1) | !((t %% b^(i-1))==0 | (t==p & i<=imax))) break()
    	cat("t =",t,"i =",i,"\n")
	    if(i==1) {
	        # First round: merge source files
    	    		outfile_prefix = paste0(prefix,"_",kmertype,kmerlen,".patternmerge.j.",i,".",t) # SGE updated
        		outfile_patternKey = paste0(batches_dir, outfile_prefix,".patternKey.txt.gz")
	    		outfile_patternKeySize = paste0(batches_dir, outfile_prefix,".patternKeySize.txt")
    	    		outfile_patternIndex = paste0(batches_dir, outfile_prefix,".patternIndex.txt.gz")
			outfile_completed = paste0(batches_dir, outfile_prefix,".patternbatch.completed.txt") # SGE added
       		beg = b*(t-1)+1
	        end = min(b*t, n)
        	if(end<beg) stop("Problem with input arguments, please check")
        	infiles = files[beg:end]
        	# *** [This should already be dealt with in the parameters function - maybe remove] ***
	        nattempts = 0
    	    while(!all(file.exists(infiles))) {
        	    nattempts=nattempts+1
            	if(nattempts>1000) stop("Could not find files",infiles)
	            Sys.sleep(1)
    	    }

			# Check files aren't empty
			infiles_size = sapply(infiles, function(x) as.numeric(system(paste0("ls -l ",x," | cut -d ' ' -f5"), intern = T)), USE.NAMES = F) # SGE added
			if(any(infiles_size==0)) stop("One or more file size is zero ",infiles,"\n") # SGE added - could make higher than zero
        	# Create patternKeySize files
        	for(j in 1:length(infiles)) {
				stopifnot(system(paste0("zcat ",infiles[j]," | wc -l > ",gsub(".patternKey.txt.gz",".patternKeySize.txt",infiles[j],fixed=TRUE)))==0)
            	cat("Created ",gsub(".patternKey.txt.gz",".patternKeySize.txt",infiles[j],fixed=TRUE),"\n")
	        }

    	    # Merge the patternKey files and redefine the patternKeySize and patternIndex files
        	if(length(infiles)==1) {
				stopifnot(system(paste0("cp ",infiles[1]," ",outfile_patternKey))==0)
				stopifnot(system(paste0("cp ",gsub(".patternKey.txt.gz",".patternKeySize.txt",infiles[1],fixed=TRUE)," ",outfile_patternKeySize))==0)
				stopifnot(system(paste0("cp ",gsub(".patternKey.txt.gz",".patternIndex.txt.gz",infiles[1],fixed=TRUE)," ",outfile_patternIndex))==0)
				stopifnot(system2("/bin/bash",paste0("-c 'touch ",outfile_completed,"'"))==0)

        	} else if(length(infiles)==2) {
            	prefixA = gsub(".patternKey.txt.gz","",infiles[1],fixed=TRUE)
            	prefixB = gsub(".patternKey.txt.gz","",infiles[2],fixed=TRUE)
				stopifnot(system(paste0(patternmergepath, " ", prefixA, " ", prefixB, " ", batches_dir,outfile_prefix))==0)
				stopifnot(system2("/bin/bash",paste0("-c 'touch ",outfile_completed,"'"))==0)

        	} else {
            	# Create temporary storage filenames
            	tmpfile_prefix = paste0(batches_dir, "tmpfile",1:2,".",outfile_prefix)
            	tmpfile_patternKey = paste0(tmpfile_prefix,".patternKey.txt.gz")
            	tmpfile_patternKeySize = paste0(tmpfile_prefix,".patternKeySize.txt")
            	tmpfile_patternIndex = paste0(tmpfile_prefix,".patternIndex.txt.gz")

            	# Do the merges
            	for(j in 2:length(infiles)) {

                	if(j==2) {
                    	prefixA = gsub(".patternKey.txt.gz","",infiles[1],fixed=TRUE)
                    	prefixB = gsub(".patternKey.txt.gz","",infiles[2],fixed=TRUE)
						stopifnot(system(paste0(patternmergepath, " ",prefixA," ",prefixB," ",tmpfile_prefix[1+(j%%2)]))==0)
                	} else {
                    	prefixA = tmpfile_prefix[1+((j-1)%%2)]
                    	prefixB = gsub(".patternKey.txt.gz","",infiles[j],fixed=TRUE)
						stopifnot(system(paste0(patternmergepath, " ",prefixA," ",prefixB," ",tmpfile_prefix[1+(j%%2)]))==0)
                	}
            	}
            	# Finally move the tmpfiles into place
            	# Behaviour: j equals the last value in the list
				stopifnot(system(paste0("mv ",tmpfile_patternKey[1+(j%%2)]," ",outfile_patternKey))==0)
				stopifnot(system(paste0("mv ",tmpfile_patternKeySize[1+(j%%2)]," ",outfile_patternKeySize))==0)
				stopifnot(system(paste0("mv ",tmpfile_patternIndex[1+(j%%2)]," ",outfile_patternIndex))==0)
            	# Delete the other tmpfiles
            	unlink(tmpfile_patternKey[1+((j-1)%%2)])
            	unlink(tmpfile_patternKeySize[1+((j-1)%%2)])
            	unlink(tmpfile_patternIndex[1+((j-1)%%2)])
				# Create file stating that the round is completed
				stopifnot(system2("/bin/bash",paste0("-c 'touch ",outfile_completed,"'"))==0)
        	}
    		
    	} else {
        	# Subsequent rounds: merge merged files
        	te = as.integer(ceiling(t/(b^(i-1)))*b^(i-1))
        	outfile_prefix = paste0(prefix,"_",kmertype,kmerlen,".patternmerge.j.",i,".",te) # SGE updated
        	outfile_patternKey = paste0(batches_dir, outfile_prefix,".patternKey.txt.gz")
        	outfile_patternKeySize = paste0(batches_dir, outfile_prefix,".patternKeySize.txt")
        	outfile_patternIndex = paste0(batches_dir, outfile_prefix,".patternIndex.txt.gz")
		outfile_completed = paste0(batches_dir, outfile_prefix,".patternbatch.completed.txt") # SGE added
        	beg = te-b^(i-1)+b^(i-2)
	        end = min(te,as.integer(ceiling(t/b^(i-2))*b^(i-2)))
        	if(end<beg) stop("Problem with input arguments, please check")
        	inc = b^(i-2)
        	
        	
        	infiles = paste0(batches_dir, prefix,"_",kmertype,kmerlen,".patternmerge.j.",i-1,".",seq(from=beg,to=end,by=inc),".patternKey.txt.gz")
		infiles_completed = paste0(batches_dir, prefix,"_",kmertype,kmerlen,".patternmerge.j.",i-1,".",seq(from=beg,to=end,by=inc),".patternbatch.completed.txt") # SGE added
		nattempts = 0
        	while(!all(file.exists(infiles_completed)) | !all(file.exists(infiles))) { # SGE updated to check completed file existence as well
            	nattempts=nattempts+1
            	if(nattempts>3600) stop("Could not find files",infiles_completed)
            	Sys.sleep(1)
        	}

			infiles_size = sapply(infiles, function(x) as.numeric(system(paste0("ls -l ",x," | cut -d ' ' -f5"), intern = T))) # SGE added
			if(any(infiles_size==0)) stop("One or more file size is zero ",paste(infiles, collapse = " "),"\n") # SGE added - could make higher than zero
        	# Merge the patternKey files and redefine the patternKeySize and patternIndex files
        	if(length(infiles)==1) {
				stopifnot(system(paste0("cp ",infiles[1]," ",outfile_patternKey))==0)
				stopifnot(system(paste0("cp ",gsub(".patternKey.txt.gz",".patternKeySize.txt",infiles[1],fixed=TRUE)," ",outfile_patternKeySize))==0)
				stopifnot(system(paste0("cp ",gsub(".patternKey.txt.gz",".patternIndex.txt.gz",infiles[1],fixed=TRUE)," ",outfile_patternIndex))==0)
				stopifnot(system2("/bin/bash",paste0("-c 'touch ",outfile_completed,"'"))==0)
        	} else if(length(infiles)==2) {
            	prefixA = gsub(".patternKey.txt.gz","",infiles[1],fixed=TRUE)
            	prefixB = gsub(".patternKey.txt.gz","",infiles[2],fixed=TRUE)
				stopifnot(system(paste0(patternmergepath, " ",prefixA," ",prefixB," ", batches_dir, outfile_prefix))==0)
				stopifnot(system2("/bin/bash",paste0("-c 'touch ",outfile_completed,"'"))==0)
        	} else {
            	# Create temporary storage filenames
            	ntmp = 2
            	if(ntmp<2) stop("ntmp must be at least 2")
            	tmpfile_prefix = paste0(batches_dir, "tmpfile",1:ntmp,".",outfile_prefix)
            	tmpfile_patternKey = paste0(tmpfile_prefix,".patternKey.txt.gz")
            	tmpfile_patternKeySize = paste0(tmpfile_prefix,".patternKeySize.txt")
            	tmpfile_patternIndex = paste0(tmpfile_prefix,".patternIndex.txt.gz")
            	# Do the merges
            	for(j in 2:length(infiles)) {
                	if(j==2) {
                    	prefixA = gsub(".patternKey.txt.gz","",infiles[1],fixed=TRUE)
                    	prefixB = gsub(".patternKey.txt.gz","",infiles[2],fixed=TRUE)
						stopifnot(system(paste0(patternmergepath, " ",prefixA," ",prefixB," ",tmpfile_prefix[1+((j-1)%%ntmp)]))==0)
                	} else {
                    	prefixA = tmpfile_prefix[1+((j-2)%%ntmp)]
                    	prefixB = gsub(".patternKey.txt.gz","",infiles[j],fixed=TRUE)
						stopifnot(system(paste0(patternmergepath, " ",prefixA," ",prefixB," ",tmpfile_prefix[1+((j-1)%%ntmp)]))==0)
                    	# Delete the previous tmpfile
                    	unlink(tmpfile_patternKey[1+((j-2)%%ntmp)])
                    	unlink(tmpfile_patternKeySize[1+((j-2)%%ntmp)])
                    	unlink(tmpfile_patternIndex[1+((j-2)%%ntmp)])
                	}
            	}
            	# Finally move the tmpfiles into place
            	# Behaviour: j equals the last value in the list
				stopifnot(system(paste0("mv ",tmpfile_patternKey[1+((j-1)%%ntmp)]," ",outfile_patternKey))==0)
				stopifnot(system(paste0("mv ",tmpfile_patternKeySize[1+((j-1)%%ntmp)]," ",outfile_patternKeySize))==0)
				stopifnot(system(paste0("mv ",tmpfile_patternIndex[1+((j-1)%%ntmp)]," ",outfile_patternIndex))==0)
				stopifnot(system2("/bin/bash",paste0("-c 'touch ",outfile_completed,"'"))==0)
        	}
    	}
	}

	if(t==p) {
		create_final_merged_pattern_files(outfile_patternKey = outfile_patternKey, outfile_patternKeySize = outfile_patternKeySize, outfile_patternIndex = outfile_patternIndex, outfile_prefix = outfile_prefix, prefix = prefix, output_dir = output_dir, kmertype = kmertype, kmerlen = kmerlen, files = files, batches_dir = batches_dir)
	}

}

create_final_merged_pattern_files = function(outfile_patternKey = NULL, outfile_patternKeySize = NULL, outfile_patternIndex = NULL, outfile_prefix = NULL, prefix = NULL, output_dir = NULL, kmertype = NULL, kmerlen = NULL, files = NULL, batches_dir = NULL){
	
	# Check final files aren't empty # SGE added
	outfiles_size = sapply(c(outfile_patternKey, outfile_patternKeySize, outfile_patternIndex), function(x) as.numeric(system(paste0("ls -l ",x," | cut -d ' ' -f5"), intern = T)), USE.NAMES = F) # SGE added - could set higher than zero
	if(any(outfiles_size==0)) stop("One or more task_id ",t, " final ",outfile_prefix, " patternKey, patternKeySize or patternIndex files are empty")
	final_file_prefix = paste0(output_dir, prefix, "_", kmertype, kmerlen)
	stopifnot(system(paste0("mv ",outfile_patternKey," ", final_file_prefix, ".patternmerge.patternKey.txt.gz"))==0)
	stopifnot(system(paste0("mv ",outfile_patternKeySize," ", final_file_prefix, ".patternmerge.patternKeySize.txt"))==0)
	stopifnot(system(paste0("mv ",outfile_patternIndex," ", final_file_prefix, ".patternmerge.patternIndex.txt.gz"))==0)
	
	cat(paste0("Final pattern files: ", final_file_prefix, ".patternmerge.patternKey.txt.gz ", final_file_prefix, ".patternmerge.patternKeySize.txt ", final_file_prefix, ".patternmerge.patternIndex.txt.gz "),"\n")
	
	# Remove all completed files - SGE added
	# completed_files = dir(pattern=glob2rx(paste0(batches_dir, prefix,"_",kmertype,kmerlen,"*.patternbatch.completed.txt")))
	completed_files = system(paste0("ls ", paste0(batches_dir, prefix,"_",kmertype,kmerlen,"*.patternbatch.completed.txt")), intern = T)
	# Remove all intermediate pattern files - SGE added
	cat("Removing intermediate files","\n")
	# cat("Length completed files:", length(completed_files), "\n")
	cmd = paste("rm",paste(completed_files,collapse=" "))
	stopifnot(system2("/bin/bash",paste0("-c '",cmd,"'"), wait = TRUE)==0)
	
	# pattern_files = dir(pattern=glob2rx(paste0(batches_dir, prefix,"_",kmertype,kmerlen,".patternmerge.j*.patternKey.txt.gz")))
	# pattern_files = c(pattern_files, dir(pattern=glob2rx(paste0(batches_dir, prefix,"_",kmertype,kmerlen,".patternmerge.j*.patternKeySize.txt"))))
	# pattern_files = c(pattern_files, dir(pattern=glob2rx(paste0(batches_dir, prefix,"_",kmertype,kmerlen,".patternmerge.j*.patternIndex.txt.gz"))))
	pattern_files = system(paste0("ls ", paste0(batches_dir, prefix,"_",kmertype,kmerlen,".patternmerge.j*.patternKey.txt.gz")), intern = T)
	pattern_files = c(pattern_files, system(paste0("ls ", paste0(batches_dir, prefix,"_",kmertype,kmerlen,".patternmerge.j*.patternKeySize.txt")), intern = T))
	pattern_files = c(pattern_files, system(paste0("ls ", paste0(batches_dir, prefix,"_",kmertype,kmerlen,".patternmerge.j*.patternIndex.txt.gz")), intern = T))

	# cmd = paste("rm",paste(pattern_files,collapse=" "))
	# pattern_files = c(pattern_files, dir(pattern=glob2rx(paste0(prefix,"_",kmertype,kmerlen,"*-*.patternKeySize.txt"))))
	pattern_files = c(pattern_files, gsub(".patternKey.txt.gz",".patternKeySize.txt",files,fixed=TRUE)) # SGE replaced above
	# cat("Length pattern files:", length(pattern_files), "\n")
	cmd = paste("rm",paste(pattern_files,collapse=" "))
	stopifnot(system2("/bin/bash",paste0("-c '",cmd,"'"), wait = TRUE)==0)

	
}

# Function created from DJW pattern2kinship.Rscript
create_kinship_batch = function(batch_prefix = NULL, batches_dir = NULL, nkmersbatch = NULL, patterncountspath = NULL, pattern2kinshippath = NULL, kmertype = NULL, kmerlen = NULL){
	
	# Calculate patternCounts
	# file_prefix = gsub(".patternKey.txt.gz","",files[od][t],fixed=TRUE)
	cmd = paste0(patterncountspath, " ", batches_dir, batch_prefix) # SGE changed file_prefix to batches_dir, batch_prefix
	stopifnot(system(cmd)==0)

	# Compute kinship matrix
	cmd = paste0(pattern2kinshippath, " ", batches_dir, batch_prefix) # SGE changed file_prefix to batches_dir, batch_prefix
	stopifnot(system(cmd)==0)
	
	# Check kinship file not empty # SGE added
	if(as.numeric(system(paste0("ls -l ", batches_dir, batch_prefix, ".kinship.txt.gz | cut -d ' ' -f5"), intern = T))==0) 	stop(paste0(batches_dir,batch_prefix, " kinship matrix file empty"),"\n")

	# Output kinship matrix weight (i.e. total count)
	# write(file.end[t]-file.beg[t]+1,paste0(file_prefix,".kinshipWeight.txt"))
	write(nkmersbatch[2]-nkmersbatch[1]+1,paste0(batches_dir, batch_prefix, ".kinshipWeight.txt")) # SGE changed from above
	# Check kinship weight file not empty # SGE added
	if(as.numeric(system(paste0("ls -l ",batches_dir, batch_prefix, ".kinshipWeight.txt | cut -d ' ' -f5"), intern = T))==0) stop(paste0(batches_dir,batch_prefix, " kinship weight file empty"),"\n")
	
	cat("Created kinship matrix batch:",paste0(batches_dir, batch_prefix, ".kinship.txt.gz"),"\n")
	# Write output file prefix to output files
	outfile_completed = paste0(batches_dir, batch_prefix, ".kinship.completed.txt")
	stopifnot(system2("/bin/bash",paste0("-c 'touch ",outfile_completed,"'"), wait = TRUE)==0)
	
}

# Function created from DJW kinshipmerge.Rscript
merge_kinship_matrices = function(t = NULL, n = NULL, b = NULL, p = NULL, imax = NULL, files = NULL, prefix = NULL, nkmersbatch = NULL, output_dir = NULL, kmertype = NULL, kmerlen = NULL, batches_dir = NULL){
	
	output_prefix_batches = apply(nkmersbatch, 1, function(b, batches_dir, prefix, kmertype, kmerlen) paste0(batches_dir, prefix, "_", kmertype, kmerlen, ".", b[1], "-", b[2]), batches_dir = batches_dir, prefix = prefix, kmertype = kmertype, kmerlen = kmerlen)
	# batches_dir = file.path(output_dir,paste0(kmertype,"kmer", kmerlen,"_patternbatches/"))
	
	# Merge
	i = 0
	nsamp = 0
	while(TRUE) {
    	i=i+1
	if((n==1 & i>1) | !((t %% b^(i-1))==0 | (t==p & i<=imax))) break()
    	cat("t=",t,"i=",i,"\n")
    	if(i==1) {
        	# First round: merge initial kinship matrices
        	kinship = matrix(as.double(0),0,0)
        	kinship.weight = as.integer(0)
        	# Define input/output files
        	outfile_prefix = paste0(batches_dir, prefix,"_", kmertype, kmerlen, ".kinshipmerge.j.",i,".",t) # SGE updated
        	beg = b*(t-1)+1
        	end = min(b*t, n)
        	if(end<beg) stop("Problem with input arguments, please check")
        	infiles = files[beg:end] # SGE removed [od] after files
		infiles_completed = paste0(output_prefix_batches, ".kinship.completed.txt")[beg:end] # SGE added
        	nattempts = 0
      	# Check input files exist and are not empty
		infiles_size = sapply(infiles, function(x) as.numeric(system(paste0("ls -l ",x," | cut -d ' ' -f5"), intern = T)), USE.NAMES = F) # SGE added - could set higher than zero, an appropriate number for completed
        	while(!all(file.exists(infiles)) | !all(infiles_size>0) | !all(file.exists(infiles_completed))) { # SGE added file size check
            	nattempts=nattempts+1
            	if(nattempts>100) stop("Could not find files",infiles)
            	Sys.sleep(60)
        	}
        	# Read one line from the first kinship matrix
        	kin1 = scan(infiles[1],what=double(0),nlines=1, quiet = TRUE) # SGE added quiet = TRUE
        	nsamp = length(kin1)
        	# Allocate memory
        	kinship = matrix(as.double(0),nsamp,nsamp)
        	# Read each file and augment the kinship matrix
        	for(j in 1:length(infiles)) {
            	kinship.j = matrix(scan(infiles[j],what=double(0),nlines=nsamp, quiet = TRUE),nsamp,byrow=TRUE) # SGE added quiet = TRUE
            	if(ncol(kinship.j)!=nsamp) {
                	stop("Expected ",nsamp," columns in ",infiles[j])
            	}
            	if(any(is.na(kinship.j))) {
                	stop("Found NAs in ",infiles[j])
            	}
            	weightfile = paste0(gsub(".kinship.txt.gz","",infiles[j],fixed=TRUE),".kinshipWeight.txt")
            	kinship.weight.j = scan(weightfile,what=integer(0), quiet = TRUE) # SGE added quiet = TRUE
            	if(length(kinship.weight.j)!=1) {
                	stop("Expected one entry in kinship weight file ",weightfile)
            	}
            	# Augment the kinship matrix and its weights
            	kinship = kinship + kinship.weight.j*kinship.j
            	kinship.weight = kinship.weight + kinship.weight.j
        	}
        	# Save the kinship matrix and its weight
        	kinship = kinship / as.double(kinship.weight)
        	write(sprintf("%.17g",kinship),gzfile(paste0(outfile_prefix,".kinship.txt.gz")),ncol=nsamp)
        	write(sprintf("%d",kinship.weight),paste0(outfile_prefix,".kinshipWeight.txt"))
			stopifnot(system2("/bin/bash",paste0("-c 'touch ",outfile_prefix,".kinship.completed.txt'"), wait = TRUE)==0)
    	} else {
        	# Subsequent rounds: merge merged kinship matrices
        	kinship = matrix(as.double(0),nsamp,nsamp)
        	kinship.weight = as.integer(0)
        	# Define input/output files
        	te = as.integer(ceiling(t/(b^(i-1)))*b^(i-1))
        	outfile_prefix = paste0(batches_dir, prefix,"_", kmertype, kmerlen, ".kinshipmerge.j.",i,".",te) # SGE updated
        	outfile_kinship = paste0(outfile_prefix,".kinship.txt.gz")
        	outfile_kinshipWeight = paste0(outfile_prefix,".kinshipWeight.txt")
		outfile_completed = paste0(outfile_prefix,".kinship.completed.txt") # SGE added
        	beg = te-b^(i-1)+b^(i-2)
        	end = min(te,as.integer(ceiling(t/b^(i-2))*b^(i-2)))
        	if(end<beg) stop("Problem with input arguments, please check")
        	inc = b^(i-2)
        	infiles = paste0(batches_dir, prefix,"_", kmertype, kmerlen, ".kinshipmerge.j.",i-1,".",seq(from=beg,to=end,by=inc),".kinship.txt.gz") # SGE updated
        	infiles.weights = paste0(batches_dir, prefix,"_", kmertype, kmerlen, ".kinshipmerge.j.",i-1,".",seq(from=beg,to=end,by=inc),".kinshipWeight.txt") # SGE updated
		infiles.completed = paste0(batches_dir, prefix,"_", kmertype, kmerlen, ".kinshipmerge.j.",i-1,".",seq(from=beg,to=end,by=inc),".kinship.completed.txt") # SGE added
        	nattempts = 0
        	while(!all(file.exists(infiles)) | !all(file.exists(infiles.weights))| !all(file.exists(infiles.completed))) { # SGE added file size check
            	nattempts=nattempts+1
            	if(nattempts>3600) stop("Could not find files",infiles)
            	Sys.sleep(1)
        	}
			infiles_size = sapply(infiles, function(x) as.numeric(system(paste0("ls -l ",x," | cut -d ' ' -f5"), intern = T)), USE.NAMES = F) # SGE added
			if(any(infiles_size==0)) stop("One or more file size is zero ",paste(infiles, collapse = " "),"\n") # SGE added - could make higher than zero
			infiles.weights_size = sapply(infiles.weights, function(x) as.numeric(system(paste0("ls -l ",x," | cut -d ' ' -f5"), intern = T)), USE.NAMES = F) # SGE added
			if(any(infiles.weights_size==0)) stop("One or more file size is zero ",paste(infiles.weights, collapse = " "),"\n") # SGE added - could make higher than zero
        	# Read each file and augment the kinship matrix
        	for(j in 1:length(infiles)) {
        		kinship.j = matrix(scan(infiles[j],what=double(0),nlines=nsamp, quiet = TRUE),nsamp,byrow=TRUE) # SGE added quiet = TRUE
            	if(ncol(kinship.j)!=nsamp) {
                	stop("Expected ",nsamp," columns in ",infiles[j])
	            }
    	        if(any(is.na(kinship.j))) {
        	        stop("Found NAs in ",infiles[j])
            	}
            	kinship.weight.j = scan(infiles.weights[j],what=integer(0), quiet = TRUE) # SGE added quiet = TRUE
	            if(length(kinship.weight.j)!=1) {
    	            stop("Expected one entry in kinship weight file ",weightfile)
        	    }
            	# Augment the kinship matrix and its weights
	            kinship = kinship + kinship.weight.j*kinship.j
    	        kinship.weight = kinship.weight + kinship.weight.j
			}
			# Save the kinship matrix and its weight
			kinship = kinship / as.double(kinship.weight)
			write(sprintf("%.17g",kinship),gzfile(outfile_kinship),ncol=nsamp)
			write(sprintf("%d",kinship.weight),outfile_kinshipWeight)
			stopifnot(system2("/bin/bash",paste0("-c 'touch ",outfile_completed,"'"), wait = TRUE)==0)
		}
	}

	if(t==p) {
		# Define input/output files: bug fix DJW 20220522
		i = i-1
		te = as.integer(ceiling(t/(b^(i-1)))*b^(i-1))
		outfile_prefix = paste0(batches_dir, prefix,"_", kmertype, kmerlen, ".kinshipmerge.j.",i,".",te) # SGE updated
		outfile_kinship = paste0(outfile_prefix,".kinship.txt.gz")
		outfile_kinshipWeight = paste0(outfile_prefix,".kinshipWeight.txt")
		create_final_kinship_file(outfile_kinship = outfile_kinship, outfile_kinshipWeight = outfile_kinshipWeight, t = t, outfile_prefix = outfile_prefix, prefix = prefix, output_dir = output_dir, kmertype = kmertype, kmerlen = kmerlen, files = files, batches_dir = batches_dir)
	}

}

# Function created from DJW kinshipmerge.Rscript
create_final_kinship_file = function(outfile_kinship = NULL, outfile_kinshipWeight = NULL, t = NULL, outfile_prefix = NULL, output_dir = NULL, prefix = NULL, kmertype = NULL, kmerlen = NULL, files = NULL, batches_dir = NULL){
	
	outfiles_size = sapply(c(outfile_kinship, outfile_kinshipWeight), function(x) as.numeric(system(paste0("ls -l ",x," | cut -d ' ' -f5"), intern = T)), USE.NAMES = F) # SGE added
	if(any(outfiles_size==0)) stop("One or more task_id ",t, " ", outfile_prefix, " kinship or kinship weight files are empty") # SGE corrected out_prefix to outfile_prefix
	final_file_prefix = paste0(output_dir, prefix, "_", kmertype, kmerlen)
	stopifnot(system(paste0("mv ",outfile_kinship," ", final_file_prefix, ".kinshipmerge.kinship.txt.gz"))==0)
	stopifnot(system(paste0("mv ",outfile_kinshipWeight," ", final_file_prefix, ".kinshipmerge.kinshipWeight.txt"))==0)
	
	cat(paste0("Final kinship files: ", final_file_prefix, ".kinshipmerge.kinship.txt.gz ", final_file_prefix, ".kinshipmerge.kinshipWeight.txt"),"\n")

	
	# Remove all completed files - SGE added
	# completed_files = dir(pattern=glob2rx(paste0(batches_dir, prefix,"_",kmertype,kmerlen,"*.kinship.completed.txt")))
	completed_files = system(paste0("ls ", paste0(batches_dir, prefix,"_",kmertype,kmerlen,"*.kinship.completed.txt")), intern = T)
	cmd = paste("rm",paste(completed_files,collapse=" "))
	stopifnot(system2("/bin/bash",paste0("-c '",cmd,"'"), wait = TRUE)==0)

	# Remove all intermediate kinship files - SGE added
	# Remove those in pwd
	# kinship_files = dir(pattern=glob2rx(paste0(batches_dir, prefix,"_",kmertype,kmerlen,".kinshipmerge.j*.kinship.txt.gz")))
	# kinship_files = c(kinship_files, dir(pattern=glob2rx(paste0(batches_dir, prefix,"_",kmertype,kmerlen,".kinshipmerge.j*.kinshipWeight.txt"))))
	kinship_files = system(paste0("ls ", paste0(batches_dir, prefix,"_",kmertype,kmerlen,".kinshipmerge.j*.kinship.txt.gz")), intern = T)
	kinship_files = c(kinship_files, system(paste0("ls ", paste0(batches_dir, prefix,"_",kmertype,kmerlen,".kinshipmerge.j*.kinshipWeight.txt")), intern = T))
	# Remove those in pattern subdirectory
	kinship_files = c(kinship_files, files)
	kinship_files = c(kinship_files, gsub(".kinship.txt.gz",".patternCounts.txt.gz", files, fixed = TRUE))
	kinship_files = c(kinship_files, gsub(".kinship.txt.gz",".kinshipWeight.txt", files, fixed = TRUE))
	cat("Removing intermediate files","\n")
	# kinship_files = c(kinship_files, dir(pattern=glob2rx(paste0(prefix,"_",kmertype,kmerlen,".*-*.kinshipWeight.txt"))))
	# kinship_files = c(kinship_files, dir(pattern=glob2rx(paste0(prefix,"_",kmertype,kmerlen,".*-*.patternCounts.txt.gz"))))
	# kinship_files = c(kinship_files, dir(pattern=glob2rx(paste0(prefix,"_",kmertype,kmerlen,".*-*.kinship.txt.gz"))))
	cmd = paste("rm",paste(kinship_files,collapse=" "))
	stopifnot(system2("/bin/bash",paste0("-c '",cmd,"'"), wait = TRUE)==0)
	
}


###################################################################################################


# Initialize variables
t = as.integer(args[1])
p = as.integer(args[2])
id_file = as.character(args[3])
fullkmerlistfile = as.character(args[4])
kmercountslistfile = as.character(args[5])
output_dir = as.character(args[6])
output_prefix = as.character(args[7])
kmertype = tolower(as.character(args[8]))
software_file = as.character(args[9])
if(length(args)>9) {
    kmerlen = as.integer(args[10])
    mincount = as.integer(args[11])
} else {
    kmerlen = as.integer(31)
    mincount = as.integer(5)
}

# Check input arguments
if(is.na(p)) stop("Error: p must be an integer","\n")
if(!file.exists(id_file)) stop("Error: sample ID file doesn't exist","\n")
if(!file.exists(fullkmerlistfile)) stop("Error: full kmer list file doesn't exist","\n")
if(!file.exists(kmercountslistfile)) stop("Error: kmer counts list file doesn't exist","\n")
if(!file.exists(output_dir)) stop("Error: output directory doesn't exist","\n")
if(unlist(strsplit(output_dir,""))[length(unlist(strsplit(output_dir,"")))]!="/") output_dir = paste0(output_dir, "/")
if(kmertype!="protein" & kmertype!="nucleotide") stop("Error: kmer type must be either 'protein' or 'nucleotide'","\n")
if(is.na(kmerlen)) stop("Error: kmer length must be an integer","\n")
if(is.na(mincount)) stop("Error: min count must be an integer","\n")


# Read in software file
software_paths = read.table(software_file, h = T, sep = "\t", quote = "")
required_software = c("R", "scriptpath")
if(any(is.na(match(required_software, as.character(software_paths$name))))) stop(paste0("Error: missing required software path in the software file - requires ",paste(required_software, collapse = ", ")),"\n")

Rpath = as.character(software_paths$path)[which(toupper(as.character(software_paths$name))=="R")]
Rscriptpath = paste0(Rpath, "script")
if(!file.exists(Rscriptpath)) stop("Error: Rscript path",Rscriptpath,"doesn't exist","\n") # SGE added

script_location = as.character(software_paths$path)[which(tolower(as.character(software_paths$name))=="scriptpath")]
if(!dir.exists(script_location)) stop("Error: script location directory specified in the software paths file doesn't exist","\n")

stringlist2patternpath = file.path(script_location, "stringlist2pattern")
if(!file.exists(stringlist2patternpath)) stop("Error: stringlist2pattern path doesn't exist - check pipeline script location in the software file","\n")

kmerlist2patternpath = file.path(script_location, "kmerlist2pattern")
if(!file.exists(kmerlist2patternpath)) stop("Error: kmerlist2pattern path doesn't exist - check pipeline script location in the software file","\n")

patternmergepath = file.path(script_location, "patternmerge")
if(!file.exists(patternmergepath)) stop("Error: patternmerge path doesn't exist - check pipeline script location in the software file","\n")

patterncountspath = file.path(script_location, "patterncounts")
if(!file.exists(patterncountspath)) stop("Error: patterncounts path doesn't exist - check pipeline script location in the software file","\n")

pattern2kinshippath = file.path(script_location, "pattern2kinship")
if(!file.exists(pattern2kinshippath)) stop("Error: pattern2kinship path doesn't exist - check pipeline script location in the software file","\n")

pattern2presencecountRscript = file.path(script_location, "pattern2presencecount.Rscript")
if(!file.exists(pattern2kinshippath)) stop("Error: pattern2presencecount.Rscript path doesn't exist - check pipeline script location in the software file","\n")







# Report variables
cat("#############################################", "\n")
# Report variables
cat("Running on host: ",system("hostname", intern=TRUE),"\n")
cat("Command line arguments","\n")
cat(args, "\n\n")
cat("Parameters:","\n")
cat("task_id:", t, "\n")
cat("p:", p,"\n")
cat("ID file path:", id_file,"\n")
cat("Full kmer list file:", fullkmerlistfile,"\n")
cat("Kmer counts list file:", kmercountslistfile,"\n")
cat("Analysis directory:", output_dir,"\n")
cat("Output prefix:", output_prefix,"\n")
cat("Kmer type:", kmertype,"\n")
cat("Software file:", software_file, "\n")
cat("Script location:", script_location, "\n")
cat("Rscript path:", Rscriptpath, "\n")
cat("Kmer length:", kmerlen,"\n")
cat("Min count:", mincount,"\n")
cat("#############################################", "\n\n")



if(t>p) {
cat("Task",t,"not required\n")
quit("no")
}

# Create pattern batch
# Function returns total number of kmers
batch = create_pattern_batch(fullkmerlistfile = fullkmerlistfile, p = p, t = t, stringlist2patternpath = stringlist2patternpath, kmerlist2patternpath = kmerlist2patternpath, kmercountslistfile = kmercountslistfile, kmerlen = kmerlen, mincount = mincount, kmertype = kmertype, output_prefix = output_prefix, output_dir = output_dir) # SGE changed to function

cat("\n")

# Create kinship matrix for the batch
create_kinship_batch(batch_prefix = batch$output_prefix_batch, batches_dir = batch$batches_dir, nkmersbatch = batch$nkmersbatch[t,], patterncountspath = patterncountspath, pattern2kinshippath = pattern2kinshippath, kmertype = kmertype, kmerlen = kmerlen)

cat("\n")


# For a subset of the processes, use to merge the batches created by the other processes
# SGE incorporated from patternmerge.Rscript

# Which processes to keep
p = floor(p/5); p = ifelse(p==0, 1, p)

if(t>p) {
	
	cat("Created pattern batch",t,"\n")
	cat("Task",t,"not required for pattern and kinship matrix merging","\n")
	cat("Finished in",(proc.time()[3]-start.time)/3600,"hours\n")
	quit("no")
	
} else {
	
	cat("Created pattern and kinship batch",t,"\n")
	cat("Running kmer pattern merge","\n")
	# Get merge parameters
	merge_parameters = merge_pattern_batches_parameters(nkmersbatch = batch$nkmersbatch, p = p, t = t, output_dir = output_dir, output_prefix = output_prefix, kmertype = kmertype, kmerlen = kmerlen)
	
	# Merge patterns
	merge_patterns(t = t, n = merge_parameters$n, b = merge_parameters$b, p = p, imax = merge_parameters$imax, prefix = output_prefix, files = merge_parameters$infiles, patternmergepath = patternmergepath, output_dir = output_dir, kmertype = kmertype, kmerlen = kmerlen, batches_dir = batch$batches_dir)
	cat("\n")
	cat("Running kinship matrix merge","\n")
	# Merge kinship matrices
	kinship_matrix_files = gsub(".patternKey.txt.gz", ".kinship.txt.gz", merge_parameters$infiles, fixed = TRUE)
	merge_kinship_matrices(t = t, n = merge_parameters$n, b = merge_parameters$b, p = p, imax = merge_parameters$imax, files = kinship_matrix_files, prefix = output_prefix, nkmersbatch = batch$nkmersbatch, output_dir = output_dir, kmertype = kmertype, kmerlen = kmerlen, batches_dir = batch$batches_dir)
	cat("\n")

	
}

# Keep the last process to calculate the MACs for the patterns
if(t==p){
	
	cat("Calculating number of genomes each pattern is present in (of the genomes with a non NA phenotype)","\n")
	stopifnot(system(paste0(Rscriptpath, " ", pattern2presencecountRscript, " ", output_dir, output_prefix, "_", kmertype, kmerlen, " ", output_dir, " ", id_file, " FALSE"))==0)
	cat("\n")
	
}

cat("Finished in",(proc.time()[3]-start.time)/3600,"hours\n")