plANI_v2.Rfunction

#' plani
#'
#' Calculates PLANI (PLasmid Average Nucleotide Identity) for metagenomic contigs.
#' @param contigs is the full path to a fasta file containing metagenomic contigs.
#' @param plasmids is the full path to a fasta file containing the reference plamsids used as database.
#' @param outdir a name for the output directory that will be created for results.
#' @param win is the window size for ANI fragments.
#' @param winst is the step for moving the window.
#' @param amin is the minimum size (percentage of fragment length) for any blast hit to be considered.
#' @param fmin is the minimum fraction of hits from a contig to match a certain plasmid and to be reported.
#' @param proc is the number of cores.
#' @keywords ANI plasmid blastn
#' @export

plani <- function(

			contigs,
			plasmids,
			
			win      = 1000,
			winst    = 200,
			amin     = 0.7,
			fmin     = 0.5,
			proc     = 2,
			
			outdir   = 'plani_result'
			
)	{

	# Options #
	
	options(getClass.msg=FALSE)
	os<-.getOS()
		
	# Select OS #
	
	if (os=='linux'){
	
		blastn <- paste(system.file('blast',package='taxxo'),'/linux/blastn',sep='')
		mkbldb <- paste(system.file('blast',package='taxxo'),'/linux/makeblastdb',sep='')
		bldbcd <- paste(system.file('blast',package='taxxo'),'/linux/blastdbcmd',sep='')
		
	} else if (os=='darwin'){
		
		blastn <- paste(system.file('blast',package='taxxo'),'/darwin/blastn',sep='')
		mkbldb <- paste(system.file('blast',package='taxxo'),'/darwin/makeblastdb',sep='')
		bldbcd <- paste(system.file('blast',package='taxxo'),'/darwin/blastdbcmd',sep='')

	} else {
		
		stop('ERROR: Unknown OS, unable to proceed.')
	}
	
	# Dependencies #

	suppressMessages(library(foreach,quietly=T))
	suppressMessages(library(doMC,quietly=T))
	suppressMessages(library(seqinr,quietly=T))
	suppressMessages(library(Biostrings,quietly=T))
		
	# Define internal function #

	# chunker

	chunker <- function(m,n) {

		s <- split(m,cut(seq_along(m),n,labels=F))
		
		return(s)
	}
	
	# spliter
	
	spliter <- function(
	
				w=win,
				s=winst,
				v
				
	) {
	
		n  <- length(v)
		x  <- 1
		x2 <- w
		o  <- 1
		l  <- list()

		if ( n>w & x2<=(n-s) ) {

			while ( x2<=(n-s) ) {
	
				ini    <- x
				fin    <- x2
				l[[o]] <- v[ini:fin]
				x      <- x+s
				x2     <- x2+s
				o      <- o+1

			}

			l[[o]] <- v[(fin+1):n]
			
		} else {
						
			l[[1]]<-v
		}
	
		return(l)
	}


	# plANI parallel #
	
	planiparallel <- function(y)
	
	
	{
		
		setwd(outdir)
	
		chun<-chunks[[y]]

			
		# Generate chunk fasta #
		
		chu <- as.vector(chun)
		fc  <- paste('chunk',y,'.txt',sep='')
		fc2 <- gsub('.txt','.fasta',fc)

		cat(chu,sep='\n',file=fc)
			
		cmd <- paste(
						
					bldbcd,
					' -entry_batch ',fc,
					' -db ./contigs_database/contigsdb > ',
					fc2,sep=''
					
					)
						
		system(cmd,ignore.stderr=T)
		
		fc3 <- read.fasta(fc2)
		fs3 <- lapply(getSequence(fc3),toupper)
		fn3 <- getName(fc3)
								
		# Create output 
			
		rmat           <- matrix(nrow=length(chun),ncol=length(nw),0)
		colnames(rmat) <- nw
		rownames(rmat) <- fn3

		for ( i in 1:length(fn3) ) {

			# BLAST
		
			blout1 <- paste('blout1',y,i,sep='.')
			query  <- paste('query1',y,i,'fasta',sep='.')
			q1     <- fs3[[i]]
				

			query1  <- spliter(v=q1,w=win,s=winst)
			headers <- paste(fn3[i],'_fragment_',seq(1:length(query1)),sep='')
					
			fragmin <- round(length(query1)*fmin,digits=0)

			write.fasta(query1,names=headers,file=query)
					
					
			blcmd<-paste(
				
						blastn,
						" -query ",query,
						" -db ./plasmids_database/plasmidsdb ",
						" -xdrop_gap_final 150 ",
						" -outfmt '6 qseqid sseqid pident length' -out ",
						blout1,
						sep=''
							
						)
						 
			system(blcmd,ignore.stderr=T,ignore.stdout=T)
			
			if ( file.info(blout1)$size>0 ) {
				
				tab1           <- read.csv(blout1,sep='\t',header=F)
				colnames(tab1) <- c('qid','sid','pid','len')
				
				amini <- round(win*amin,digits=0)
					
				tab2 <- tab1[which(tab1$len>=amini),]
								
				if ( dim(tab2)[1]>0 ) {
						
					tb <- table(tab2$sid)
					tf <- which(tb>=fmin)

					if ( length(tf)>0 ) {
							
						pla <- names(tb)[tf]
							
						for (p in pla) {
								
							pids <- tab2[which(tab2$sid==p),'pid']
							pmn <- mean(as.numeric(as.vector(pids)))
							
							rmat[i,which(nw==p)] <- pmn
						}
					}
				}
			}
			
			system(paste('rm -rf chunk*.txt chunk*.fasta',blout1,query))
		}
		
		return(rmat)
	}
	
	

	# Make output directory #
	
	system(paste('mkdir',outdir))


	# Parse plasmid file #
	
	system(paste('mkdir ',outdir,'/plasmids_database',sep=''))
		
	system(
			paste(	
				'cp ',
				plasmids,
				' ',outdir,
				'/plasmids_database/plasmids.fasta',
				sep=''
				)
			)
	
	setwd(outdir)
	setwd('plasmids_database')
	
	pl <- readDNAStringSet('plasmids.fasta')
	nm <- gsub('>','',system('grep ">" plasmids.fasta',intern=T))
	nw <- paste('plasmid',seq(1,length(nm)),sep='_')
	
	names(pl) <- nw
	
	writeXStringSet(pl,file='new_plasmids.fasta',format='fasta')
	
	rm(pl)
	system('rm -rf plasmids.fasta')
	
	mkdb <- paste(
	
				  mkbldb,
				  ' -in new_plasmids.fasta',
				  ' -dbtype nucl',
				  ' -hash_index -parse_seqids',
				  ' -title plasmidsdb',
				  ' -out plasmidsdb',
				  sep=''
				  
				  )
				  
	system(mkdb,ignore.stderr=T,ignore.stdout=T)
	
	pldf           <- as.data.frame(cbind(nw,nm))
	colnames(pldf) <- c('Identifier','Original')

	write.table(pldf,
				sep='\t',
				quote=F,
				row.names=F,
				file='../plasmids_names.csv'
				
				)
	
	setwd('../../')
	

	# Parse contigs #
	
	system(paste('mkdir ',outdir,'/contigs_database',sep=''))

	
	system(
			paste(	
				'cp ',
				contigs,
				' ',outdir,
				'/contigs_database/contigs.fasta',
				sep=''
				)
			)

	setwd(outdir)
	setwd('contigs_database')
	
	ct <- readDNAStringSet('contigs.fasta')
	nc <- gsub('>','',system('grep ">" contigs.fasta',intern=T))
	wc <- paste('contig',seq(1,length(nc)),sep='_')
	
	names(ct) <- wc
	
	writeXStringSet(ct,file='new_contigs.fasta',format='fasta')

	rm(ct)
	system('rm -rf contigs.fasta')

	mkdb2 <- paste(
	
			  mkbldb,
			  ' -in new_contigs.fasta',
			  ' -dbtype nucl',
			  ' -hash_index -parse_seqids',
			  ' -title contigsdb',
			  ' -out contigsdb',
			  sep=''
				  
			  )
				  
	system(mkdb2,ignore.stderr=T,ignore.stdout=T)


	ctdf           <- as.data.frame(cbind(wc,nc))
	colnames(ctdf) <- c('Identifier','Original')

	write.table(ctdf,
				sep='\t',
				quote=F,
				row.names=F,
				file='../contig_names.csv'
				
				)

	setwd('../../')

	
	# Check cores #
	
	lc <- length(wc)
	
	if (proc > lc) {
		
		proc <- lc
	}

	
	# Run #
		
	chunks<-chunker(wc,proc)
		
	registerDoMC(proc)

	final<-foreach(y=1:proc) %dopar% {
		
		planiparallel(y)	
	}
		
	plani_result <- do.call('rbind',final)
	
	setwd(outdir)
	
	saveRDS(plani_result,file='plani_result.rds')
	
	setwd('../')
}