########################################################################################################
# This is software for SDWT, Stationary Discrete Wavelet Transformation,by Jing Wei 2276887838@qq.com, Jie Lin at linjie891@163.com,or Yue Jiang at yueljiang@163.com.
#
# SDWT implements an alignment-free approach which is a k-mer based frequency profile feature method.
# There are 4 steps in SDWT. 
# (1)kmer Extraction: it breaks a sequence (DNA or protein sequence) into k-mers.
# (2) NumberTransform: it turn a k-mer into a complex number, then, a sequence can be represented as a series of complex numbers
# (3) FeatureTransform: a series of complex numbers is transformed by using stationary discrete wavelet transformation equation Haar into a feature vector.
# (4) Application: The feature vector is used to classification and/or cluster 
######################################################################################################
#Examples for usage:

#Classification example:
#> source("SDWTfunction.R")  #execute SDWT scripts which is in a file named SDWTfunction.R
#> seqsAndLabel = readFile("traindata/file.lst")  # read in file 
#> testseqs = testReadFile("testdata/file.lst")
#> traindata = TrainFeature(3,seqsAndLabel$seqs,3,4,DNAalphabet)
#> testfeature = TestFeature(3,testseqs,3,traindata$means,traindata$SD,4,DNAalphabet)
#> library(class)
#> predictResult = knn(traindata$trainseqfeature,testfeature,seqsAndLabel$seqLabel,k=1)
#
#
#-------------------------------------------------------------------
#
#Clustering example:
#> source("SDWTfunction.R")
#> seqsAndLabel = readFile("cluster/file.lst")
#> feature = Feature(3,seqsAndLabel$seqs,3,4,DNAalphabet)
#> KmeansResult = kmeans(feature,100)
######################################################################################################


DNAalphabet="ACGT"  # DNA alphabet 
proteinalphabet="ACDEFGHIKLMNPQRSTVWY"   # protein alphabet 


#########################################################
#function: Connect，Read in a family and stitching.
#   input: 
#          x: a DNA/protein family.
#          n1: The number of sequences that a family contains.
#  output:  
#          y: Splicing the complete sequence in the current family.
#   return:	the same as output		
#   Usage: oneSeq_connet = Connect(x,n1).
#    Date: March 27, 2018
#######################################################
Connect=function(x,n1)
{
	y=array()
	m=0
	i=1
	for(j in 1:n1)
	{
		temp_x=NULL
		if(substr(x[i],1,1)==">")
		{#If the first character read is ">", i adds 1 and reads the next line
			m=m+1
			i=i+1
		}	
		while(m==j)
		{#When m=j, control the spliced ​​newline
			if(substr(x[i],1,1)!=">"&i!=(length(x)+1))
			{#Splicing if the read first character is not ">"
				temp_x=paste(temp_x,x[i],sep="")
				i=i+1
			}
			else break#Otherwise jump out of the loop without stitching, start the next stitching
		}
		y[j]=temp_x	
	}
	return(y)
}

#########################################################
#function: readFile，Read in all the family and stitching.
#   input: 
#			filename: the input file name, which contains a list of file names. Each file contains the corresponding (DNA/protein) sequences in FST(DNA) or FASTA(protein).
#  output(return):  
#           seqsAndLabel: a matrix named seqsAndLabel is generated. Each line of matrix contains two elements:  sequence and its corresponding family label.			
#   Usage: seqsAndLabel = readFile("familylist").
#    Date: March 27, 2018
#######################################################
readFile = function(filename)
{
	read_file = read.table(filename,header = FALSE)
	family_nums = nrow(read_file)
	seqs = NULL
	seqLabel = NULL
	label = 0
	path = strsplit(filename,split="/")[[1]][1]
	for(i in 1:family_nums)
	{
		temp_a=readLines(paste(path,"/",read_file[i,1],sep=""),encoding = "UTF-8")#read all files in a specific path one by one which is indicated in the file.
		n=0
		label = label+1
		for(j in 1:length(temp_a))
		{
			if(substr(temp_a[j],1,1)==">")
			n=n+1
		}
		temp_label = rep(label,n)
		seqLabel = c(seqLabel,temp_label)
		oneSeq_connet = NULL
		oneSeq_connet = Connect(temp_a,n)#Splicing the sequence in the current family	
		seqs=c(seqs,oneSeq_connet)	#Combine the sequence from each family read in one file
	}
	outputSeqLabel = list(seqs=seqs,seqLabel=seqLabel)
	return (outputSeqLabel)
}

#########################################################
#function: testreadFile，Read all test set family files and stitch
#   input: 
#			filename: the input test file name, which contains a list of test file names. Each file contains the corresponding (DNA/protein) sequences in FST(DNA) or FASTA(protein).
#  output(return):  
#			allHOGseqs: test data sequence
#   Usage: testseqs = testReadFile("familylist")
#    Date: March 27, 2018
#######################################################
testReadFile = function(filename)
{
	read_file = read.table(filename,header = FALSE)
	family_nums = nrow(read_file)
	allHOGseqs = NULL
	path = strsplit(filename,split="/")[[1]][1]
	for(i in 1:family_nums)
	{
		temp_a=readLines(paste(path,"/",read_file[i,1],sep=""),encoding = "UTF-8")#read all files in a specific path one by one which is indicated in the file.
		n=0
		for(j in 1:length(temp_a))
		{
			if(substr(temp_a[j],1,1)==">")
			n=n+1
		}
		oneSeq_connet = NULL
		oneSeq_connet = Connect(temp_a,n)#Splicing the sequence in the current family
		allHOGseqs=c(allHOGseqs,oneSeq_connet)	#Combine the sequence from each family read in one file
	}
	return (allHOGseqs)
}

#########################################################
#function: kmer，calculate  k-mers frequencies in the input sequences
#   input: 
#			seq: input DNA/protein sequences
#			K: the length of k in integer number
#			alphabet: the actual alphabet for DNA/protein. For DNA, it is the set "ACGT". For protein, it is the set "ACDEFGHIKLMNPQRSTVWY"
#			length(dat): The number of input sequences
#  output(return):  
#			frequency of each k-mer 
#   Usage:  testseqs = testReadFile("familylist")
#    Date:  March 27, 2018
#######################################################
kmer<-function(seq, k, alphabet, len=length(seq))
{
    dyn.load("Rkmer.so")
	alphasize=nchar(alphabet)
	matrix(.C("Rkmer",seq,as.integer(k),as.integer(len),alphabet,retArray=integer(len*alphasize^k))$retArray,nrow=len,ncol=alphasize^k,byrow=T)
}

#########################################################
#function: TrainFeature，Get the feature vector of the training set
#   input: 
#			K:  the length of k in integer number. 
#   		dat:  DNA/protein sequences .
#    		decomLevel: decomposition level of stationary wavelet transform in integer.  
#			char: the alphabet size where DNA is 4 and protein is 20.
#			alphabet: the actual alphabet for DNA/protein. For DNA, it is the set "ACGT". For protein, it is the set "ACDEFGHIKLMNPQRSTVWY".
#	output(return):  
#			trainseqfeature: a matrix contains feature vector.
#      		means: mean of a specific k-mer.
#           SD: Standard Deviation of a specific k-mer.
#   Usage:  trainfeature = TrainFeature(3,seqsAndLabel$seqs,3,4,DNAalphabet)
#    Date:  March 27, 2018
#######################################################
TrainFeature =function(K,dat,decomLevel,char,alphabet)
{
	angle=360/char^K
	K_word=array(0:0,c(char^K,2))
	angle_sum = 0
	for(i in 1:char^K)
	{
		angle_sum = angle_sum + angle
		K_word[i,1]=sinpi(angle_sum /180)
		K_word[i,2]=cospi(angle_sum/180)	
	}
	##########
	wordFrequence = kmer(dat,K,alphabet,length(dat))
	StandFrequency = scale(wordFrequence)
	means = attr(StandFrequency ,"scaled:center")
	SD = attr(StandFrequency ,"scaled:scale")
	fre_real=array(0:0,c(nrow(StandFrequency),char^K))
	fre_imag=array(0:0,c(nrow(StandFrequency),char^K))
	for (i in 1:nrow(StandFrequency)) 
	{  
		#Word Frequency Multiplied by Complex Numbers
		for(j in 1:char^K)
		{
			fre_real[i,j]=StandFrequency[i,j]*(K_word[j,1])#Real
			fre_imag[i,j]=StandFrequency[i,j]*(K_word[j,2])#Imaginary
		}
	}
	complexFileName = cbind(fre_real,fre_imag)#Real and imaginary mergers
	if (nrow(complexFileName)%%(2^decomLevel)!=0)
	{
		add = array(0:0,c((2^decomLevel)-nrow(complexFileName)%%(2^decomLevel),ncol(complexFileName)))
		complexFileName = rbind(complexFileName,add)
	}
	write.table(complexFileName,file="traincomplexFileName.csv",append=TRUE,row.names = FALSE, col.names=FALSE, sep=",")
	mat=paste("matlab -nojvm -nodisplay -nosplash -nodesktop -r \"trainSDWTgetFeature(\'traincomplexFileName.csv\'",decomLevel,K,sep=",")
	mat=paste(mat,")\"",sep="")
	system(mat)
	trainseqfeature = read.csv("trainfeatureFileName.csv",header=FALSE)[1:length(dat),]
	output = list(trainseqfeature=trainseqfeature,means=means,SD=SD)
	return (output)
}

######################################################################################
#function: TestFeature，Get the feature vector of the test set
#   input: 
#			K:  the length of k in integer number. 
#		    dat:  DNA/protein sequences.
#  			decomLevel: decomposition level of stationary wavelet transform in integer.  
#			char: the alphabet size where DNA is 4 and protein is 20. 
#			alphabet: the actual alphabet for DNA/protein. For DNA, it is the set "ACGT". For protein, it is the set "ACDEFGHIKLMNPQRSTVWY".
#	output(return):  
#			testfeature: a matrix contains feature vector. 
#   Usage:  testfeature = TestFeature(3,testseqs,3,traindata$means,traindata$SD,4,DNAalphabet)
#    Date:  March 27, 2018
####################################################################################
TestFeature = function(K,dat,decomLevel,center,sde,char,alphabet)
{
	angle=360/char^K
	K_word=array(0:0,c(char^K,2))
	angle_sum = 0
	for(i in 1:char^K)
	{
		angle_sum = angle_sum + angle
		K_word[i,1]=sinpi(angle_sum /180)
		K_word[i,2]=cospi(angle_sum/180)	
	}
	##########
	wordFrequence = kmer(dat,K,alphabet,length(dat))
	StandFrequency=array(0:0,c(nrow(wordFrequence),char^K))
	for (i in 1:nrow(wordFrequence))
	{
		for(j in 1:char^K)
		{
			StandFrequency[i,j]=(wordFrequence[i,j]-center[j])/sde[j]
		}
	}
	fre_real=array(0:0,c(nrow(StandFrequency),char^K))
	fre_imag=array(0:0,c(nrow(StandFrequency),char^K))
	for (i in 1:nrow(StandFrequency)) 
	{  
		for(j in 1:char^K)
		{
			fre_real[i,j]=StandFrequency[i,j]*(K_word[j,1])#
			fre_imag[i,j]=StandFrequency[i,j]*(K_word[j,2])#
		}
	}
	complexFileName = cbind(fre_real,fre_imag)#
	if (nrow(complexFileName)%%(2^decomLevel)!=0)
	{
		add = array(0:0,c((2^decomLevel)-nrow(complexFileName)%%(2^decomLevel),ncol(complexFileName)))
		complexFileName = rbind(complexFileName,add)
	}
	write.table(complexFileName,file="testcomplexFileName.csv",append=TRUE,row.names = FALSE, col.names=FALSE, sep=",")
	mat=paste("matlab -nojvm -nodisplay -nosplash -nodesktop -r \"testSDWTgetFeature(\'testcomplexFileName.csv\'",decomLevel,K,sep=",")
	mat=paste(mat,")\"",sep="")
	system(mat)
	testfeature = read.csv("testfeatureFileName.csv",header=FALSE)[1:length(dat),]
	return (testfeature)
}

#########################################################
#function: Feature，Get the feature vector of the input sequence characteristics
#   input: 
#			K:  the length of k in integer number. 
#		    dat:  DNA/protein sequences.
#  			decomLevel: decomposition level of stationary wavelet transform in integer.  
#			char: the alphabet size where DNA is 4 and protein is 20. 
#			alphabet: the actual alphabet for DNA/protein. For DNA, it is the set "ACGT". For protein, it is the set "ACDEFGHIKLMNPQRSTVWY".
#	output(return):  
#			feature: a matrix contains feature vector. 
#   Usage:  feature = Feature(3,testseqs,3,traindata$means,traindata$SD,4,DNAalphabet).
#    Date:  March 27, 2018
#######################################################
Feature = function(K,dat,decomLevel,char,alphabet)
{
	angle=360/char^K
	K_word=array(0:0,c(char^K,2))
	angle_sum = 0
	for(i in 1:char^K)
	{
		angle_sum = angle_sum + angle
		K_word[i,1]=sinpi(angle_sum /180)
		K_word[i,2]=cospi(angle_sum/180)	
	}
	##########
	wordFrequence = kmer(dat,K,alphabet,length(dat))
	StandFrequency = scale(wordFrequence)
	
	fre_real=array(0:0,c(nrow(StandFrequency),char^K))
	fre_imag=array(0:0,c(nrow(StandFrequency),char^K))
	for (i in 1:nrow(StandFrequency)) 
	{ 
		for(j in 1:char^K)
		{
			fre_real[i,j]=StandFrequency[i,j]*(K_word[j,1])#
			fre_imag[i,j]=StandFrequency[i,j]*(K_word[j,2])#
		}
	}
	complexFileName = cbind(fre_real,fre_imag)#
	if (nrow(complexFileName)%%(2^decomLevel)!=0)
	{
		add = array(0:0,c((2^decomLevel)-nrow(complexFileName)%%(2^decomLevel),ncol(complexFileName)))
		complexFileName = rbind(complexFileName,add)
	}
	write.table(complexFileName,file="complexFileName.csv",append=TRUE,row.names = FALSE, col.names=FALSE, sep=",")
	mat=paste("matlab -nojvm -nodisplay -nosplash -nodesktop -r \"SDWTgetFeature(\'complexFileName.csv\'",decomLevel,K,char,sep=",")
	mat=paste(mat,")\"",sep="")
	system(mat)
	feature = read.csv("featureFileName.csv",header=FALSE)[1:length(dat),]
	return (feature)
}


