<commit_before>package japsadev.bio.hts.barcode;
import java.io.File;
import java.io.IOException;
import java.text.DecimalFormat;
import java.util.ArrayList;
import java.util.Comparator;

import jaligner.matrix.MatrixLoaderException;
import japsa.seq.Alphabet;
import japsa.seq.Sequence;
import japsa.seq.SequenceOutputStream;
import japsa.seq.SequenceReader;
import japsa.util.Logging;

public class BarCodeAnalysis {
	int 	SCAN_WINDOW; 
	double	DIST_THRES,
			SCORE_THRES;
	public static boolean 	print=false,
							strict=false; // both-ends-matching 
	ArrayList<Sequence> barCodesLeft = new ArrayList<Sequence>(); //barcode sequences from left end
	ArrayList<Sequence> barCodesRight = new ArrayList<Sequence>(); //barcode from right end
	Process[] processes;
	int nSamples;
	int barcodeLen;
	SequenceOutputStream[] streamToScaffolder, streamToFile;

	public BarCodeAnalysis(String barcodeFile, String scriptFile) throws IOException{
		ArrayList<Sequence> allSeq = SequenceReader.readAll(barcodeFile, Alphabet.DNA());
		if(strict){
			allSeq.sort(Comparator.comparing(Sequence::getName));
			int i = 0;
			while(i < allSeq.size()){
				barCodesLeft.add(allSeq.get(i++));
				barCodesRight.add(allSeq.get(i++));
			}
		}else{
			barCodesLeft = allSeq;
			for(Sequence seq:barCodesLeft)
				barCodesRight.add(Alphabet.DNA.complement(seq));

		}
		
		nSamples = barCodesLeft.size();

		processes = new Process[nSamples];
		streamToScaffolder = new SequenceOutputStream[nSamples];
		if(print)
			streamToFile = new SequenceOutputStream[nSamples+1]; //unknown sequences included

		String id;
		for(int i=0;i<nSamples;i++){		
			Sequence barCode = barCodesLeft.get(i);

			id = barCode.getName();
			//System.out.println(i + " >" + id + ":" + barCode);

			ProcessBuilder pb = new ProcessBuilder(scriptFile, id)
					.redirectError(new File("log_" + id + ".err"))
					.redirectOutput(new File("log_" + id + ".out"));
			pb.directory(new File(System.getProperty("user.dir")));
			
			processes[i]  = pb.start();

			Logging.info("Job for " + id  + " started");
			streamToScaffolder[i] = new SequenceOutputStream(processes[i].getOutputStream());
			barcodeLen += barCodesLeft.get(i).length();
		}
		
		barcodeLen /= nSamples;
		
		if(print){
			streamToFile = new SequenceOutputStream[nSamples+1]; // plus unknown
			for(int i=0;i<nSamples;i++){		
				streamToFile[i] = SequenceOutputStream.makeOutputStream(barCodesLeft.get(i).getName()+"_clustered.fastq");
			}
			streamToFile[nSamples] = SequenceOutputStream.makeOutputStream("unknown_clustered.fastq");

		}
		
		
		SCAN_WINDOW = barcodeLen * 3;
		SCORE_THRES = .7;
		DIST_THRES = .1;
	}
	
	public void setThreshold(double ident){
		SCORE_THRES = ident;
	}
	/*
	 * Trying to clustering MinION read data into different samples based on the barcode
	 */
	public void clustering(String dataFile) throws IOException, InterruptedException, MatrixLoaderException{
		SequenceReader reader;
		if(dataFile.equals("-"))
			reader = SequenceReader.getReader(System.in);
		else
			reader = SequenceReader.getReader(dataFile);
		Sequence seq;

		Sequence s5, s3;
		final double[] 	lf = new double[nSamples], //left-forward
				lr = new double[nSamples],	//left-reversed
				rr = new double[nSamples], //right-reversed
				rf = new double[nSamples]; //right-forward
		
		jaligner.Alignment 	alignmentLF = new jaligner.Alignment(),
							alignmentLR = new jaligner.Alignment(),
							alignmentRF = new jaligner.Alignment(),
							alignmentRR = new jaligner.Alignment();
		jaligner.Sequence js5,js3, jBarcodeLeft, jBarcodeRight;
		jaligner.Alignment 	bestLeftAlignment = new jaligner.Alignment(),
							bestRightAlignment = new jaligner.Alignment();

//		Sequence barcodeSeq = new Sequence(Alphabet.DNA4(),barcodeLen,"barcode");
//		Sequence tipSeq = new Sequence(Alphabet.DNA4(),SCAN_WINDOW,"tip");
//
//		BarcodeAlignment barcodeAlignment = new BarcodeAlignment(barcodeSeq, tipSeq);

		while ((seq = reader.nextSequence(Alphabet.DNA())) != null){
			if(seq.length() < barcodeLen*4){
				System.err.println("Ignore short sequence " + seq.getName());
				continue;
			}
			//alignment algorithm is applied here. For the beginning, Smith-Waterman local pairwise alignment is used

			s5 = seq.subSequence(0, SCAN_WINDOW);
//			s3=seq.subSequence(seq.length()-SCAN_WINDOW,seq.length());
			s3 = Alphabet.DNA.complement(seq.subSequence(seq.length()-SCAN_WINDOW,seq.length()));



			double bestScore = 0.0;
			double distance = 0.0; //distance between bestscore and the runner-up
			
			int bestIndex = nSamples;

			for(int i=0;i<nSamples; i++){
				Sequence barcodeLeft = barCodesLeft.get(i);
				Sequence barcodeRight = barCodesRight.get(i); //rc of right barcode sequence

//				barcodeAlignment.setBarcodeSequence(barcodeLeft);				
//				barcodeAlignment.setReadSequence(s5);				
//				lf[i]=barcodeAlignment.align();				
//
//				barcodeAlignment.setReadSequence(Alphabet.DNA.complement(s3));
//				lr[i]=barcodeAlignment.align();
//
//				barcodeAlignment.setBarcodeSequence(barcodeRight);
//				barcodeAlignment.setReadSequence(Alphabet.DNA.complement(s3));
//				rr[i]=barcodeAlignment.align();
//				
//				barcodeAlignment.setReadSequence(s5);
//				rf[i]=barcodeAlignment.align();
				js5 = new jaligner.Sequence(s5.toString());
				js3 = new jaligner.Sequence(s3.toString());
				jBarcodeLeft = new jaligner.Sequence(barcodeLeft.toString());
				jBarcodeRight = new jaligner.Sequence(barcodeRight.toString());

				alignmentLF = jaligner.SmithWatermanGotoh.align(js5, jBarcodeLeft, jaligner.matrix.MatrixLoader.load("BLOSUM62"), 10f, 0.5f);
				alignmentLR = jaligner.SmithWatermanGotoh.align(js3, jBarcodeLeft, jaligner.matrix.MatrixLoader.load("BLOSUM62"), 10f, 0.5f);
				alignmentRF = jaligner.SmithWatermanGotoh.align(js5, jBarcodeRight, jaligner.matrix.MatrixLoader.load("BLOSUM62"), 10f, 0.5f);
				alignmentRR = jaligner.SmithWatermanGotoh.align(js3, jBarcodeRight, jaligner.matrix.MatrixLoader.load("BLOSUM62"), 10f, 0.5f);
				
				lf[i] = alignmentLF.getIdentity()/(float)barcodeLeft.length();
				lr[i] = alignmentLR.getIdentity()/(float)barcodeLeft.length();
				rf[i] = alignmentRF.getIdentity()/(float)barcodeRight.length();
				rr[i] = alignmentRR.getIdentity()/(float)barcodeRight.length();
				

				double myScore = 0.0;
				if(strict){
					myScore = Math.max(lf[i] + rr[i] , lr[i] + rf[i])/2;
				}
				else{	
					myScore = Math.max(Math.max(lf[i], lr[i]), Math.max(rf[i], rr[i]));
				}
				
				if (myScore > bestScore){
					//Logging.info("Better score=" + myScore);
					distance = myScore-bestScore;
					bestScore = myScore;		
					bestIndex = i;
					if(strict){
						if(lf[i] + rr[i] > lr[i] + rf[i]){
							bestLeftAlignment = alignmentLF;
							bestRightAlignment = alignmentRR;
						}else{
							bestLeftAlignment = alignmentLR;
							bestRightAlignment = alignmentRF;
						}
						
					}else{
						if(myScore==lf[i] || myScore==rr[i]){
							bestLeftAlignment = alignmentLF;
							bestRightAlignment = alignmentRR;
						}else{
							bestLeftAlignment = alignmentLR;
							bestRightAlignment = alignmentRF;
						}
					}
					
				} else if((bestScore-myScore) < distance){
					distance=bestScore-myScore;
				}
					
			}
			
			
			String retval="";
			DecimalFormat twoDForm =  new DecimalFormat("#.##");
			if(bestScore < SCORE_THRES || distance < DIST_THRES 
//					|| bestLeftAlignment.getLength() < 0.8 * bestLeftAlignment.getOriginalSequence2().length()
//					|| bestRightAlignment.getLength() < 0.8 * bestRightAlignment.getOriginalSequence2().length()
					){
				//Logging.info("Unknown sequence " + seq.getName());
				retval = "unknown:"+Double.valueOf(twoDForm.format(bestScore))+":"+Double.valueOf(twoDForm.format(distance))+"|0-0:0-0|";
				seq.setName(retval + seq.getName());

				if(print)
					seq.print(streamToFile[nSamples]);
			}
			//if the best (sum of both ends) alignment in template sequence is greater than in complement
			else {
//				Logging.info("Sequence " + seq.getName() + " might belongs to sample " + barCodesLeft.get(bestIndex).getName() + " with score=" + bestScore);<commit_msg>primer+tail+barcode work<commit_after>package japsadev.bio.hts.barcode;
import java.io.File;
import java.io.IOException;
import java.text.DecimalFormat;
import java.util.ArrayList;
import java.util.Comparator;

import jaligner.matrix.MatrixLoaderException;
import japsa.seq.Alphabet;
import japsa.seq.Sequence;
import japsa.seq.SequenceOutputStream;
import japsa.seq.SequenceReader;
import japsa.util.Logging;

public class BarCodeAnalysis {
	int 	SCAN_WINDOW; 
	double	DIST_THRES,
			SCORE_THRES;
	public static boolean 	print=false,
							strict=false; // both-ends-matching 
	ArrayList<Sequence> barCodesLeft = new ArrayList<Sequence>(); //barcode sequences from left end
	ArrayList<Sequence> barCodesRight = new ArrayList<Sequence>(); //barcode from right end
	Process[] processes;
	int nSamples;
	int barcodeLen;
	SequenceOutputStream[] streamToScaffolder, streamToFile;

	public BarCodeAnalysis(String barcodeFile, String scriptFile) throws IOException{
		ArrayList<Sequence> allSeq = SequenceReader.readAll(barcodeFile, Alphabet.DNA());
		if(strict){
			allSeq.sort(Comparator.comparing(Sequence::getName));
			int i = 0;
			while(i < allSeq.size()){
				barCodesLeft.add(allSeq.get(i++));
				barCodesRight.add(allSeq.get(i++));
			}
		}else{
			barCodesLeft = allSeq;
			for(Sequence seq:barCodesLeft)
				barCodesRight.add(Alphabet.DNA.complement(seq));

		}
		
		nSamples = barCodesLeft.size();

		processes = new Process[nSamples];
		streamToScaffolder = new SequenceOutputStream[nSamples];
		if(print)
			streamToFile = new SequenceOutputStream[nSamples+1]; //unknown sequences included

		String id;
		for(int i=0;i<nSamples;i++){		
			Sequence barCode = barCodesLeft.get(i);

			id = barCode.getName();
			//System.out.println(i + " >" + id + ":" + barCode);

			ProcessBuilder pb = new ProcessBuilder(scriptFile, id)
					.redirectError(new File("/dev/null"))
					.redirectOutput(new File("/dev/null"));
			pb.directory(new File(System.getProperty("user.dir")));
			
			processes[i]  = pb.start();

			Logging.info("Job for " + id  + " started");
			streamToScaffolder[i] = new SequenceOutputStream(processes[i].getOutputStream());
			barcodeLen += barCodesLeft.get(i).length();
		}
		
		barcodeLen /= nSamples;
		
		if(print){
			streamToFile = new SequenceOutputStream[nSamples+1]; // plus unknown
			for(int i=0;i<nSamples;i++){		
				streamToFile[i] = SequenceOutputStream.makeOutputStream(barCodesLeft.get(i).getName()+"_clustered.fastq");
			}
			streamToFile[nSamples] = SequenceOutputStream.makeOutputStream("unknown_clustered.fastq");

		}
		
		
		SCAN_WINDOW = barcodeLen * 3;
		SCORE_THRES = .7;
		DIST_THRES = .1;
	}
	
	public void setThreshold(double ident){
		SCORE_THRES = ident;
	}
	
	//alignment matrix
	float[][] scores = {
		{  2.7f, -4.5f, -4.5f, -4.5f},
		{ -4.5f,  2.7f, -4.5f, -4.5f},
		{ -4.5f, -4.5f,  2.7f, -4.5f},
		{ -4.5f, -4.5f, -4.5f,  2.7f}			
	};
	
	private jaligner.matrix.Matrix 	matrix = jaligner.matrix.MatrixGenerator.generate(2.7f, -4.5f);
	float openPenalty = 4.7f, extendPenalty = 1.6f;
	/*
	 * Trying to clustering MinION read data into different samples based on the barcode
	 */
	public void clustering(String dataFile) throws IOException, InterruptedException, MatrixLoaderException{
		SequenceReader reader;
		if(dataFile.equals("-"))
			reader = SequenceReader.getReader(System.in);
		else
			reader = SequenceReader.getReader(dataFile);
		Sequence seq;

		Sequence s5, s3;
		final double[] 	lf = new double[nSamples], //left-forward
				lr = new double[nSamples],	//left-reversed
				rr = new double[nSamples], //right-reversed
				rf = new double[nSamples]; //right-forward
		
		jaligner.Alignment 	alignmentLF = new jaligner.Alignment(),
							alignmentLR = new jaligner.Alignment(),
							alignmentRF = new jaligner.Alignment(),
							alignmentRR = new jaligner.Alignment();
		jaligner.Sequence js5,js3, jBarcodeLeft, jBarcodeRight;
		jaligner.Alignment 	bestLeftAlignment = new jaligner.Alignment(),
							bestRightAlignment = new jaligner.Alignment();

//		Sequence barcodeSeq = new Sequence(Alphabet.DNA4(),barcodeLen,"barcode");
//		Sequence tipSeq = new Sequence(Alphabet.DNA4(),SCAN_WINDOW,"tip");
//
//		BarcodeAlignment barcodeAlignment = new BarcodeAlignment(barcodeSeq, tipSeq);

		while ((seq = reader.nextSequence(Alphabet.DNA())) != null){
			if(seq.length() < barcodeLen*4){
				System.err.println("Ignore short sequence " + seq.getName());
				continue;
			}
			//alignment algorithm is applied here. For the beginning, Smith-Waterman local pairwise alignment is used

			s5 = seq.subSequence(0, SCAN_WINDOW);
//			s3=seq.subSequence(seq.length()-SCAN_WINDOW,seq.length());
			s3 = Alphabet.DNA.complement(seq.subSequence(seq.length()-SCAN_WINDOW,seq.length()));



			double bestScore = 0.0;
			double distance = 0.0; //distance between bestscore and the runner-up
			
			int bestIndex = nSamples;

			for(int i=0;i<nSamples; i++){
				Sequence barcodeLeft = barCodesLeft.get(i);
				Sequence barcodeRight = barCodesRight.get(i); //rc of right barcode sequence

//				barcodeAlignment.setBarcodeSequence(barcodeLeft);				
//				barcodeAlignment.setReadSequence(s5);				
//				lf[i]=barcodeAlignment.align();				
//
//				barcodeAlignment.setReadSequence(Alphabet.DNA.complement(s3));
//				lr[i]=barcodeAlignment.align();
//
//				barcodeAlignment.setBarcodeSequence(barcodeRight);
//				barcodeAlignment.setReadSequence(Alphabet.DNA.complement(s3));
//				rr[i]=barcodeAlignment.align();
//				
//				barcodeAlignment.setReadSequence(s5);
//				rf[i]=barcodeAlignment.align();
				js5 = new jaligner.Sequence(s5.toString());
				js3 = new jaligner.Sequence(s3.toString());
				jBarcodeLeft = new jaligner.Sequence(barcodeLeft.toString());
				jBarcodeRight = new jaligner.Sequence(barcodeRight.toString());

			
				alignmentLF = jaligner.SmithWatermanGotoh.align(js5, jBarcodeLeft, matrix, openPenalty, extendPenalty);
				alignmentLR = jaligner.SmithWatermanGotoh.align(js3, jBarcodeLeft, matrix, openPenalty, extendPenalty);
				alignmentRF = jaligner.SmithWatermanGotoh.align(js5, jBarcodeRight, matrix, openPenalty, extendPenalty);
				alignmentRR = jaligner.SmithWatermanGotoh.align(js3, jBarcodeRight, matrix, openPenalty, extendPenalty);
				
				lf[i] = alignmentLF.getIdentity()/(float)Math.max(barcodeLeft.length(),alignmentLF.getLength());
				lr[i] = alignmentLR.getIdentity()/(float)Math.max(barcodeLeft.length(),alignmentLR.getLength());
				rf[i] = alignmentRF.getIdentity()/(float)Math.max(barcodeRight.length(),alignmentRF.getLength());
				rr[i] = alignmentRR.getIdentity()/(float)Math.max(barcodeRight.length(),alignmentRR.getLength());
				

				double myScore = 0.0;
				if(strict){
					myScore = Math.max(lf[i] + rr[i] , lr[i] + rf[i])/2;
				}
				else{	
					myScore = Math.max(Math.max(lf[i], lr[i]), Math.max(rf[i], rr[i]));
				}
				
				if (myScore > bestScore){
					//Logging.info("Better score=" + myScore);
					distance = myScore-bestScore;
					bestScore = myScore;		
					bestIndex = i;
					if(strict){
						if(lf[i] + rr[i] > lr[i] + rf[i]){
							bestLeftAlignment = alignmentLF;
							bestRightAlignment = alignmentRR;
						}else{
							bestLeftAlignment = alignmentLR;
							bestRightAlignment = alignmentRF;
						}
						
					}else{
						if(myScore==lf[i] || myScore==rr[i]){
							bestLeftAlignment = alignmentLF;
							bestRightAlignment = alignmentRR;
						}else if(myScore==lr[i] || myScore==rf[i]){
							bestLeftAlignment = alignmentLR;
							bestRightAlignment = alignmentRF;
						}
					}
					
				} else if((bestScore-myScore) < distance){
					distance=bestScore-myScore;
				}
					
			}
			
			
			String retval="";
			DecimalFormat twoDForm =  new DecimalFormat("#.##");
			if(bestScore < SCORE_THRES || distance < DIST_THRES 
//					|| bestLeftAlignment.getLength() < 0.8 * bestLeftAlignment.getOriginalSequence2().length()
//					|| bestRightAlignment.getLength() < 0.8 * bestRightAlignment.getOriginalSequence2().length()
					){
				//Logging.info("Unknown sequence " + seq.getName());
				retval = "unknown:"+Double.valueOf(twoDForm.format(bestScore))+":"+Double.valueOf(twoDForm.format(distance))+"|0-0:0-0|";
				seq.setName(retval + seq.getName());

				if(print)
					seq.print(streamToFile[nSamples]);
			}
			//if the best (sum of both ends) alignment in template sequence is greater than in complement
			else {
//				Logging.info("Sequence " + seq.getName() + " might belongs to sample " + barCodesLeft.get(bestIndex).getName() + " with score=" + bestScore);
{'fraction_non_alphanumeric': 0.14024975984630164, 'fraction_numerical': 0.009726224783861672, 'mean_word_length': 3.9544913741820347, 'pattern_counts': {'":': 4, '<': 33, '<?xml version=': 0, '>': 19, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3869505', 'n_tokens_mistral': 6650, 'n_tokens_neox': 5495, 'n_words': 992}