#include <cstdio>
#include <cstdlib>
#include <iostream>
#include <string>
#include <fstream>
#include <map>
#include <sys/stat.h>



using namespace std;

//subroutine for converting sequence to base 4 and the output to base 10 representation
//returns base 4 value converted to base 10

int strval(string str){
	int f=0;
	int sum=0;
	string seq (str);
   	for (int i=0; i<seq.length(); i++){
        	int power=1;
        	if (seq[i]==('a')){
        	        f=1;
        	        for (int j=0; j<i; j++){
        	        	power=power*4;
        	        }
        	        sum=sum+power*f;
        	        power=1;
        	}
        	if (seq[i]==('t')){
        	        f=2;
        	        for (int j=0; j<i; j++){
        	        	power=power*4;
        	        }
        	        sum=sum+power*f;
        	        power=1;
        	}
        	if (seq[i]==('c')){
        	        f=3;
        	        for (int j=0; j<i; j++){
        	        	power=power*4;
        	        }
        	        sum=sum+power*f;
        	        power=1;               
        	}
    	}
    	return sum;
}


//subroutine for checking whether middle input fragment has only 1 mutation
//returns mutation position relative to fragment

int baseconverter(int number){
	int power=1;
	int n;
	for(int i=0;;i++){
		if(number-power<0){
        		n=i-1;
        		break;
		}
		power=power*4;
	}
	int array4[n];

	for(int kk=n;kk>=0;kk--){
		power=power/4;
        	for(int j=0;j<4;j++){
                	if(number-(j+1)*power<0){
                        	array4[kk]=j;
                        	number=number-j*power;
                        	break;
                	}
        	}
	}

	int zz=0;
	int tag=0;
	for(int i=0;i<=n;i++){
		if(array4[i]!=0){
        		tag=i;
        		zz++;
		}
	}

	if(zz==1){
		return tag;
	}
	if(zz!=1){
		return -1;
	}
}


//This subroutine returns the base 4 value of the genomic region corresponding to the third fragment of an input sequence (8 nts for k=13)

int baseconverter2(int number, int end){
	int power=1;
	int n;
	for(int i=0;;i++){
		if(number-power<0){
        		n=i-1;
        		break;
		}
		power=power*4;
	}
	
	int array4[n];

	for(int kk=n;kk>=0;kk--){
		power=power/4;
        	for(int j=0;j<4;j++){
                	if(number-(j+1)*power<0){
                        	array4[kk]=j;
                        	number=number-j*power;
                        	break;
                	}
        	}
	}

	
	int revisedsum=0;
	int revisedpower=1;
	for(int i=0;i<end;i++){
		revisedsum=revisedsum+array4[i]*revisedpower;
		revisedpower=revisedpower*4;
	}

	return revisedsum;
}


//subroutine for converting base 4 sequence representation to nucleotides
//returns 13-mer sequences for k=13

string baseconverter3(int number){
	int power=1;
	int n;
	string wtseq="";
	string base;
	for(int i=0;;i++){
		if(number-power<0){
			n=i-1;
			break;
		}
		power=power*4;
	}

	int array4[n];

	for(int kk=n;kk>=0;kk--){
		power=power/4;
        	for(int j=0;j<4;j++){
        	        if(number-(j+1)*power<0){
                        	array4[kk]=j;
                        	number=number-j*power;
                        	break;
                	}
        	}
	}

	for(int i=0;i<=n;i++){
		if(array4[i]==0){
			base="g";
		}
		if(array4[i]==1){
			base="a";
		}
		if(array4[i]==2){
			base="t";
		}
		if(array4[i]==3){
			base="c";
		}

		wtseq.append(base);
	}

	int check=13-wtseq.length();
	if (check!=0){
		for (int i=0;i<check;i++){
			wtseq.append("g");
		}
	}
	
	return wtseq;
}


//begin main program here

int main(int nNumberofArgs, char* pszArgs[]){

	system("clear");

//	char buffer[1000];
//	std::cout << "Current directory is: " << getcwd(buffer, 1000) << "\n";

	cout<<"			'galign', a SNP detecting application"<<endl;
	cout<<"			 ------------------------------------"<<endl<<endl;
	cout<<"This program aligns short reads (27-39 bps) obtained from genome-wide"<<endl;
	cout<<"sequencing to a genome of interest."<<endl<<endl;
	cout<<"You will need to input two file names.  The first is the genome sequence"<<endl;
	cout<<"file, which can be found in the 'Genome_sequences' folder.  The second"<<endl;
	cout<<"file is the sequence reads file, which can be found in the 'Sequence_reads'"<<endl;
	cout<<"folder."<<endl<<endl;
	cout<<"The output of this program will be placed in a file in the folder"<<endl;
	cout<<"'Alignment_results' as 'yourname_galign.txt'"<<endl<<endl;

//input relevant parameters
	string gseq;
	string gseq1;

	int r=0;

	while(r==0){
		gseq1="Desktop/galign/Genome_sequences/";
		cout<<"Genome sequence file: ";
		cin>>gseq;
		gseq1.append(gseq);
	
		FILE *fp;
		fp=fopen(gseq1.c_str(),"r");
		if(fp!=NULL){
			r=1;
		}
	}

	string seqread;
	string seqread1;

	int rr=0;
	while(rr==0){
		seqread1="Desktop/galign/Sequence_reads/";
		cout<<"Sequence reads file: ";
		cin>>seqread;
		seqread1.append(seqread);
	
		FILE *fpp;
		fpp=fopen(seqread1.c_str(),"r");
		if(fpp!=NULL){
			rr=1;
		}
	}

	string yourname;
	string yourname1="Desktop/galign/Alignment_results/";
	cout<<"Your first name: ";
	cin>>yourname;
	yourname1.append(yourname);
	yourname1.append("_galign.txt");

	int k=13;

	int arraysize=1;
	for (int i=1; i<=k;i++){
		arraysize=4*arraysize;
	}

	int *array;
	array=new int[arraysize];
	for (int i=0;i<arraysize;i++){
		array[i]=0;
	}

	string line;
	int length=0;

	ifstream handle (gseq1.c_str());
	while (! handle.eof()){
		getline (handle, line);
		length=length+line.length();
	}
	handle.close();

	int *genomearray;
	genomearray=new int[length];
	for (int i=0; i<length;i++){
		genomearray[i]=0;
	}

	string test1;
	string test;
	string bit;
	int t=0;
	int j=0;

	cout<<"Building genome index..."<<endl;

	ifstream handle1 (gseq1.c_str());
	while (! handle1.eof()){
		getline (handle1, test1);
		test.append(test1);
		if(test.length()>=k){
			for (int i=0; i<=test.length()-k; i++){
                bit=test.substr(i,k);
                t=strval(bit);
                genomearray[j]=t;
                array[t]=array[t]+1;
//                cout<<j<<"\t"<<genomearray[j]<<"\t"<<array[t]<<endl;
                j++;
			}
			test=test.substr(test.length()-k+1,k-1);
		}
	}
	handle1.close();

	int g=0;

	for (int i=0;i<arraysize;i++){
		if (array[i]==0){
			g++;
		}
	}

	for (int i=0;i<arraysize;i++){
		if (array[i]==0){
			array[i]=1;
		}
	}

	int pointerarraysize=length+g;

	int *array2;
	array2=new int[pointerarraysize];

	for (int i=0;i<pointerarraysize;i++){
		array2[i]='\0';
	}

	int gg=0;

	int* *pArray;
	pArray=new int*[arraysize+1];
	pArray[0]=&array2[0];


	for(int i=1;i<arraysize+1;i++){
		gg=gg+array[i-1];
		pArray[i]=&array2[gg];
	}

	delete[] array;

	int *array3;
	array3=new int[pointerarraysize];
	for (int i=0;i<pointerarraysize;i++){
		array3[i]=0;
	}


	cout<<"Assigning pointers..."<<endl;

	string testt;
	j=0;
	ifstream handle2 (gseq1.c_str());
	while (! handle2.eof()){
		getline (handle2, test1);
		testt.append(test1);
		if(testt.length()>=k){
			for (int i=0; i<=testt.length()-k; i++){
               bit=testt.substr(i,k);
               t=strval(bit);
               *(pArray[t]+array3[t])=j+1;
               array3[t]=array3[t]+1;
               j++;
			}
			testt=testt.substr(testt.length()-k+1,k-1);
		}
	}
	handle2.close();

	delete[] array3; //deleting array3


//Now input test read to see if can find it

	cout<<"Performing alignments..."<<endl;

	int *wtarray;
	wtarray=new int[length];

	int *mutarray;
	mutarray=new int[length];

	string *sequencearray;
	sequencearray=new string[length];

	int *mutationposition;
	mutationposition=new int[length];

	for (int i=0;i<length;i++){
		wtarray[i]=0;
		mutarray[i]=0;
		sequencearray[i]="";
		mutationposition[i]=-1;
	}

	ifstream myfile1 (seqread1.c_str());
	string input1;
	getline (myfile1,input1);
	int esize=input1.length()-2*k;
	myfile1.close();


	cout<<endl;
	
	ifstream myfile (seqread1.c_str());
	string input;
	int q=0;
	while (! myfile.eof()){
		getline (myfile, input);
		q++;
		cout<<"\033[1A";
		cout<<q<<endl;

		string frag;
		int nt, nt1;
		int dif, dif1, dif2, dif3;
		int ind1, ind2;
		int mutationsite;
		int bcreturn;
		int difference;
		int positive=0;

		frag=input.substr(0,k);
		nt=strval(frag);
		dif=pArray[nt+1]-pArray[nt];
		if (dif>15){
			continue;
		}


		for (int j=1;j<=dif;j++){
			ind1=*(pArray[nt]+j-1);
			if (ind1>length-4*k){
				break;
			}
			difference=strval(input.substr(2*k,esize))-baseconverter2(genomearray[ind1+2*k-1],esize);
			if(difference<0){
				difference=-difference;
			}
			if(genomearray[ind1+k-1]==strval(input.substr(k,k))){
				if (difference==0){
					for(int pp=ind1;pp<=ind1+input.length()-1;pp++){
						wtarray[pp]=wtarray[pp]+1;
					}
					positive=1;
					break;
				}else{
					bcreturn=baseconverter(difference);
					if (bcreturn>=0){
						mutationsite=ind1+2*k+bcreturn;
//		cout<<ind1<<"\t"<<ind1+input.length()-1<<"\t"<<mutationsite<<endl;
						for(int pp=ind1;pp<mutationsite;pp++){
							wtarray[pp]++;
						}
						for(int pp=mutationsite+1;pp<=ind1+input.length()-1;pp++){
							wtarray[pp]++;
						}
						mutarray[mutationsite]++;
						sequencearray[mutationsite]=input;
						mutationposition[mutationsite]=mutationsite-ind1;
						positive=1;
						break;
					}
				}
			}
			
			if(difference==0){
				dif1=genomearray[ind1+k-1]-strval(input.substr(k,k));
				if(dif1<0){
					dif1=-dif1;
				}
				bcreturn=baseconverter(dif1);
				if (bcreturn>=0){
					mutationsite=ind1+k+bcreturn;
//			cout<<ind1<<"\t"<<ind1+input.length()-1<<"\t"<<mutationsite<<endl;
					for(int pp=ind1;pp<mutationsite;pp++){
						wtarray[pp]++;
					}
					for(int pp=mutationsite+1;pp<=ind1+input.length()-1;pp++){
						wtarray[pp]++;
					}
					mutarray[mutationsite]++;
					sequencearray[mutationsite]=input;
					mutationposition[mutationsite]=mutationsite-ind1;
					positive=1;
					break;
				}
			}
			frag=input.substr(k,k);
			nt1=strval(frag);
			dif2=pArray[nt1+1]-pArray[nt1];
			for (int jj=1;jj<=dif2;jj++){
				ind2=*(pArray[nt1]+jj-1);
				if (ind2<k+1){
					continue;
				}
				difference=strval(input.substr(2*k,esize))-baseconverter2(genomearray[ind2+k-1],esize);
				if(difference<0){
					difference=-difference;
				}
		
				if(difference==0){
					dif3=genomearray[ind2-k-1]-strval(input.substr(0,k));
					if(dif3<0){
						dif3=-dif3;
					}
					bcreturn=baseconverter(dif3);
					if (bcreturn>=0){
						mutationsite=ind2-k+bcreturn;
//				cout<<ind2-k<<"\t"<<ind2+k-1+esize<<"\t"<<mutationsite<<endl;
						for(int pp=ind2-k;pp<mutationsite;pp++){
							wtarray[pp]++;
						}
						for(int pp=mutationsite+1;pp<=ind2+k-1+esize;pp++){
							wtarray[pp]++;
						}
						mutarray[mutationsite]++;
						sequencearray[mutationsite]=input;
						mutationposition[mutationsite]=mutationsite-ind2+k;
						positive=1;
						break;
					}
				}
			}
			if(positive==1){
				break;
			}
		}



//if alignment failed try using reverse complement of sequence read

		string reverse;

		if(positive==0){
			for(int i=input.length()-1;i>=0;i--){
				if (input.substr(i,1)=="g"){
					reverse=reverse.append("c");
					continue;
				}
				if(input.substr(i,1)=="a"){
					reverse=reverse.append("t");
					continue;
				}
				if(input.substr(i,1)=="t"){
					reverse=reverse.append("a");
					continue;
				}
				if(input.substr(i,1)=="c"){
					reverse=reverse.append("g");
					continue;
				}
				reverse=reverse.append("N");
			}

			frag=reverse.substr(0,k);
			nt=strval(frag);
			dif=pArray[nt+1]-pArray[nt];
			if (dif>15){
				continue;
			}

			for (int j=1;j<=dif;j++){
				ind1=*(pArray[nt]+j-1);
				if (ind1>length-4*k){
					break;
				}
				difference=strval(reverse.substr(2*k,esize))-baseconverter2(genomearray[ind1+2*k-1],esize);
				if(difference<0){
					difference=-difference;
				}
				if(genomearray[ind1+k-1]==strval(reverse.substr(k,k))){
					if (difference==0){
//			cout<<ind1<<"\t"<<ind1+reverse.length()-1<<endl;
						for(int pp=ind1;pp<=ind1+reverse.length()-1;pp++){
							wtarray[pp]++;
						}
						positive=1;
						break;
					}else{
						bcreturn=baseconverter(difference);
						if (bcreturn>=0){
							mutationsite=ind1+2*k-1+bcreturn;
//		cout<<ind1<<"\t"<<ind1+reverse.length()-1<<"\t"<<mutationsite<<endl;
							for(int pp=ind1;pp<mutationsite;pp++){
								wtarray[pp]++;
							}
							for(int pp=mutationsite+1;pp<=ind1+reverse.length()-1;pp++){
								wtarray[pp]++;
							}
							mutarray[mutationsite]++;
							sequencearray[mutationsite]=reverse;
							mutationposition[mutationsite]=mutationsite-ind1;
							positive=1;
							break;
						}
					}
				}
		
				if(difference==0){
					dif1=genomearray[ind1+k-1]-strval(reverse.substr(k,k));
					if(dif1<0){
						dif1=-dif1;
					}
					bcreturn=baseconverter(dif1);
					if (bcreturn>=0){
						mutationsite=ind1+k+bcreturn;
//			cout<<ind1<<"\t"<<ind1+reverse.length()-1<<"\t"<<mutationsite<<endl;
						for(int pp=ind1;pp<mutationsite;pp++){
							wtarray[pp]++;
						}
						for(int pp=mutationsite+1;pp<=ind1+reverse.length()-1;pp++){
							wtarray[pp]++;
						}
						mutarray[mutationsite]++;
						sequencearray[mutationsite]=reverse;
						mutationposition[mutationsite]=mutationsite-ind1;
						positive=1;
							break;
					}
				}
				
				frag=reverse.substr(k,k);
				nt1=strval(frag);
				dif2=pArray[nt1+1]-pArray[nt1];
				for (int jj=1;jj<=dif2;jj++){
					ind2=*(pArray[nt1]+jj-1);
					if (ind2<k+1){
						continue;
					}
					difference=strval(reverse.substr(2*k,esize))-baseconverter2(genomearray[ind2+k-1],esize);
					if(difference<0){
						difference=-difference;
					}
		
					if(difference==0){
						dif3=genomearray[ind2-k-1]-strval(reverse.substr(0,k));
						if(dif3<0){
							dif3=-dif3;
						}
						bcreturn=baseconverter(dif3);
						if (bcreturn>=0){
							mutationsite=ind2-k+bcreturn;
//				cout<<ind2-k<<"\t"<<ind2+k-1+esize<<"\t"<<mutationsite<<endl;
							for(int pp=ind2-k;pp<mutationsite;pp++){
								wtarray[pp]++;
							}
							for(int pp=mutationsite+1;pp<=ind2+k-1+esize;pp++){
								wtarray[pp]++;
							}
							mutarray[mutationsite]++;
							sequencearray[mutationsite]=reverse;
							mutationposition[mutationsite]=mutationsite-ind2+k;
							positive=1;
							break;
						}
					}
				}
				if(positive==1){
					break;
				}
			}
		}
	}

	myfile.close();


// print results

	ofstream file (yourname1.c_str());
	for (int i=0;i<length;i++){
		if(mutarray[i]==0){
			file<<wtarray[i]<<"\n"<<mutarray[i]<<endl;
		}
		if(mutarray[i]!=0 && mutationposition[i]<16){
			file<<wtarray[i]<<"\n"<<mutarray[i]<<"\n"<<sequencearray[i].substr(mutationposition[i],k)<<"\n"<<baseconverter3(genomearray[i-1])<<endl;
		}
		if(mutarray[i]!=0 && mutationposition[i]>=16){
			file<<wtarray[i]<<"\n"<<mutarray[i]<<"\n"<<sequencearray[i].substr(mutationposition[i]-k+1,k)<<"\n"<<baseconverter3(genomearray[i-k])<<endl;
		}
	}
	file.close();


//ofstream file5 (yourname2.c_str());
//for (int i=0;i<length;i++){
//	if(mutarray[i]==0){
//		file5<<wtarray[i]<<"\t"<<mutarray[i]<<endl;
//	}
//	if(mutarray[i]!=0 && mutationposition[i]<16){
//	file5<<wtarray[i]<<"\t"<<mutarray[i]<<"\t"<<sequencearray[i].substr(mutationposition[i],k)<<"\t"<<baseconverter3(genomearray[i-1])<<endl;
//	}
//	if(mutarray[i]!=0 && mutationposition[i]>=16){
//	file5<<wtarray[i]<<"\t"<<mutarray[i]<<"\t"<<sequencearray[i].substr(mutationposition[i]-k+1,k)<<"\t"<<baseconverter3(genomearray[i-k])<<endl;
//	}
//}
//file5.close();

	delete[] genomearray;
	delete[] array2;
	delete[] pArray;
	delete[] wtarray;
	delete[] mutarray;

	cout<<"Done"<<endl;

	system("PAUSE");
	return 0; 
}