<commit_before>/*****************************************************************************
 * Copyright (c) <NAME>, Monash Uni & UQ, All rights reserved.         *
 *                                                                           *
 * Redistribution and use in source and binary forms, with or without        *
 * modification, are permitted provided that the following conditions        *
 * are met:                                                                  * 
 *                                                                           *
 * 1. Redistributions of source code must retain the above copyright notice, *
 *    this list of conditions and the following disclaimer.                  *
 * 2. Redistributions in binary form must reproduce the above copyright      *
 *    notice, this list of conditions and the following disclaimer in the    *
 *    documentation and/or other materials provided with the distribution.   *
 * 3. Neither the names of the institutions nor the names of the contributors*
 *    may be used to endorse or promote products derived from this software  *
 *    without specific prior written permission.                             *
 *                                                                           *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS   *
 * IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, *
 * THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR    *
 * PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR         *
 * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,     *
 * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,       *
 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR        *
 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF    *
 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING      *
 * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS        *
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.              *
 ****************************************************************************/

/**************************     REVISION HISTORY    **************************
 * 18/10/2013 - <NAME>: Created                                        
 *  
 ****************************************************************************/

package japsa.tools.bio.np;

import java.io.File;
import java.io.IOException;
import java.util.ArrayList;
import java.util.HashMap;

import org.slf4j.Logger;
import org.slf4j.LoggerFactory;

import htsjdk.samtools.SAMRecord;
import htsjdk.samtools.SAMRecordIterator;
import htsjdk.samtools.SamReader;
import htsjdk.samtools.SamReaderFactory;
import htsjdk.samtools.ValidationStringency;
import japsa.seq.Alphabet;
import japsa.seq.Sequence;
import japsa.seq.SequenceReader;
import japsa.util.CommandLine;
import japsa.util.deploy.Deployable;
import japsa.bio.hts.scaffold.AlignmentRecord;
import japsa.bio.hts.scaffold.Contig;



@Deployable(
		scriptName = "jsa.np.flankDetect",
		scriptDesc = "Detect flanking sequences from both ends of nanopore reads"
		)
public class FlankSeqsDetectorCmd extends CommandLine{
    private static final Logger LOG = LoggerFactory.getLogger(FlankSeqsDetectorCmd.class);

	public FlankSeqsDetectorCmd(){
		super();
		Deployable annotation = getClass().getAnnotation(Deployable.class);		
		setUsage(annotation.scriptName() + " [options]");
		setDesc(annotation.scriptDesc());

		addString("flankFile",null,"Flank sequences file, maximum 2 sequences",true);
		addString("refFile","","Reference sequences");
		addString("bamFile",null,"Bam file",true);
		addDouble("qual", 1, "Mininum quality");
		addInt("insert", 10, "Minimum length of insert sequence in-between 2 flanking sequences");
		addInt("tips", 20, "Maximum percentage of the overhangs compared to the corresponding flanking sequence");
		addDouble("cover", 80, "Mininum percentage of flank sequence coverage for a valid alignment");

		addStdHelp();	
	}
	/**
	 * @param args
	 */
	public static void main(String[] args) throws IOException{
		/*********************** Setting up script ****************************/		 
		/*********************** Setting up script ****************************/		
		CommandLine cmdLine = new FlankSeqsDetectorCmd();		
		args = cmdLine.stdParseLine(args);
		/**********************************************************************/
		String 	flankSeqsFile= cmdLine.getStringVal("flankFile"),
				refSeqsFile= cmdLine.getStringVal("refFile");
		String bamFile = cmdLine.getStringVal("bamFile");
		double 	qual = cmdLine.getDoubleVal("qual"),
				flkCov = cmdLine.getDoubleVal("cover");
		int insertLength = cmdLine.getIntVal("insert"),
			tipsPercentage = cmdLine.getIntVal("tips");
		
		SequenceReader seqReader = SequenceReader.getReader(flankSeqsFile);
		Sequence seq;
		HashMap<String, Contig> refSeqs = new HashMap<>();
		ArrayList<Contig> flankSeqs=new ArrayList<>();
							
		int index=0;
		while ((seq = seqReader.nextSequence(Alphabet.DNA())) != null)
			flankSeqs.add(new Contig(index++,seq));
		seqReader.close();
		
		if(!refSeqsFile.isEmpty()){
			seqReader = SequenceReader.getReader(refSeqsFile);
			index=0;
			while ((seq = seqReader.nextSequence(Alphabet.DNA())) != null)
				refSeqs.put(seq.getName(), new Contig(index++,seq));
			seqReader.close();
		}
		
		if(flankSeqs.size() > 2){
			System.err.println("More than 2 sequences!");
			System.exit(1);
		}
		
		SamReaderFactory.setDefaultValidationStringency(ValidationStringency.SILENT);
		SamReader reader = SamReaderFactory.makeDefault().open(new File(bamFile));
		SAMRecordIterator iter = reader.iterator();

		SAMRecord curSAMRecord=null;
		AlignmentRecord curAlnRecord=null;
		FlankRecord fr=null;
		String readID = "";

		while (iter.hasNext()) {
	
			try {
				curSAMRecord = iter.next();
			}catch(Exception e) {
				LOG.warn("Ignore one faulty SAM record: \n {}", e.getMessage());
				continue;
			}
			
			if (curSAMRecord.getReadUnmappedFlag() || curSAMRecord.getMappingQuality() < qual || curSAMRecord.isSecondaryAlignment()){		
				LOG.info("Ignore unmapped, low-quality or not primarily map record from {}...", curSAMRecord.getReadName());
				continue;		
			}
								
			if (!readID.equals(curSAMRecord.getReadName())){
				//output prev
				if(fr!=null)
					System.out.println(fr.printJunctionOnRef());
					
				//update for next
				readID = curSAMRecord.getReadName();
				fr=new FlankRecord(readID);
			}
			//adding info to current FlankRecord
			Contig ctg=null;
			String refName = curSAMRecord.getReferenceName();
			if(refSeqs.containsKey(refName)){
				ctg=refSeqs.get(refName);
				curAlnRecord=new AlignmentRecord(curSAMRecord, ctg);
				//increase score if primary
				if(!curSAMRecord.getSupplementaryAlignmentFlag())
					curAlnRecord.score+=100;

				if(curAlnRecord.readAlignmentEnd()-curAlnRecord.readAlignmentStart() < insertLength){
					LOG.info("Ignore integration size too short: {}", curAlnRecord.toString());
					continue;
				}
				if(fr.refRec==null||fr.refRec.score < curAlnRecord.score)
					fr.refRec=curAlnRecord;
				
			}else{
				
				if(flankSeqs.get(0).getName().equals(refName)){
					ctg=flankSeqs.get(0);
					curAlnRecord=new AlignmentRecord(curSAMRecord, ctg);
					if(curAlnRecord.refEnd-curAlnRecord.refStart < (double)flkCov*ctg.length()/100.0)
						continue;
					//not too far from the tip of read
					else if(Math.min(-curAlnRecord.readAlignmentEnd()+curAlnRecord.readLength, curAlnRecord.readAlignmentStart()) > (double)ctg.length()*tipsPercentage/100.0){
						continue;
					}
					if(fr.f0Rec==null||fr.f0Rec.score < curAlnRecord.score)
						fr.f0Rec=curAlnRecord;
					
				}else if(flankSeqs.size()>1 && flankSeqs.get(1).getName().equals(refName)){
					ctg=flankSeqs.get(1);
					curAlnRecord=new AlignmentRecord(curSAMRecord, ctg);
					if(curAlnRecord.refEnd-curAlnRecord.refStart < (double)flkCov*ctg.length()/100.0)
						continue;
					//not too far from the tip of read
					else if(Math.min(-curAlnRecord.readAlignmentEnd()+curAlnRecord.readLength, curAlnRecord.readAlignmentStart()) > (double)ctg.length()*tipsPercentage/100.0){
						continue;
					}
					if(fr.f1Rec==null||fr.f1Rec.score < curAlnRecord.score)
						fr.f1Rec=curAlnRecord;
				}else{
					LOG.info("Flank not found: {} != {}!", refName, flankSeqs.get(0).getName());
					System.exit(1);
				}
					
			}

			
		}// while
		iter.close();
		/**********************************************************************/
	}
	

//	public static void main(String[] args) throws IOException{
//		/*********************** Setting up script ****************************/		 
//		/*********************** Setting up script ****************************/		
//		CommandLine cmdLine = new FlankSeqsDetectorCmd();		
//		args = cmdLine.stdParseLine(args);
//		/**********************************************************************/
//		String flankSeqsFile= cmdLine.getStringVal("flankFile");
//		String bamFile = cmdLine.getStringVal("bamFile");
//		double 	qual = cmdLine.getDoubleVal("qual"),
//				flkCov = cmdLine.getDoubleVal("cover");
//		int insertLength = cmdLine.getIntVal("insert"),
//			tipsPercentage = cmdLine.getIntVal("tips");
//		
//		SequenceReader seqReader = SequenceReader.getReader(flankSeqsFile);
//		Sequence seq;
//		ArrayList<Contig> flankSeqs = new ArrayList<>(); 
//		int index=0;<commit_msg>another change in flankdetect: remove confusing alignments<commit_after>/*****************************************************************************
 * Copyright (c) <NAME>, Monash Uni & UQ, All rights reserved.         *
 *                                                                           *
 * Redistribution and use in source and binary forms, with or without        *
 * modification, are permitted provided that the following conditions        *
 * are met:                                                                  * 
 *                                                                           *
 * 1. Redistributions of source code must retain the above copyright notice, *
 *    this list of conditions and the following disclaimer.                  *
 * 2. Redistributions in binary form must reproduce the above copyright      *
 *    notice, this list of conditions and the following disclaimer in the    *
 *    documentation and/or other materials provided with the distribution.   *
 * 3. Neither the names of the institutions nor the names of the contributors*
 *    may be used to endorse or promote products derived from this software  *
 *    without specific prior written permission.                             *
 *                                                                           *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS   *
 * IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, *
 * THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR    *
 * PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR         *
 * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,     *
 * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,       *
 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR        *
 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF    *
 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING      *
 * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS        *
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.              *
 ****************************************************************************/

/**************************     REVISION HISTORY    **************************
 * 18/10/2013 - <NAME>: Created                                        
 *  
 ****************************************************************************/

package japsa.tools.bio.np;

import java.io.File;
import java.io.IOException;
import java.util.ArrayList;
import java.util.HashMap;

import org.slf4j.Logger;
import org.slf4j.LoggerFactory;

import htsjdk.samtools.SAMRecord;
import htsjdk.samtools.SAMRecordIterator;
import htsjdk.samtools.SamReader;
import htsjdk.samtools.SamReaderFactory;
import htsjdk.samtools.ValidationStringency;
import japsa.seq.Alphabet;
import japsa.seq.Sequence;
import japsa.seq.SequenceReader;
import japsa.util.CommandLine;
import japsa.util.deploy.Deployable;
import japsa.bio.hts.scaffold.AlignmentRecord;
import japsa.bio.hts.scaffold.Contig;



@Deployable(
		scriptName = "jsa.np.flankDetect",
		scriptDesc = "Detect flanking sequences from both ends of nanopore reads"
		)
public class FlankSeqsDetectorCmd extends CommandLine{
    private static final Logger LOG = LoggerFactory.getLogger(FlankSeqsDetectorCmd.class);

	public FlankSeqsDetectorCmd(){
		super();
		Deployable annotation = getClass().getAnnotation(Deployable.class);		
		setUsage(annotation.scriptName() + " [options]");
		setDesc(annotation.scriptDesc());

		addString("flankFile",null,"Flank sequences file, maximum 2 sequences",true);
		addString("refFile","","Reference sequences");
		addString("bamFile",null,"Bam file",true);
		addDouble("qual", 1, "Mininum quality");
		addInt("insert", 10, "Minimum length of insert sequence in-between 2 flanking sequences");
		addInt("tips", 20, "Maximum percentage of the overhangs compared to the corresponding flanking sequence");
		addDouble("cover", 80, "Mininum percentage of flank sequence coverage for a valid alignment");

		addStdHelp();	
	}
	/**
	 * @param args
	 */
	public static void main(String[] args) throws IOException{
		/*********************** Setting up script ****************************/		 
		/*********************** Setting up script ****************************/		
		CommandLine cmdLine = new FlankSeqsDetectorCmd();		
		args = cmdLine.stdParseLine(args);
		/**********************************************************************/
		String 	flankSeqsFile= cmdLine.getStringVal("flankFile"),
				refSeqsFile= cmdLine.getStringVal("refFile");
		String bamFile = cmdLine.getStringVal("bamFile");
		double 	qual = cmdLine.getDoubleVal("qual"),
				flkCov = cmdLine.getDoubleVal("cover");
		int insertLength = cmdLine.getIntVal("insert"),
			tipsPercentage = cmdLine.getIntVal("tips");
		
		SequenceReader seqReader = SequenceReader.getReader(flankSeqsFile);
		Sequence seq;
		HashMap<String, Contig> refSeqs = new HashMap<>();
		ArrayList<Contig> flankSeqs=new ArrayList<>();
							
		int index=0;
		while ((seq = seqReader.nextSequence(Alphabet.DNA())) != null)
			flankSeqs.add(new Contig(index++,seq));
		seqReader.close();
		
		if(!refSeqsFile.isEmpty()){
			seqReader = SequenceReader.getReader(refSeqsFile);
			index=0;
			while ((seq = seqReader.nextSequence(Alphabet.DNA())) != null)
				refSeqs.put(seq.getName(), new Contig(index++,seq));
			seqReader.close();
		}
		
		if(flankSeqs.size() > 2){
			System.err.println("More than 2 sequences!");
			System.exit(1);
		}
		
		SamReaderFactory.setDefaultValidationStringency(ValidationStringency.SILENT);
		SamReader reader = SamReaderFactory.makeDefault().open(new File(bamFile));
		SAMRecordIterator iter = reader.iterator();

		SAMRecord curSAMRecord=null;
		AlignmentRecord curAlnRecord=null;
		FlankRecord fr=null;
		String readID = "";

		while (iter.hasNext()) {
	
			try {
				curSAMRecord = iter.next();
			}catch(Exception e) {
				LOG.warn("Ignore one faulty SAM record: \n {}", e.getMessage());
				continue;
			}
			
			if (curSAMRecord.getReadUnmappedFlag() || curSAMRecord.getMappingQuality() < qual || curSAMRecord.isSecondaryAlignment()){		
				LOG.info("Ignore unmapped, low-quality or not primarily map record from {}...", curSAMRecord.getReadName());
				continue;		
			}
								
			if (!readID.equals(curSAMRecord.getReadName())){
				//output prev
				if(fr!=null)
					System.out.println(fr.printJunctionOnRef());
					
				//update for next
				readID = curSAMRecord.getReadName();
				fr=new FlankRecord(readID);
			}
			//adding info to current FlankRecord
			Contig ctg=null;
			String refName = curSAMRecord.getReferenceName();
			if(refSeqs.containsKey(refName)){
				ctg=refSeqs.get(refName);
				curAlnRecord=new AlignmentRecord(curSAMRecord, ctg);

				if(curAlnRecord.readAlignmentEnd()-curAlnRecord.readAlignmentStart() < insertLength){
					LOG.info("Ignore integration size too short: {}", curAlnRecord.toString());
					continue;
				}
				if(fr.refRec==null||fr.refRec.qual < curAlnRecord.qual)
					fr.refRec=curAlnRecord;
				else if(fr.refRec.qual == curAlnRecord.qual){
					LOG.info("Ignore read with confusing alignment to {}:\n {} \nvs\n {}", refName, curAlnRecord.toString(), fr.refRec);
					fr.refRec=null;
					continue;
				}
				
			}else{
				
				if(flankSeqs.get(0).getName().equals(refName)){
					ctg=flankSeqs.get(0);
					curAlnRecord=new AlignmentRecord(curSAMRecord, ctg);
					if(curAlnRecord.refEnd-curAlnRecord.refStart < (double)flkCov*ctg.length()/100.0)
						continue;
					//not too far from the tip of read
					else if(Math.min(-curAlnRecord.readAlignmentEnd()+curAlnRecord.readLength, curAlnRecord.readAlignmentStart()) > (double)ctg.length()*tipsPercentage/100.0){
						continue;
					}
					if(fr.f0Rec==null||fr.f0Rec.qual < curAlnRecord.qual)
						fr.f0Rec=curAlnRecord;
					else if(fr.f0Rec.qual == curAlnRecord.qual){
						LOG.info("Ignore read with confusing alignment to {}:\n {} \nvs\n {}", refName, curAlnRecord.toString(), fr.refRec);
						fr.f0Rec=null;
						continue;
					}
					
				}else if(flankSeqs.size()>1 && flankSeqs.get(1).getName().equals(refName)){
					ctg=flankSeqs.get(1);
					curAlnRecord=new AlignmentRecord(curSAMRecord, ctg);
					if(curAlnRecord.refEnd-curAlnRecord.refStart < (double)flkCov*ctg.length()/100.0)
						continue;
					//not too far from the tip of read
					else if(Math.min(-curAlnRecord.readAlignmentEnd()+curAlnRecord.readLength, curAlnRecord.readAlignmentStart()) > (double)ctg.length()*tipsPercentage/100.0){
						continue;
					}
					if(fr.f1Rec==null||fr.f1Rec.qual < curAlnRecord.qual)
						fr.f1Rec=curAlnRecord;
					else if(fr.f1Rec.qual == curAlnRecord.qual){
						LOG.info("Ignore read with confusing alignment to {}:\n {} \nvs\n {}", refName, curAlnRecord.toString(), fr.refRec);
						fr.f1Rec=null;
						continue;
					}
				}else{
					LOG.info("Flank not found: {} != {}!", refName, flankSeqs.get(0).getName());
					System.exit(1);
				}
					
			}

			
		}// while
		iter.close();
		/**********************************************************************/
	}
	

//	public static void main(String[] args) throws IOException{
//		/*********************** Setting up script ****************************/		 
//		/*********************** Setting up script ****************************/		
//		CommandLine cmdLine = new FlankSeqsDetectorCmd();		
//		args = cmdLine.stdParseLine(args);
//		/**********************************************************************/
//		String flankSeqsFile= cmdLine.getStringVal("flankFile");
//		String bamFile = cmdLine.getStringVal("bamFile");
//		double 	qual = cmdLine.getDoubleVal("qual"),
//				flkCov = cmdLine.getDoubleVal("cover");
//		int insertLength = cmdLine.getIntVal("insert"),
//			tipsPercentage = cmdLine.getIntVal("tips");
//		
//		SequenceReader seqReader = SequenceReader.getReader(flankSeqsFile);
//		Sequence seq;
//		ArrayList<Contig> flankSeqs = new ArrayList<>(); 
//		int index=0;
{'fraction_non_alphanumeric': 0.1804300086818855, 'fraction_numerical': 0.0061794596803023336, 'mean_word_length': 3.490254528777803, 'pattern_counts': {'":': 0, '<': 33, '<?xml version=': 0, '>': 27, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 24, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3399805', 'n_tokens_mistral': 6596, 'n_tokens_neox': 5544, 'n_words': 1309}