<commit_before>/*
    FASTX-toolkit - FASTA/FASTQ preprocessing tools.
    Copyright (C) 2009  <NAME> (<EMAIL>)

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU Affero General Public License as
    published by the Free Software Foundation, either version 3 of the
    License, or (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU Affero General Public License for more details.

    You should have received a copy of the GNU Affero General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/
#include <cstddef>
#include <cstdlib>
#include <algorithm>
#include <ostream>
#include <iostream>
#include <string>
#include <vector>
#include <string.h>
#include <stdio.h>

#include "sequence_alignment.h"

#include <errno.h>
#include <err.h>

#include <config.h>

#include "fastx.h"
#include "fastx_args.h"


#define MAX_ADAPTER_LEN 100

const char* usage=
"usage: fastx_clipper [-h] [-a ADAPTER] [-D] [-l N] [-n] [-d N] [-c] [-C] [-o] [-v] [-z] [-i INFILE] [-o OUTFILE]\n" \
"Part of " PACKAGE_STRING " by <NAME> (<EMAIL>)\n" \
"\n" \
"   [-h]         = This helpful help screen.\n" \
"   [-a ADAPTER] = ADAPTER string. default is CCTTAAGG (dummy adapter).\n" \
"   [-l N]       = discard sequences shorter than N nucleotides. default is 5.\n" \
"   [-d N]       = Keep the adapter and N bases after it.\n" \
"                  (using '-d 0' is the same as not using '-d' at all. which is the default).\n" \
"   [-c]         = Discard non-clipped sequences (i.e. - keep only sequences which contained the adapter).\n" \
"   [-C]         = Discard clipped sequences (i.e. - keep only sequences which did not contained the adapter).\n" \
"   [-k]         = Report Adapter-Only sequences.\n" \
"   [-n]         = keep sequences with unknown (N) nucleotides. default is to discard such sequences.\n" \
"   [-v]         = Verbose - report number of sequences.\n" \
"                  If [-o] is specified,  report will be printed to STDOUT.\n" \
"                  If [-o] is not specified (and output goes to STDOUT),\n" \
"                  report will be printed to STDERR.\n" \
"   [-z]         = Compress output with GZIP.\n" \
"   [-D]	 = DEBUG output.\n" \
"   [-i INFILE]  = FASTA/Q input file. default is STDIN.\n" \
"   [-o OUTFILE] = FASTA/Q output file. default is STDOUT.\n" \
"\n";

//Default adapter - Dummy sequence
char adapter[MAX_ADAPTER_LEN]="CCTTAAGG";
unsigned int min_length=5;
int discard_unknown_bases=1;
int keep_delta=0;
int discard_non_clipped=0;
int discard_clipped=0;
int show_adapter_only=0;
int debug = 0 ;


//Statistics for verbose report
unsigned int count_input=0 ;
unsigned int count_discarded_too_short=0; // see [-l N] option
unsigned int count_discarded_adapter_at_index_zero=0;  //empty sequences (after clipping)
unsigned int count_discarded_no_adapter_found=0; // see [-c] option
unsigned int count_discarded_adapter_found=0; // see [-C] option
unsigned int count_discarded_N=0; // see [-n]

FASTX fastx;
HalfLocalSequenceAlignment align;

int parse_program_args(int __attribute__((unused)) optind, int optc, char* optarg)
{
	switch(optc) {
		case 'k':
			show_adapter_only=1;
			break;

		case 'D':
			debug++;
			break ;

		case 'c':
			discard_non_clipped = 1;
			break;

		case 'C':
			discard_clipped = 1 ;
			break ;
		case 'd':
			if (optarg==NULL) 
				errx(1, "[-d] parameter requires an argument value");
			keep_delta = strtoul(optarg,NULL,10);
			if (keep_delta<0) 
				errx(1,"Invalid number bases to keep (-d %s)", optarg);
			break;
		case 'a':
			strncpy(adapter,optarg,sizeof(adapter)-1);
			//TODO:
			//if (!valid_sequence_string(adapter)) 
			//	errx(1,"Invalid adapter string (-a %s)", adapter);
			break ;
			
		case 'l':
			if (optarg==NULL) 
				errx(1,"[-l] parameter requires an argument value");
			
			min_length = strtoul(optarg, NULL, 10);
			break;
			
		case 'n':
			discard_unknown_bases = 0 ;
			break;

		default:
			errx(1,"Unknown argument (%c)", optc ) ;

	}
	return 1;
}

int parse_commandline(int argc, char* argv[])
{

	fastx_parse_cmdline(argc, argv, "kDCcd:a:s:l:n", parse_program_args);

	if (keep_delta>0) 
		keep_delta += strlen(adapter);
	return 1;
}

int adapter_cutoff_index ( const SequenceAlignmentResults& alignment_results ) __attribute__ ((const));
int adapter_cutoff_index ( const SequenceAlignmentResults& alignment_results )
{
	#if 0
	int mismatches = alignment_results.mismatches ;
	
	//The adapter(=target) is expected to align from the first base.
	//If the start is not zero (=not aligned from first base),
	//count each skipped base as a mismatch
	mismatches += alignment_results.target_start ;

	//The adapter is expected to align up to the end
	//of the adapter(=target), or the end of the query.
	//If it doesn't, count the un-aligned bases as mismatches
	int missing_from_query_end = (alignment_results.query_size - alignment_results.query_end-1);
	int missing_from_target_end = (alignment_results.target_size - alignment_results.target_end-1);

	int missing_from_end = std::min(missing_from_query_end, missing_from_target_end);
	
	mismatches += missing_from_end ;
	

	 
	std::cout << "Missing from start = " << alignment_results.target_start
		  << " Missing from end = " << missing_from_end
		  << " mismatches = " << mismatches 
		  << std::endl;
	
	if (mismatches > max_mismatches)
		return -1;

	return alignment_results.query_start;
	#endif

	int alignment_size = alignment_results.neutral_matches +
			     alignment_results.matches + 
			     alignment_results.mismatches +
			     alignment_results.gaps ;

	//No alignment at all?
	if (alignment_size==0)
		return -1;

	//Any good alignment at the end of the query
	//(even only a single nucleotide)
	//Example:
	//  The adapter starts with CTGTAG, The Query ends with CT - it's a match.
	if ( alignment_results.query_end == alignment_results.query_size-1
	     &&
	     alignment_results.mismatches == 0 ) {
	     	//printf("--1\n");
		return alignment_results.query_start ;
	}

	if ( alignment_size > 5
	     &&
	     alignment_results.target_start == 0
	     &&
	     (alignment_results.matches * 100 / alignment_size ) >= 75 ) {
	     	//printf("--2\n");
		return alignment_results.query_start ;
	}

	if ( alignment_size > 11 
	     &&
	     (alignment_results.matches * 100 / alignment_size ) >= 80 ) {
	     	//printf("--2\n");
		return alignment_results.query_start ;
	}
<commit_msg>fastx-clipper: add -M option: minimum adapter alignment length.<commit_after>/*
    FASTX-toolkit - FASTA/FASTQ preprocessing tools.
    Copyright (C) 2009  <NAME> (<EMAIL>)

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU Affero General Public License as
    published by the Free Software Foundation, either version 3 of the
    License, or (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU Affero General Public License for more details.

    You should have received a copy of the GNU Affero General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/
#include <cstddef>
#include <cstdlib>
#include <algorithm>
#include <ostream>
#include <iostream>
#include <string>
#include <vector>
#include <string.h>
#include <stdio.h>

#include "sequence_alignment.h"

#include <errno.h>
#include <err.h>

#include <config.h>

#include "fastx.h"
#include "fastx_args.h"


#define MAX_ADAPTER_LEN 100

const char* usage=
"usage: fastx_clipper [-h] [-a ADAPTER] [-D] [-l N] [-n] [-d N] [-c] [-C] [-o] [-v] [-z] [-i INFILE] [-o OUTFILE]\n" \
"Part of " PACKAGE_STRING " by <NAME> (<EMAIL>)\n" \
"\n" \
"   [-h]         = This helpful help screen.\n" \
"   [-a ADAPTER] = ADAPTER string. default is CCTTAAGG (dummy adapter).\n" \
"   [-l N]       = discard sequences shorter than N nucleotides. default is 5.\n" \
"   [-d N]       = Keep the adapter and N bases after it.\n" \
"                  (using '-d 0' is the same as not using '-d' at all. which is the default).\n" \
"   [-c]         = Discard non-clipped sequences (i.e. - keep only sequences which contained the adapter).\n" \
"   [-C]         = Discard clipped sequences (i.e. - keep only sequences which did not contained the adapter).\n" \
"   [-k]         = Report Adapter-Only sequences.\n" \
"   [-n]         = keep sequences with unknown (N) nucleotides. default is to discard such sequences.\n" \
"   [-v]         = Verbose - report number of sequences.\n" \
"                  If [-o] is specified,  report will be printed to STDOUT.\n" \
"                  If [-o] is not specified (and output goes to STDOUT),\n" \
"                  report will be printed to STDERR.\n" \
"   [-z]         = Compress output with GZIP.\n" \
"   [-D]	 = DEBUG output.\n" \
"   [-M N]       = require minimum adapter alignment length of N.\n" \
"                  If less than N nucleotides aligned with the adapter - don't clip it." \
"   [-i INFILE]  = FASTA/Q input file. default is STDIN.\n" \
"   [-o OUTFILE] = FASTA/Q output file. default is STDOUT.\n" \
"\n";

//Default adapter - Dummy sequence
char adapter[MAX_ADAPTER_LEN]="CCTTAAGG";
unsigned int min_length=5;
int discard_unknown_bases=1;
int keep_delta=0;
int discard_non_clipped=0;
int discard_clipped=0;
int show_adapter_only=0;
int debug = 0 ;
int minimum_adapter_length = 0;


//Statistics for verbose report
unsigned int count_input=0 ;
unsigned int count_discarded_too_short=0; // see [-l N] option
unsigned int count_discarded_adapter_at_index_zero=0;  //empty sequences (after clipping)
unsigned int count_discarded_no_adapter_found=0; // see [-c] option
unsigned int count_discarded_adapter_found=0; // see [-C] option
unsigned int count_discarded_N=0; // see [-n]

FASTX fastx;
HalfLocalSequenceAlignment align;

int parse_program_args(int __attribute__((unused)) optind, int optc, char* optarg)
{
	switch(optc) {
		case 'M':
			if (optarg==NULL) 
				errx(1, "[-M] parameter requires an argument value");
			minimum_adapter_length = atoi(optarg);
			if (minimum_adapter_length<=0) 
				errx(1,"Invalid minimum adapter length (-M %s)", optarg);
			break;

		case 'k':
			show_adapter_only=1;
			break;

		case 'D':
			debug++;
			break ;

		case 'c':
			discard_non_clipped = 1;
			break;

		case 'C':
			discard_clipped = 1 ;
			break ;
		case 'd':
			if (optarg==NULL) 
				errx(1, "[-d] parameter requires an argument value");
			keep_delta = strtoul(optarg,NULL,10);
			if (keep_delta<0) 
				errx(1,"Invalid number bases to keep (-d %s)", optarg);
			break;
		case 'a':
			strncpy(adapter,optarg,sizeof(adapter)-1);
			//TODO:
			//if (!valid_sequence_string(adapter)) 
			//	errx(1,"Invalid adapter string (-a %s)", adapter);
			break ;
			
		case 'l':
			if (optarg==NULL) 
				errx(1,"[-l] parameter requires an argument value");
			
			min_length = strtoul(optarg, NULL, 10);
			break;
			
		case 'n':
			discard_unknown_bases = 0 ;
			break;

		default:
			errx(1,"Unknown argument (%c)", optc ) ;

	}
	return 1;
}

int parse_commandline(int argc, char* argv[])
{

	fastx_parse_cmdline(argc, argv, "M:kDCcd:a:s:l:n", parse_program_args);

	if (keep_delta>0) 
		keep_delta += strlen(adapter);
	return 1;
}

int adapter_cutoff_index ( const SequenceAlignmentResults& alignment_results ) __attribute__ ((const));
int adapter_cutoff_index ( const SequenceAlignmentResults& alignment_results )
{
	#if 0
	int mismatches = alignment_results.mismatches ;
	
	//The adapter(=target) is expected to align from the first base.
	//If the start is not zero (=not aligned from first base),
	//count each skipped base as a mismatch
	mismatches += alignment_results.target_start ;

	//The adapter is expected to align up to the end
	//of the adapter(=target), or the end of the query.
	//If it doesn't, count the un-aligned bases as mismatches
	int missing_from_query_end = (alignment_results.query_size - alignment_results.query_end-1);
	int missing_from_target_end = (alignment_results.target_size - alignment_results.target_end-1);

	int missing_from_end = std::min(missing_from_query_end, missing_from_target_end);
	
	mismatches += missing_from_end ;
	

	 
	std::cout << "Missing from start = " << alignment_results.target_start
		  << " Missing from end = " << missing_from_end
		  << " mismatches = " << mismatches 
		  << std::endl;
	
	if (mismatches > max_mismatches)
		return -1;

	return alignment_results.query_start;
	#endif

	int alignment_size = alignment_results.neutral_matches +
			     alignment_results.matches + 
			     alignment_results.mismatches +
			     alignment_results.gaps ;

	//No alignment at all?
	if (alignment_size==0)
		return -1;

	if (minimum_adapter_length>0 && alignment_size<minimum_adapter_length)
		return -1;

	//Any good alignment at the end of the query
	//(even only a single nucleotide)
	//Example:
	//  The adapter starts with CTGTAG, The Query ends with CT - it's a match.
	if ( alignment_results.query_end == alignment_results.query_size-1
	     &&
	     alignment_results.mismatches == 0 ) {
	     	//printf("--1\n");
		return alignment_results.query_start ;
	}

	if ( alignment_size > 5
	     &&
	     alignment_results.target_start == 0
	     &&
	     (alignment_results.matches * 100 / alignment_size ) >= 75 ) {
	     	//printf("--2\n");
		return alignment_results.query_start ;
	}

	if ( alignment_size > 11 
	     &&
	     (alignment_results.matches * 100 / alignment_size ) >= 80 ) {
	     	//printf("--2\n");
		return alignment_results.query_start ;
	}

{'fraction_non_alphanumeric': 0.1221324484201414, 'fraction_numerical': 0.009955273409320444, 'mean_word_length': 3.1580683863227352, 'pattern_counts': {'":': 0, '<': 69, '<?xml version=': 0, '>': 50, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, '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': '2119789', 'n_tokens_mistral': 4988, 'n_tokens_neox': 4498, 'n_words': 1413}