#!/usr/bin/perl
use strict;

#          Fastq adapters clipper/reads renamer/matepair splitter
#          developed by Dr. Markiyan Samborskyy,
#          University of Cambridge, Department of Biochemistry, 
#          DNA Sequencing Facility.
#          ms587 [a.t_mole.bio.cam.ac.uk or v32sw .a_t._litech.lviv.ua
#          27 February 2012 - 30 November 2015
   
#the task is to trim illumina(solexa) fastq format, (as of 2012/02) from miseq or older ones
#in order to allow proper read pairs merging

unless(defined $ARGV[0]){die "Illumina fastq adapters trimmer/matepair linker splitter
	for Truseq, Nextera and Nextera Matepair Libraries
	ver: t16 151130
USAGE: $0 F.fastq_in] [R_fastq_in]
Arguments: 
	-idin=[read_name_begin]   - read ID prefix to search for
	-idout=[read_name_begin]  - read ID prefix to replace with
	-iin	- Interleaved inut
	-iout	- Interleaved output
	-base=[output files base name]
	-out1=[output F file file full name]
	-outlq=[output Low Quality/singleton file file full name (one of the mates had been clipped]
	-out2=[output R file file full name]
	-ends=[ends_naming_style] 1 - _F/_R; 2 - .f/.r; 3 - /1 or /2; 4- no change
	-minl=[minimum non-debris read length]
	-nohproc - do not process fastq headers
	-noclip - do not clip solexa adapter or quality
	-noadclip - do not clip solexa adapter
	-nohpclip - do not clip long homopolymers
	-mpestd -standard matepair processing -> only outer matepair are marked as pe, rest goest to single end files (to reduce some asm tools confusion)
	-ts - truseq adapters clipping
	-nx - nextera shotgun reads ends adapters trimming (Shotgun only!), othervice true seq ones are trimmed.
	-vs=[seq_screen.fasta] remove reads matching screening database
Also it has end adapter autodetection for shotgun, enabled by default - count both adapter types, and report the outcome
Also check for matepair ones too.
autodect id_in if id_out is specified
t10: improved chewed start adapters detection and removal....
t11: add kmer-based phiX/etc? removal
t12: trim qual BEFORE adapter check, select the closest trimpoint
t14: trim oligo C{8} or G{12}  on the R1/R2 (esp if it is more than 2X the # of G's in the read) - remove unballanced crap reads + delete adapter dimer poly A reads
t15: -derep using first/last 64 chr's;
t16: -vs=[screen] - now we can remove phiX (default)/etc reads;
\n";}

my %argv_def; #so we set defaults:
#$argv_def{rd_prefix}="SLXA-HG_2_";#[SLXA run_ID]

$argv_def{drq}=5; #- minimal difference in Q1max and Q2max "double read" filtering tresshold
#           (Q1max/Q2max - qualities from two chanels with the highest signal per each base) 
$argv_def{drbp}=2; #max allowed "dr" basepairs per read
$argv_def{errors_max}=8;#max number of cumulative errors allowed in the read

my $argv=get_ARGV2(\%argv_def);
#general opt
my $intl_in=0;#2 files as the input if 0 (interleaved in)
my $intl_out=0;#2 files on the output (interleaved out)
my $paired=0;
if(defined $argv->{files_list}[1] or $argv->{iin}){#or $argv->{mpe} or $argv->{mpestd} or $argv->{mpe454}
	$paired=1;#paired means there is a second input file or interleaved input sequence
}
if($argv->{iin}){
	$intl_in=1;
}
if($argv->{iout}){
	$intl_out=1;
}

#input files
my $fastq_in_F=$argv->{files_list}[0];
my $fastq_in_R=$argv->{files_list}[1];
my $file_base1=$fastq_in_F;$file_base1=~s/\.fast.+//i;
my $file_base2=$fastq_in_R;$file_base2=~s/\.fast.+//i;
unless(defined $argv->{base}){#set input basename
	if($fastq_in_F=~m/(\S+)\.fas/i){
		$argv->{base}=$1;
	}elsif($fastq_in_F=~m/(\S+)\./i){
		$argv->{base}=$1;
	}else{
		$argv->{base}=$fastq_in_F;
	}
}

if(($paired or $argv->{mpestd}) and not $intl_out){#interleaved output
	$file_base1=$argv->{base}."_F";
}else{#interleaved output
	$file_base1=$argv->{base};
}

#try autodetect library type (if not specified)
unless($argv->{nx} or $argv->{ts} or $argv->{mpe} or $argv->{mpestd}){
	my $lib_config=check_adapter_types($fastq_in_F);
	if($lib_config->{is_nx}){$argv->{nx}=1};
	if($argv->{mpeauto}){
		if($lib_config->{is_mpe}){$argv->{mpestd}=1};
	}
}

#output file names config
my $fout_suffix="_adtrim.fastq";
my $fout_suffix_lq="_adtrim_lq.fastq";
if($argv->{mpe} or $argv->{mpestd}){
	$fout_suffix="_mpes.fastq";
	$fout_suffix_lq="_mpes_se.fastq";#single ended mate pairs
}
my $fastaq_out1=$file_base1.$fout_suffix;
my $fastaq_lq_out=$file_base1.$fout_suffix_lq;
my $fastaq_out2;
if($argv->{out1}){
	$fastaq_out1=$argv->{out1};
}
if($argv->{outlq}){
	$fastaq_lq_out=$argv->{outlq};
}else{#set default for out1 defined case
	$fastaq_lq_out=$argv->{base}.$fout_suffix_lq;
}
#ok, now sort out the second fastq output
if(($paired or $argv->{mpestd}) and not $intl_out){
	if($argv->{base}){
		$file_base2=$argv->{base}."_R";
	}
	$fastaq_out2=$file_base2.$fout_suffix;
	if($argv->{out2}){
		$fastaq_out2=$argv->{out2};
	}
}

my $min_read_size=25;
if($argv->{minl}){
	$min_read_size=$argv->{minl};
}
my $hproc=1;#process fastq headers
if($argv->{nohproc}){
	$hproc=0;
}
my $qual_chr_offset=33;#can be 64 in some older fastaq files
if($argv->{Q64}){$qual_chr_offset=64};
my $fastq = 1;
if(defined $argv->{fastaq}){
	$fastq = 0;
}
my $epc_max=1;#filtering on by default
if(defined $argv->{epc}){
    $epc_max=$argv->{epc};
}
my $id_search;#="M00235:3:000000000-A0H5T";#input read name section


my $id_replace;#="DSM4137_slx_2";#replace with

if(defined $argv->{idin}){
	$id_search = $argv->{idin};
}
if(defined $argv->{idout}){
	$id_replace = $argv->{idout};
}
my $read_end_naming=3;#1 - _F/_R;
#2 .f/.r	- mira original
#3 /1 /2	- illumina /newbler/velvet
#4 1:N:0:1
#5 no change 
if(defined $argv->{ends}){#ends naming type
	$read_end_naming = $argv->{ends};
}

my @fq_fh_in;#input filehandles
my @fq_fh_out;#outout fh
my @fq_lq_fh_out;#outout fh

unless(-f $fastq_in_F){die "Usage: $0 [fastq_in file1] [fastq_in file2]* \n";}

if(defined $id_search){
	print STDERR "Replacing in headers: $id_search by $id_replace\n";
}
my $selfmatch_trim_length=16;#adapter selfmatch length
#make the linker_db - ref to hash of words
my @linkers;
unless($argv->{word_size}){$argv->{word_size}=10};
if($argv->{mpe454}){#454 PE fastq file
	@linkers=("GTTGGAACCGAAAGGGTTTGAATTCAAACCCTTTCGGTTCCAAC",
			 "TCGTATAACTTCGTATAATGTATGCTATACGAAGTTATTACG",
			 "CGTAATAACTTCGTATAGCATACATTATACGAAGTTATACGA");
}else{@linkers=("CTGTCTCTTATACACATCTAGATGTGTATAAGAGACAG")}#CTGTCTCTTATACACATCT, AGATGTGTATAAGAGACAG
#rd_end adapters:GATCGGAAGAGCACACGTCTGAACTCCAGTCAC, GATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT
my $linker_db = make_seq_words(\@linkers, $argv->{word_size});
my $vector_db;
if(defined $argv->{vs}){
	$vector_db = make_vec_words($argv->{vs}, 31);
}
	
my %fastq_stat;#statistics about the fastq files and mate pairs(TODO)

if($fastq_in_F=~m/\.gz$/i or $fastq_in_F=~m/\.Z$/i){
	open (FASTQ_IN1, "gzip -c -d $fastq_in_F |") or die "\nUnable to open gzipped input file: $fastq_in_F\n";
}elsif($fastq_in_F=~m/\.bz2$/i){
	open (FASTQ_IN1, "bzip2 -c -d $fastq_in_F |") or die "\nUnable to open bzipped input file: $fastq_in_F\n";
}else{
	open (FASTQ_IN1, $fastq_in_F) or die "\nUnable to open fastq input file: $fastq_in_F\n";
}
if($intl_in){#same file input
#	FASTQ_IN2=FASTQ_IN1;
}elsif($paired and not $intl_in){
	if($fastq_in_R=~m/\.gz$/i or $fastq_in_R=~m/\.Z$/i){
		open (FASTQ_IN2, "gzip -c -d $fastq_in_R |") or die "\nUnable to open gzipped input file: $fastq_in_R\n";
	}elsif($fastq_in_R=~m/\.bz2$/i){
		open (FASTQ_IN2, "bzip2 -c -d $fastq_in_R |") or die "\nUnable to open bzipped input file: $fastq_in_R\n";
	}elsif(-f $fastq_in_R){
		open (FASTQ_IN2, $fastq_in_R) or die "\nUnable to open fastq input file: $fastq_in_R\n";
	}
}



open (FASTQ_OUT1, ">".$fastaq_out1) or die "Unable to open fastaQ_R1 output file: $fastaq_out1\n";
open (FASTQ_LQ_OUT, ">".$fastaq_lq_out) or die "Unable to open fastaQLQ output file: $fastaq_lq_out\n";
if(($paired or $argv->{mpestd}) and not ($argv->{mpe} or $intl_out)){
	open (FASTQ_OUT2, ">".$fastaq_out2) or die "Unable to open fastaQ_R2 output file: $fastaq_out2\n";
}

my %derep1; my %derep2;#for dereplication;
my ($template1, $head1, $template2, $head2, $trim1, $trim2, $trim1Q, $trim2Q, $seq1PA, $seq1PAC,
	$seq1, $qual_header, $qual1, $seq2, $qual2, $seq1t, $seq2t, $qual1t, $qual2t,
	$seq2startRC, $seq2startRC2);
#1. get 2 reads
#2. if a couple - write both to the headers, if not - only one of them.
#fasta header format if paired:
#>[head1] [seq1] # [seqqual1] # # [head2] [seq2] # [seqqual] # #
#[seq1]

while(<FASTQ_IN1>){
	#print "READING: $_";
	chomp;
	if(m/^@(.+)/){
		$head1=$1;
		chomp($seq1=uc(readline(FASTQ_IN1)));
		$qual_header=readline(FASTQ_IN1);
		chomp($qual1=readline(FASTQ_IN1));
		($head1, $template1)= proc_fastq_header($head1, $id_search, $id_replace, $read_end_naming, $hproc);
		if($paired){
			#get the other end:
			if($intl_in){#interleaved input -read the same file
				chomp($head2=substr(readline(FASTQ_IN1), 1));
				chomp($seq2=uc(readline(FASTQ_IN1)));
				$qual_header=readline(FASTQ_IN1);
				chomp($qual2=readline(FASTQ_IN1));
			}else{#read other file
				chomp($head2=substr(readline(FASTQ_IN2), 1));
				chomp($seq2=uc(readline(FASTQ_IN2)));
				$qual_header=readline(FASTQ_IN2);
				chomp($qual2=readline(FASTQ_IN2));
			}
			if($hproc){
				($head2, $template2)= proc_fastq_header($head2, $id_search, $id_replace, $read_end_naming, $hproc);
			}
			#now some processing for paired ends:
		}
		#copy to trimming buffer
		$seq1t=$seq1;
		$qual1t=$qual1;
		$trim1Q=length($seq1);
		if($paired){
			$seq2t=$seq2;
			$qual2t=$qual2;
			$trim2Q=length($seq2);
		}
		#do the qual check - epc>1 - cut off @that
		
		unless($argv->{noqclip} or $argv->{noclip}){
			$trim1Q = get_fastq_epc_trim(\$seq1t, \$qual1t, $qual_chr_offset, $epc_max, 10);
			if($paired){
				$trim2Q = get_fastq_epc_trim(\$seq2t, \$qual2t, $qual_chr_offset, $epc_max, 10);
			}
		}
		
		
		#get the sequencing adapter matches&trim
		unless($argv->{noadclip} or $argv->{noclip}){
			if($paired){#for adapter trimming
				#if there is a RC match of the first 16 bp of the reverse read in the forward read - set the trimpoint to it
				$seq2startRC=reverse(substr($seq2,0,$selfmatch_trim_length));$seq2startRC=~tr/GATC/CTAG/;
				#and if the read start contains errors - try next 16 bp
				$seq2startRC2=reverse(substr($seq2,$selfmatch_trim_length,$selfmatch_trim_length));$seq2startRC2=~tr/GATC/CTAG/;
			}
			
			if($argv->{nx} and not ($argv->{mpe} or $argv->{mpestd})){#for nextera CTGTCTCTTATACACATCT
				if($seq2=~m/A{16,}$/ or $seq2=~m/C{16,}$/){
					$trim2=$-[0];
				}
				if($seq1=~m/CTGTCTCTTATACACATCT/ or $seq2=~m/CTGTCTCTTATACACATCT/){#we have adapter
					$trim1=$-[0]; 
					if($trim2>$-[0] or not defined $trim2){$trim2=$-[0]};
				}elsif($paired and $seq1=~m/$seq2startRC/){ #thats for unknown adapters
						$trim1=$-[0]+$selfmatch_trim_length; $trim2=$-[0]+$selfmatch_trim_length;
				}elsif($paired and $seq1=~m/$seq2startRC2/){ #thats for unknown adapters if the first word has errors...
					$trim1=$-[0]+($selfmatch_trim_length*2); $trim2=$-[0]+($selfmatch_trim_length*2);
				}elsif($seq1t=~m/CTGTCTCTTATACACATC$/ or $seq1t=~m/CTGTCTCTTATACACAT$/ or $seq1t=~m/CTGTCTCTTATACACA$/ or $seq1t=~m/CTGTCTCTTATACAC$/ or $seq1t=~m/CTGTCTCTTATACA$/ or $seq1t=~m/CTGTCTCTTATAC$/ or
					$seq1t=~m/CTGTCTCTTATA$/ or $seq1t=~m/CTGTCTCTTAT$/ or $seq1t=~m/CTGTCTCTTA$/ or $seq1t=~m/CTGTCTCTT$/ or $seq1t=~m/CTGTCTCT$/ or $seq1t=~m/CTGTCTC$/ or $seq1t=~m/CTGTCT$/){
					$trim1=$-[0]; $trim2=$-[0];
				}elsif($paired and ($seq2t=~m/CTGTCTCTTATACACATC$/ or $seq2t=~m/CTGTCTCTTATACACAT$/ or $seq2t=~m/CTGTCTCTTATACACA$/ or $seq2t=~m/CTGTCTCTTATACAC$/ or $seq2t=~m/CTGTCTCTTATACA$/ or $seq2t=~m/CTGTCTCTTATAC$/ or
					$seq2t=~m/CTGTCTCTTATA$/ or $seq2t=~m/CTGTCTCTTAT$/ or $seq2t=~m/CTGTCTCTTA$/ or $seq2t=~m/CTGTCTCTT$/ or $seq2t=~m/CTGTCTCT$/ or $seq2t=~m/CTGTCTC$/ or $seq2t=~m/CTGTCT$/)){
					$trim1=$-[0]; $trim2=$-[0];
				}elsif($seq1=~m/TCTCCGAGCCCA/){#fwd w 16bp offset
					if($-[0]>16){$trim1=$-[0]-16; $trim2=$-[0]-16;}else{$trim1=0;$trim2=0;}
				}elsif($seq2=~m/TCTGACGCTGCCGA/){#rev w 16bp offset
					if($-[0]>16){$trim1=$-[0]-16; $trim2=$-[0]-16;}else{$trim1=0;$trim2=0;}
				}
			}else{#not nextera PE adapters $argv->{ts}
				#check for adapter
				if($seq1=~m/AGATCGGAAGAGC/ or $seq2=~m/AGATCGGAAGAGC/){#we have adapter
					$trim1=$-[0]; $trim2=$-[0];
				}elsif($seq1=~m/GGAAGAGCACACGTCTGAAC/ or $seq2=~m/GGAAGAGCACACGTCTGAAC/){#we have truseq with a few bp clipped from the front...
					if($-[0]<5){$trim1=0; $trim2=0;#fully clipped
					}else{$trim1=$-[0] - 5; $trim2=$-[0] - 5}
				}elsif($seq1=~m/TCGTATGCCGTCTTCTGC/ or $seq2=~m/TCGTATGCCGTCTTCTGC/){#we have truseq adapter tail - estimate the start...
					if($-[0]<41){$trim1=0; $trim2=0;#fully clipped
					}else{$trim1=$-[0] - 41; $trim2=$-[0] - 41}
				}elsif($argv->{ex} and ($seq1=~m/GATCGGAAGAGCA/ or $seq1=~m/ATCGGAAGAGCAC/ or $seq2=~m/GATCGGAAGAGCGTC/ or $seq2=~m/ATCGGAAGAGCGTCG/)){#we have adapter of 1-2 bases clipped from 5' end by some exonuclease #extra clipping
					$trim1=$-[0]; $trim2=$-[0];
				}elsif($paired and $seq1=~m/$seq2startRC/){ #thats for 100/150 bp H/miSEQ
					$trim1=$-[0]+$selfmatch_trim_length; $trim2=$-[0]+$selfmatch_trim_length;
				}elsif($paired and $seq1=~m/$seq2startRC2/){ #thats for unknown adapters if the first word has errors...
					$trim1=$-[0]+($selfmatch_trim_length*2); $trim2=$-[0]+($selfmatch_trim_length*2);
				}elsif($seq1t=~m/AGATCGGAAGAG$/ or $seq1t=~m/AGATCGGAAGA$/ or $seq1t=~m/AGATCGGAAG$/ or $seq1t=~m/AGATCGGAA$/ or $seq1t=~m/AGATCGGA$/ or $seq1t=~m/AGATCGG$/ or $seq1t=~m/AGATCG$/){
					$trim1=$-[0]; $trim2=$-[0];
				}elsif($seq1=~m/^.{53,65}TTGAAAAAAAA/ or $seq2=~m/^.{45,150}CATTAAAAAAAAA/){
					
					$trim1=0; $trim2=0;#adapter dimer/errors (no insert) - delete
				}  
			}
		}
		#get the G/C homopolymers if the LQ...before trimpoint
		unless($argv->{nohpclip} or $argv->{noclip}){
			#now LQ/HP trim (if needed);
			if($seq1t=~m/G{12}/ or $seq1t=~m/C{8}/){
				if($qual1t=~m/\#{6}/){#if -q33
					if($trim1>$-[0]){$trim1=$-[0]}
					#$seq1t=substr($seq1, 0, $trim1);
					#$qual1t=substr($qual1, 0, $trim1);
				}
				#check the remainer for an unballanced LQ read
				LQ_GC_clip(\$seq1t, \$qual1t, \$trim1, $qual_chr_offset);
			}
			if($seq2t=~m/G{12}/ or $seq2t=~m/C{8}/){
				if($qual2t=~m/\#{6}/){#if -q35 or so
					if($trim2>$-[0]){$trim2=$-[0]}
					#$seq2t=substr($seq2, 0, $trim2);
					#$qual2t=substr($qual2, 0, $trim2);
				}
				LQ_GC_clip(\$seq2t, \$qual2t, \$trim2, $qual_chr_offset);
			}
		}
		#adapter trimming done... -> check results
		if(defined $trim1 and $trim1<$trim1Q){
			$seq1t=substr($seq1, 0, $trim1);
			$qual1t=substr($qual1, 0, $trim1);
		}
		if(defined $trim2 and $trim2<$trim2Q){
			$seq2t=substr($seq2, 0, $trim2);
			$qual2t=substr($qual2, 0, $trim2);
		}
		
		if($argv->{vs}){#control/vector screen
			check_remove_vec($vector_db, 31, \$seq1t, \$qual1t, \$seq2t, \$qual2t);
		}
		#dereplicate (needs quite a bit of RAM for the larger fastq files
		if($argv->{derep}){
			my $seq1h=substr($seq1t,0,64);
			if($paired){
				my $seq2h=substr($seq2t,0,64);
				if($derep1{$seq1h} and $derep2{$seq2h}){
					$seq1t=""; $qual1t="";$seq2t="";$qual2t="";
				}else{$derep1{$seq1h}=1; $derep2{$seq2h}=1;}
			}else{
				if($derep1{$seq1h}){
					$seq1t=""; $qual1t="";
				}else{$derep1{$seq1h}=1};
			}
		}

		#get the length of the remains...
		$trim1=length($seq1t);
		if($paired){$trim2=length($seq2t)}else{$trim2=0};
		
		#detect and split the mate pais if told to do so
		if($argv->{mpe}){#we have mate pair reads - split them
			my ($names,$seqs,$quals,$si,$rds,$fqoutbuf);#refs and index
			if($paired){#both ends reads to check for linker
				#linker_db, wordsz, min_rd, paired, template, #seqs, #quals, #names
				($names,$seqs,$quals)=split_mates($linker_db, $argv->{word_size}, $min_read_size, $paired, $template1, $read_end_naming, $seq1t, $qual1t,$seq2t, $qual2t);
			}else{#single read/or already joined
				($names,$seqs,$quals)=split_mates($linker_db, $argv->{word_size}, $min_read_size, $paired, $template1, $read_end_naming, $seq1t, $qual1t);
			}
			$rds=0;#reads obtained
			$fqoutbuf="";#output buffer
			#dump to fasq, if 1x end -> goes to LQ
			#also have a option of outputing secondary F2/R2 reads to separate file
			
			if($names->[1] and $names->[3]){#set F2 
			#actually better to set it to the read end bit, due to LQ issues
			#PS: add option to dell them
				$names->[1].="2";
			}elsif($names->[0] and $names->[2]){#set R2
				$names->[2].="2";
			}
			foreach $si (3,1,0,2){#first F, than R's
			#for($si=4;$si>-1;$si--){
				if($seqs->[$si]){
					$fqoutbuf.="\@".$names->[$si]."\n$seqs->[$si]\n+\n$quals->[$si]\n";
					$rds++;
				}
			}
			if($rds>1){
				syswrite(FASTQ_OUT1, $fqoutbuf) or die "Unable to write to out1: $fastaq_out1\n";
			}elsif($rds==1){
				syswrite(FASTQ_LQ_OUT, $fqoutbuf) or die "Unable to write to outlq: $fastaq_lq_out\n";
			}
		}elsif($argv->{mpestd}){#standard output into 3x files, with two paired reads and one unpaired one
			my ($names,$seqs,$quals,$si,$rds,$fqoutbuf1,$fqoutbuf2,$fqoutbuflq);#refs and index
			if($paired){#both ends reads to check for linker
				#linker_db, wordsz, min_rd, paired, template, #seqs, #quals, #names
				($names,$seqs,$quals)=split_mates($linker_db, $argv->{word_size}, $min_read_size, $paired, $template1, $read_end_naming, $seq1t, $qual1t,$seq2t, $qual2t);
			}else{#single read/or already joined
				($names,$seqs,$quals)=split_mates($linker_db, $argv->{word_size}, $min_read_size, $paired, $template1, $read_end_naming, $seq1t, $qual1t);
			}
			$rds=0;#reads obtained
			#$fqoutbuf1="";$fqoutbuf2="";#output buffers
			#dump to fasq, if 1x end -> goes to LQ
			#outputing secondary F2/R2 reads to separate file
			if(($seqs->[1] or $seqs->[3]) and ($seqs->[0] or $seqs->[2])){#we have both ends
				if($seqs->[3] and not $seqs->[1]){#F1 only
					$fqoutbuf1.="\@".$names->[3]."\n$seqs->[3]\n+\n$quals->[3]\n";
				}elsif($seqs->[1] and not $seqs->[3]){#F2 only
					$fqoutbuf1.="\@".$names->[1]."\n$seqs->[1]\n+\n$quals->[1]\n";
				}else{#F1&F2
					$fqoutbuf1.="\@".$names->[3]."\n$seqs->[3]\n+\n$quals->[3]\n";
					$fqoutbuflq.="\@".$template1."_F2\n$seqs->[1]\n+\n$quals->[1]\n";
				}
				if($seqs->[0] and not $seqs->[2]){#R1 only
					$fqoutbuf2.="\@".$names->[0]."\n$seqs->[0]\n+\n$quals->[0]\n";
				}elsif($seqs->[2] and not $seqs->[0]){#F2 only
					$fqoutbuf2.="\@".$names->[2]."\n$seqs->[2]\n+\n$quals->[2]\n";
				}else{#R1&R2
					$fqoutbuf2.="\@".$names->[0]."\n$seqs->[0]\n+\n$quals->[0]\n";
					$fqoutbuflq.="\@".$template1."_R2\n$seqs->[2]\n+\n$quals->[2]\n";
				}
			}else{#we have only one end -> all to se
				if(defined $seqs->[3] and length $seqs->[3]){
					$fqoutbuflq.="\@".$template1."_F\n$seqs->[3]\n+\n$quals->[3]\n";
				}
				if(defined $seqs->[1] and length $seqs->[1]){
					$fqoutbuflq.="\@".$template1."_F2\n$seqs->[1]\n+\n$quals->[1]\n";
				}
				if(defined $seqs->[0] and length $seqs->[0]){
					$fqoutbuflq.="\@".$template1."_R\n$seqs->[0]\n+\n$quals->[0]\n";
				}
				if(defined $seqs->[2] and length $seqs->[2]){
					$fqoutbuflq.="\@".$template1."_R2\n$seqs->[2]\n+\n$quals->[2]\n";
				}
			}
			
			if($fqoutbuf1 and $fqoutbuf2){#write normal pair
				syswrite(FASTQ_OUT1, $fqoutbuf1) or die "Unable to write to out1: $fastaq_out1\n";
				if($intl_out){#interleved output
					syswrite(FASTQ_OUT1, $fqoutbuf2) or die "Unable to write to out1(R2): $fastaq_out1\n";
				}else{
					syswrite(FASTQ_OUT2, $fqoutbuf2) or die "Unable to write to out2: $fastaq_out2\n";
				}
			}
			if($fqoutbuflq){
				syswrite(FASTQ_LQ_OUT, $fqoutbuflq) or die "Unable to write to outlq: $fastaq_lq_out\n";
			}
		}elsif($argv->{glue} and $paired){ #"non-homologus ends joining" (resque non-overlaping amp reads)
			my $j_seq="NNNNN"; my $j_qual="#####";
				$seq2t=reverse($seq2t);
				$seq2t=~tr/ATCGatcg/TAGCtagc/;
				$qual2t=reverse($qual2t);
			if($trim1<$min_read_size and $trim2<$min_read_size){
			#$template1
				syswrite(FASTQ_OUT1, "\@$template1\n$seq1t$j_seq$seq2t\n+\n$qual1t$j_qual$qual2t\n") or die "Unable to write to out1: $fastaq_out1\n";
			}else{
				syswrite(FASTQ_LQ_OUT, "\@$template1\n$seq1t$j_seq$seq2t\n+\n$qual1t$j_qual$qual2t\n") or die "Unable to write to outlq: $fastaq_lq_out\n";
			}
		}else{#dump single ended
			#if one of the bits that remains is less, than rd minsize - than output to debris file, which would be interleaved partially
			if(($trim1<$min_read_size) or ($paired and $trim2<$min_read_size)){
				if($trim1>=$min_read_size){syswrite(FASTQ_LQ_OUT, "\@$head1\n$seq1t\n+\n$qual1t\n") or die "Unable to write to outlq: $fastaq_lq_out\n";}
				if($paired and $trim2>=$min_read_size){syswrite(FASTQ_LQ_OUT, "\@$head2\n$seq2t\n+\n$qual2t\n") or die "Unable to write to: $fastaq_lq_out\n"};
			}else{#and in the normal case
				syswrite(FASTQ_OUT1, "\@$head1\n$seq1t\n+\n$qual1t\n") or die "Unable to write to out1: $fastaq_out1\n";
				if($paired){
					if($intl_out){syswrite(FASTQ_OUT1, "\@$head2\n$seq2t\n+\n$qual2t\n") or die "Unable to write to out1: $fastaq_out1\n";
					}else{
						syswrite(FASTQ_OUT2, "\@$head2\n$seq2t\n+\n$qual2t\n") or die "Unable to write to out2: $fastaq_out2\n";
					}
				}
			}
		}
	}
	
	undef($template1);undef($head1);undef($seq1);undef($qual1);undef($seq1t);undef($qual1t);undef($trim1);undef($trim1Q);undef($seq1PA);
	undef($template2);undef($head2);undef($seq2);undef($qual2);undef($seq2t);undef($qual2t);undef($trim2);undef($trim2Q);undef($seq1PAC);
}
close (FASTQ_IN1);
if($paired and not $intl_in){close (FASTQ_IN2)};
close(FASTQ_LQ_OUT);

#close (FASTA_OUT);
#close (FASTAQUAL_OUT);
close (FASTAQ_OUT);

sub check_adapter_types {#check for the adapters in the input dataset
	my $fastq_in_F=$_[0];
	my %lib_config;
	$lib_config{nx_warn_treshold}=0.01;#if above warn about truseq adapters in the nx dataset (or vice-versa); 
	$lib_config{mpe_treshold}=0.01;#if above - consider data to be mpe
	$lib_config{is_mpe}=0;$lib_config{is_nx}=0;#library is assumed to be truseq shotgun by default
	##0-truseq, #1-nexterasg, #2-nextera_mpe_full
	my @adapters = ("AGATCGGAAGAGC","CTGTCTCTTATACACA","CTGTCTCTTATACACATCTAGATGTGTATAAGAGACAG"); #AGATGTGTATAAGAGACAG
	my @adapters_count = (0,0,0);
	my $reads_count=0;
	if($fastq_in_F=~m/\.gz$/i or $fastq_in_F=~m/\.Z$/i){
		open (FASTQ_IN1, "gzip -c -d $fastq_in_F |") or die "\nUnable to open gzipped input file: $fastq_in_F\n";
	}elsif($fastq_in_F=~m/\.bz2$/i){
		open (FASTQ_IN1, "bzip2 -c -d $fastq_in_F |") or die "\nUnable to open bzipped input file: $fastq_in_F\n";
	}else{
		open (FASTQ_IN1, $fastq_in_F) or die "\nUnable to open fastq input file: $fastq_in_F\n";
	}
	print "Trying to detect the input file adapters type: $fastq_in_F\n";
	my ($head1, $seq1, $qual_header, $qual1);
	while(<FASTQ_IN1>){
		if(m/^@(.+)/){
			$head1=$1;
			chomp($seq1=uc(readline(FASTQ_IN1)));
			$qual_header=readline(FASTQ_IN1);
			chomp($qual1=readline(FASTQ_IN1));
			for(my $i=0; $i<=$#adapters; $i++){
				if($seq1=~m/$adapters[$i]/){
					$adapters_count[$i]++;
				}
			}
			$reads_count++;
		}
	}#EOF
	$lib_config{mpe_linkers}=$adapters_count[2];
	
	print "File: $fastq_in_F\tReads: $reads_count\tTruseq: $adapters_count[0]\tNextera_sg: $adapters_count[1]\tMpe: $adapters_count[2]\n";
	if($reads_count){
		$lib_config{mpe_linkers_ratio}=$lib_config{mpe_linkers}/$reads_count;
		$lib_config{nx_adapter_ratio}=$adapters_count[1]/$reads_count;
		$lib_config{ts_adapter_ratio}=$adapters_count[0]/$reads_count;
		print "Truseq_ratio: $lib_config{ts_adapter_ratio}\tNextera_sg_ratio: $lib_config{nx_adapter_ratio}\tMatepair_ratio: $lib_config{mpe_linkers_ratio}\n";
		if($lib_config{mpe_linkers_ratio}>$lib_config{mpe_treshold}){
			$lib_config{is_mpe}=1;
			print "Matepair nextera linkers detected!\n";
		}elsif($lib_config{nx_adapter_ratio}>$lib_config{mpe_treshold}){
			if($lib_config{nx_adapter_ratio}>$lib_config{ts_adapter_ratio}){
				$lib_config{is_nx}=1;
				print "Nextera shotgun adapters detected!\n";
			}
		}
	}
	return(\%lib_config);
}

#separate template calling from the header proc or have switch to existing sub?
sub proc_fastq_header { #0 - header, 1 - $id_search, 2 - id_replace, 3 - end naming type, 
	#(\$head, $id_search, $id_replace, $read_end_naming);
	#header preffix susbtitution
	my $template;
	my $head=$_[0];#ref to head
	my $idin=$_[1];
	my $idout=$_[2];
	if(defined $idout){
		if(defined $idin and $idin=~m/\S/){
			#print STDERR "Replacing in $$head_ref: $idin by $idout\n";
			$head=~s/$idin/$idout/s;
			$head=~s/\:/_/g;
		}elsif($head=~m/(\S+?:\d+:\S+?:\d+):\d+:\d+:\d+/){ #convert read preffix to something userfull
			my $match=quotemeta($1);
			$head=~s/$match/$idout/s;
			$head=~s/\:/_/g;
		}elsif($head=~m/(\S+?)_\d+_\d+_\d+_\d+_\d+/){ #convert read preffix to something userfull
			my $match=quotemeta($1);
			$head=~s/$match/$idout/s;
		};
		
	}
	my $headp=$head;#deref head
	my $end;
	my $rd_ns=$_[3];#reads ouput naming schema
	my $hproc=$_[4];
	#read ends (re)naming according to desired schema
	#detect schema/get the direction of the read
	if($headp=~m|^(\S+?)/1|){
		$template=$1;$end=1;
	}elsif($headp=~m|^(\S+?)/2|){
		$template=$1;$end=2;
	}elsif($headp=~m|^(\S+?) 1_|){
		$template=$1;$end=1;
	}elsif($headp=~m|^(\S+?) 2_|){
		$template=$1;$end=2;
	}elsif($headp=~m|^(\S+?)_F$| or $headp=~m|^(\S+?)_F |){
		$template=$1;$end=1;
	}elsif($headp=~m|^(\S+?)_R$| or $headp=~m|^(\S+?)_R |){
		$template=$1;$end=2;
	}elsif($headp=~m|^(\S+?)\.f$|){
		$template=$1;$end=1;
	}elsif($headp=~m|^(\S+?)\.r$|){
		$template=$1;$end=2;
	}else{#no pair end info
		$template=$head;
	}
	#output suffixes:
	my @FWD_suffixes=("", "_F", ".f", "/1", " 1:0:N:1");
	my @REV_suffixes=("", "_R", ".r", "/2", " 2:0:N:1");
	$headp=$template;
	if($end==1){#forward read
		$headp.=$FWD_suffixes[$rd_ns];
	}elsif($end==2){#reverse read
		$headp.=$REV_suffixes[$rd_ns];
	}
	if($hproc){
		$head=$headp;
		#putting back head: $
	}
# 	unless(defined $template){
# 		$template=$head;
# 	}
	return($head, $template);
}
# sub get_fastq_seqqual {#we need fastq data from the filehandle
#     my $fh=$_[0];
#     my $seq=
# }

sub get_fastq_epc_trim {
    #we need to get the trimpoint in the fastq quality using epc sliding window
    #if LQ, cut in the middle of the window
    my $seq_ref=$_[0];
    my $chr_q_ref=$_[1];
    my @epcs;
	my $qual_chr_offset=$_[2];
    my $epc_max=$_[3];#maximum error probability
	my $epc_win_sz=$_[4];#in this bp window
	my $epc_win_step=int($epc_win_sz/2);
	if($epc_win_step<1){$epc_win_step=1};
	
	unless($epc_max==0){ #if we are exactly at 0 - clipping/filtering off
		foreach my $qual (split ("", $$chr_q_ref)){
			push @epcs,	10**(-(ord($qual)-$qual_chr_offset)/10);
		}
		my $i=0;
		my $epc=0;
		if($epc_win_sz==10){#w/o foreach, speed optimised.
			while($i+$epc_win_sz<=$#epcs and $epc<=$epc_max){
				$epc=$epcs[$i]+$epcs[$i+1]+$epcs[$i+2]+$epcs[$i+3]+$epcs[$i+4]+$epcs[$i+5]+$epcs[$i+6]+$epcs[$i+7]+$epcs[$i+8]+$epcs[$i+9];
				$i+=$epc_win_step;
			}
		}else{#slower version of the sub
			my @win_qual;
			while($i+$epc_win_sz<=$#epcs and $epc<=$epc_max){
				$epc=0;
				foreach my $ep (@epcs[$i, ($i+$epc_win_sz-1)]){
					$epc+=$ep;
				}
				$i+=$epc_win_step;
			}
		}
		if($epc>$epc_max){
			#trim the seq
			$$seq_ref=substr($$seq_ref, 0, ($i+1));
			$$chr_q_ref=substr($$chr_q_ref, 0, ($i+1));
		}
	}
	return(length($$chr_q_ref));
}

sub LQ_GC_clip {# \$seq1t, \$qual1t, \$trim1); #qual offset
	#detect low quality unballanced G/C read
	my $seq1tr=$_[0];
	my $qual1tr=$_[1];
	my $trim1r=$_[2];
	my $q_offset=$_[3];
	my $epc_read=0;
	my $seqb=$$seq1tr;
	my $seq_l=length($seqb);
	unless($seq_l){return(0)};
	my $Gs=($seqb =~ tr/Gg/__/);
	my $Cs=($seqb =~ tr/Cc/__/);
	#foreach my $qual (split ("", $$seq1tr)){
	#	$epc_read+=	10**(-(ord($qual)-$q_offset)/10);
	#}
	my $epc_avg=$epc_read/$seq_l;
	#print STDERR "AVG_EPC: $epc_avg, $Gs, $Cs Read: $$seq1tr\n";
	if(($Gs>($Cs*3) or $Cs>($Gs*2)) and (($Gs+$Cs)/$seq_l)>0.7){
		#print STDERR "Deleting GC unballanced read: $$seq1tr $$qual1tr\n";
		#LQ and unballanced (artifacts) - delete read
		$$trim1r=0;$$seq1tr="";$$qual1tr="";
	}
}

sub make_seq_words {#for making local words database from our subject sequence
	#words score 1 for unicore/complement strand match
	my @seqs= @{$_[0]};
	my $word_size=$_[1];
	my %words_db;
	foreach my $seq (@seqs){
		$seq=uc($seq);#input db sequence
		my $seq_rc=reverse($seq);$seq_rc=~tr/ATCG/TAGC/;
		#my @seq_u=split("",$seq);
		#my @seq_c=split("",$seq_rc);
		my $i; my $word;
		my $seql=length($seq);
		#making forward words list
		
		for($i=0;$i<=($seql-$word_size); $i++){
			$word=substr($seq,$i,$word_size);
			$words_db{$word}=1;
			#		print ("DB word: $word_u\n");
		}
		#making reverse words list
		for($i=0;$i<=($seql-$word_size); $i++){
			$word=substr($seq_rc,$i,$word_size);
			$words_db{$word}=1;
		}
	}
	return(\%words_db)
}


sub fastq2fastaq_ref {#put the quality in the fastaq format: (CHR+64) from fastq: (CHR+33)
    my @chr_q=split("",${$_[0]});
    my @chr_q64;
	foreach my $q (@chr_q){
#		print $q." ".chr((ord($q))+31)."\n";
		push @chr_q64, chr((ord($q))+31);
	}
    ${$_[0]}=join("",@chr_q64);
}
sub get_ARGV2 {	#for getting command line arguments - now supports default presets
	my $argv_defaults=$_[0];
	my %argv_config;
	if($argv_defaults=~m/HASH/){
		foreach my $key (keys %{$argv_defaults}){
				$argv_config{$key}=$argv_defaults->{$key};
		}
	}
	my @files_list;
	foreach my $argument (@ARGV){
		if($argument=~m/^-+(\S+?)=(\S+)/){
			$argv_config{$1}=$2;
		}elsif($argument=~m/^-+(\S+)/){
			$argv_config{$1}=1;
		}else{
			push @files_list, $argument;
		}
	}
	if(@files_list){
		$argv_config{files_list}=\@files_list;
	}
	return(\%argv_config);
}

sub format_number { #some string formatting for nice ordering of filenames (pad the beggining with 0's)
	my $input_number = $_[0];
	my $decimal_places;
	if ($_[1]){
		$decimal_places = $_[1];
	}else{$decimal_places=4}
	while(length($input_number)<$decimal_places){
		$input_number = "0"."$input_number";
	}
	return($input_number);
}

sub split_mates {#separate two halfs of the mate pair read
	#linkerdb, word_sz, min_readl, template
	my $linker_db=$_[0];
	my $word_size=$_[1];
	my $min_rdl=$_[2];#minimal length of the read, should be 25 in production version
	my $paired=$_[3]; #if there are 2 reads, than=1;
	my $template=$_[4];#$read_in->{BEGIN_SEQUENCE};
	my $ends=$_[5]; #the naming schema for the reads ends like:_F or .f or /1
	#now 1-st read data
	my $seq1 = $_[6];
	my $qual1= $_[7]; 
# 	unless($paired){
# 		clip_linker_edge(\$seq1, \$qual1);
# 	}
	my $seql1= length($seq1);
	my $seq2; my $qual2; my $seql2;#pos are 0-based
	my ($lnk1s, $lnk1e, $lnk1scr)=get_linker_match($linker_db, $word_size, $seq1, $qual1);
	my ($lnk2s, $lnk2e, $lnk2scr);
	unless(defined $_[8] and defined $_[9]){$paired=0};#only 1x input
	if($paired){#both reads to do
		$seq2  = $_[8];
		$qual2 = $_[9];
		$seql2= length($seq2);
		($lnk2s, $lnk2e, $lnk2scr)=get_linker_match($linker_db, $word_size, $seq2, $qual2);
	}
	#1 - _F/_R;
	#2 .f/.r	- mira original
	#3 /1 /2	- illumina /newbler/velvet
	#4 1:N:0:1
	my $endF="_F";
	my $endR="_R";
	if($ends==2){
		$endF=".f";$endR=".r";
	}elsif($ends==3){
		$endF="/1";$endR="/2";
	}elsif($ends==4){
		$endF=" 1:0:N:1";$endR=" 2:0:N:1";
	}
	#now we have linker positions,
	#get linker ends (splitting points) for $match_group_selected
	#read seq/qual/names indices: 0-1R,1-1F,2-2R,3-2F
	my @sr_names;my @sr_seq;my @sr_qual;
	#	print ("Read: $read_in->{BEGIN_SEQUENCE},\nbest match is from $match_grp_selected->{begin} to $match_grp_selected->{end}, score $match_grp_selected->{score}, to end: $match_grp_selected->{toend}\n");
	my @split_reads;
	#lets decide whether we are going to split the read onto parts
	#figure out which seq goes where
	unless($paired){#lonely or merged read
		#print "Single read: $template, lnk1e:$lnk1e\n";
		if($lnk1e){
			if($lnk1s>=$min_rdl){#rev read
				$sr_seq[0]=reverse(substr($seq1,0,$lnk1s));
				$sr_seq[0]=~tr/ATCGatcg/TAGCtagc/;
				$sr_qual[0]=reverse(substr($qual1,0,$lnk1s));
				$sr_names[0]=$template.$endR;
				#print "seq_r_in: $seq1\nout: $sr_seq[0]\n";
			}
			if(($seql1-$lnk1e)>=$min_rdl){#fwd read
				$sr_seq[3]=substr($seq1,$lnk1e);
				$sr_qual[3]=substr($qual1,$lnk1e);
				$sr_names[3]=$template.$endF;
				#print "seq_f_out: $seq1\nout: $sr_seq[3]\n";
			}
		}elsif($min_rdl<=$seql1){#no link, but still useable - search for linker remains, if not found - leave it unchanged
			#add checks for partial nextera liker start/end:
			clip_linker_edge(\$seq1, \$qual1);
			#print "Unpaired read, no linker: $template\n";			
			#if th read became too short...
			if($min_rdl<=length($seq1)){
				$sr_seq[3]=$seq1; $sr_qual[3]=$qual1; $sr_names[3]=$template;
			};
		}
#		print "Unpaired reads: ".join(", ", @sr_names)."\n";
	}else{#a sweet couple (but still separate reads)
		if($lnk1e){
			if($lnk1s>=$min_rdl){#rev read
				$sr_seq[0]=reverse(substr($seq1,0,$lnk1s));
				$sr_seq[0]=~tr/ATCGatcg/TAGCtagc/;
				$sr_qual[0]=reverse(substr($qual1,0,$lnk1s));
				$sr_names[0]=$template.$endR;
			}
			if(($seql1-$lnk1e)>=$min_rdl){#fwd read
				$sr_seq[1]=substr($seq1,$lnk1e);
				$sr_qual[1]=substr($qual1,$lnk1e);
				$sr_names[1]=$template.$endF;
			}
		}else{#no linker in this half - just flip it
			clip_linker_edge(\$seq1, \$qual1);
			if(length($seq1)>=$min_rdl){#rev read
				$sr_seq[0]=reverse($seq1);
				$sr_seq[0]=~tr/ATCGatcg/TAGCtagc/;
				$sr_qual[0]=reverse($qual1);
				$sr_names[0]=$template.$endR;
			}
		}
		if($lnk2e){
			if($lnk2s>=$min_rdl){
				$sr_seq[2]=substr($seq2,0,$lnk2s);
				$sr_qual[2]=substr($qual2,0,$lnk2s);
				$sr_names[2]=$template.$endR;
			}
			if(($seql2-$lnk2e)>=$min_rdl){#fwd read
				$sr_seq[3]=reverse(substr($seq2,$lnk2e));
				$sr_seq[3]=~tr/ATCGatcg/TAGCtagc/;
				$sr_qual[3]=reverse(substr($qual2,$lnk2e));
				$sr_names[3]=$template.$endF;
			}
		}else{#no linker in this half - just flip it
			clip_linker_edge(\$seq2, \$qual2);
			if(length($seq2)>=$min_rdl){#rev read
				$sr_seq[3]=reverse($seq2);
				$sr_seq[3]=~tr/ATCGatcg/TAGCtagc/;
				$sr_qual[3]=reverse($qual2);
				$sr_names[3]=$template.$endF;
			}
		}
		#now names - where we add 2's?
		#4 reads - take out internal ones
		if($lnk1e and $lnk2e){ #both reads have linkers! Discard stuff in between.
			$sr_seq[1]=undef;$sr_seq[2]=undef;
			$sr_qual[1]=undef;$sr_qual[2]=undef;
			$sr_names[1]=undef;$sr_names[2]=undef;
		}
# 		elsif($sr_seq[1] and $sr_seq[3]){#set F2 
# 		#actually better to set it to the read end bit, due to LQ issues
# 		#PS: add option to dell them
# 			$sr_names[1].="2";
# 		}elsif($sr_seq[0] and $sr_seq[2]){#set R2
# 			$sr_names[2].="2";
# 		}
#		print "Paired reads: ".join(", ", @sr_names)."\n";
	}
	#check for any signifficant linker remains (chimeric multi linker read), and discard it (can also put it to se)
	my $i;
	for($i=0; $i<4;$i++){
		if($sr_seq[$i]=~m/AGATGTGTATAAGAGACAG/ or $sr_seq[$i]=~m/CTGTCTCTTATACACATCT/){
			$sr_seq[$i]=undef;$sr_qual[$i]=undef;$sr_names[$i]=undef;
		}
	}
	#and now return the results
	return(\@sr_names, \@sr_seq, \@sr_qual);
#	$fastq_stat{sff_out_reads_paired}++;
#	$fastq_stat{sff_out_reads_single}++;
#	$fastq_stat{sff_out_reads_empty}++;
}
sub clip_linker_edge {#if no linker found - try detecting the remains of it and removing near read start/end
	#\@sr_seq, \@sr_qual
	my $seq=$_[0];#ref!
	my $qual=$_[1];#ref!
#	my $length_in=length($$seq);
#	my $i;
#	for($i=0; $i<4;$i++){
		if(defined $$seq){
			if($$seq=~m/^TAAGAGACAG/){
				$$seq=substr($$seq,10); $$qual=substr($$qual,10);
			}elsif($$seq=~m/CTGTCTCTTA$/){
				$$seq=substr($$seq,0,-10); $$qual=substr($$qual,0,-10);
			}elsif($$seq=~m/^AAGAGACAG/){
				$$seq=substr($$seq,9); $$qual=substr($$qual,9);
			}elsif($$seq=~m/CTGTCTCTT$/){
				$$seq=substr($$seq,0,-9); $$qual=substr($$qual,0,-9);
			}elsif($$seq=~m/^AGAGACAG/){
				$$seq=substr($$seq,8); $$qual=substr($$qual,8);
			}elsif($$seq=~m/CTGTCTCT$/){
				$$seq=substr($$seq,0,-8); $$qual=substr($$qual,0,-8);
			}elsif($$seq=~m/^GAGACAG/){
				$$seq=substr($$seq,7); $$qual=substr($$qual,7);
			}elsif($$seq=~m/CTGTCTC$/){
				$$seq=substr($$seq,0,-7); $$qual=substr($$qual,0,-7);
			}elsif($$seq=~m/^AGACAG/){
				$$seq=substr($$seq,6); $$qual=substr($$qual,6);
			}elsif($$seq=~m/CTGTCT$/){
				$$seq=substr($$seq,0,-6); $$qual=substr($$qual,0,-6);
			}elsif($$seq=~m/^GACAG/){
				$$seq=substr($$seq,5); $$qual=substr($$qual,5);
			}elsif($$seq=~m/CTGTC$/){
				$$seq=substr($$seq,0,-5); $$qual=substr($$qual,0,-5);
			}
#			my $length_out=length($$seq);
#			if($length_in != $length_out){
#				print "before clip: $length_in after clip: $length_out $$seq\n";
#			}
		}else{die "Undef linker edge clip input!\n"}
#	}
}

sub get_linker_match {#get linker words and sequence, return linker position.
#$linker_db, $word_size, $seq1, $qual1
	my $linker_db=$_[0];
	my $word_size=$_[1];
	my $seq=$_[2];
	my $seqquc=uc($seq);#uppercqse query sequence
	my $qual=$_[3];
	my $seql=length($seq);
	#my $read_in=$_[];
	my $minscore=4;#min number of words matching in the middle of the read.
	my $minscore_edge=2;#minimal number of words matching near the end.
	my $max_toend=15;#maximal number of bp's before end, to use minscore edge;
	my $maxgap=20;#max gap in linker match (better do some qual chks)
	#adjust maxgap/minscore depending on the quality of the input dataset and trimming
	#first we find sequence matches to our db HQ defaults - maxgap=5, minscore=6
	
	#proc first read
	my $i;
	my $nm=0;
	my $word_query;
	my @m_begin;my @m_end;my @m_score;
	#print "Splitting mate paired reads: $seqquc words size: $word_size\n";
	#print "db: ".join("\n", keys(%{$linker_db})) ."\n";
#	print STDERR "Read: $seqquc\n";
	for($i=0; $i<=($seql-$word_size);$i++){
		$word_query=substr($seqquc,$i,$word_size);
		#join("",@{${$read_in}{seq}}[$i..($i+$argv_config->{word_size}-1)]);
		#		print STDERR "Query word: $word_query\n";
		if($linker_db->{$word_query}){#we have mach to word in the db

			$m_begin[$nm]=$i;
			$m_score[$nm]=0;#increment it in the cycle
			while($linker_db->{$word_query}){
				#print STDERR "Found match to: $word_query\n";
				$m_end[$nm]=$i+$word_size;
				$i++;
				$word_query=substr($seqquc,$i,$word_size);
				$m_score[$nm]++;
			}
			#$match{length}=$match{end}-$match{begin};
			#push @matches, \%match;
			if($argv->{debug}){
				print STDERR "Match from $m_begin[$nm] to $m_end[$nm], score: $m_score[$nm]\n";
			}
			$nm++;
		}
	}
	if($nm==0){
#		die "No matches. debug stop.\n";
		#No matches. - Return original data;
		$fastq_stat{sff_in_reads_wo_linker}++;
		#next; #do next read
		return(undef,undef,undef);
	}else{
#		print "\nFound $nm match(es).\n";
		#group matches
		$i=0;#match index
		my @matches_grp;
		my $max_score=0;
		my $mgi=0;#match group index
		my ($mgt_begin, $mgt_end, $mgt_score, $mgt_toend);#temp match group vars
		my @mg_begin; my @mg_end; my @mg_score;
		while($i<=$#m_begin){
			#my %match_grp;
			$mgt_begin=$m_begin[$i];
			#print "First match from $match_grp{begin} to $matches[$i]->{end}\n";
			do{
				$mgt_end=$m_end[$i];
				#print "Match end: $matches[$i]->{end}\t";
				$mgt_score+=$m_score[$i];
				$i++;
				#print "Working!";
			}while(defined $m_begin[$i] and ($m_begin[$i]-$mgt_end)<=$maxgap);
			#$i--;
			$mgt_toend=$mgt_begin>($seql-$mgt_end)?($seql-$mgt_end):$mgt_begin;
			#before pushing in we do some second level filtering:
			if(($mgt_toend<=$max_toend and $mgt_score>=$minscore_edge) 
				or($mgt_score>=$minscore)){
				$max_score=$max_score<$mgt_score?$mgt_score:$max_score;
				#add group to the candidates:
				$mg_begin[$mgi]=$mgt_begin;$mg_end[$mgi]=$mgt_end;$mg_score[$mgi]=$mgt_score;
				#print "Match group: $mg_begin[$mgi] to $mg_end[$mgi], toend: $mgt_toend, $mg_score[$mgi], max: $max_score\n";
				#print "Match group: $match_grp{begin} to $match_grp{end}, toend: $match_grp{toend}, $match_grp{score}, max: $max_score\n";
				$mgi++;
				$fastq_stat{sff_in_reads_w_linker}++;
			}else{
				$fastq_stat{sff_in_reads_wo_linker}++;
			}
		}
		#return if no significant hits:
		unless(defined $mg_begin[0]){
			#print "Exiting\n";
			return(undef,undef,undef)};
		#select best linker match
		for($mgi=0;$mgi<=$#mg_begin;$mgi++){
			if($mg_score[$mgi] == $max_score){
				$mgt_begin=$mg_begin[$mgi];
				$mgt_end=$mg_end[$mgi];
				#print "Best match: $mg_begin[$mgi]; to $mg_end[$mgi]; $mgi\n";
			}
			#print "Match: $mgt_begin to $mgt_end\n";
		}
		#here comes the best linker match
		return($mgt_begin, $mgt_end, $max_score);
	}
}
sub check_remove_vec {#check wheather (any) reads hit the vector sequence, if yes - remove both
	my $vector_db = $_[0];
	my $kmer = $_[1];
	my $seq1t = $_[2];
	my $qual1t = $_[3];
	my $seq2t = $_[4];
	my $qual2t =$_[5];
	my $is_vector=0;
	my $seed_step=int($kmer+1)/2;
	my $total_seeds=0;
	my $match_seeds=0;
	my $last_seed_pos;
	my $i; my $queryw;
	#Check read1;
	if(ref $seq1t){
		if(length($$seq1t)=>$kmer){#check for matches/count seeds
			$total_seeds+=int((length($$seq1t)-$kmer)/$seed_step);
			$last_seed_pos=length($$seq1t)-$kmer;
			for($i=0; $i<=$last_seed_pos;$i+=$seed_step){
				$queryw=substr($$seq1t,$i,$kmer);
				if($vector_db->{$queryw}){$match_seeds++}
			}
		}
	}
	#Check read2;
	if(ref $seq2t){
		if(length($$seq2t)=>$kmer){#check for matches/count seeds
			$total_seeds+=int((length($$seq2t)-$kmer)/$seed_step);
			$last_seed_pos=length($$seq2t)-$kmer;
			for($i=0; $i<=$last_seed_pos;$i+=$seed_step){
				$queryw=substr($$seq2t,$i,$kmer);
				if($vector_db->{$queryw}){$match_seeds++}
			}
		}
	}
	#make the conclusion - remove if more than 25% of kmers hit the vector
	if($match_seeds and $match_seeds*4>$total_seeds){
		if(ref $seq1t){
			$$seq1t=""; $$qual1t="";
		}
		if(ref $seq2t){
			$$seq2t=""; $$qual2t="";
		}
	}
}


sub make_vec_words {#add an extra kmer-1 start to the end, in order to cover the circular phiX
	#take input seq and chop it onto kmers storing them in a hash
	my $vec_opt=$_[0];#can be file
	my $kmer=$_[1];#wordsize 
	#check for file...(for future versions)
	#default contaminants go bellow (phiX):
	my @seqs=("TGAATAGCAAAGCCTCTACGCGATTTCATAGTGGAGGCCTCCAGCAATCTTGAACACTCATCCTTAATACCTTTCTTTTTGGGGTAATTATACTCATCGC
	GAATATCCTTAAGAGGGCGTTCAGCAGCCAGCTTGCGGCAAAACTGCGTAACCGTCTTCTCGTTCTCTAAAAACCATTTTTCGTCCCCTTCGGGGCGGTG
	GTCTATAGTGTTATTAATATCAAGTTGGGGGAGCACATTGTAGCATTGTGCCAATTCATCCATTAACTTCTCAGTAACAGATACAAACTCATCACGAACG
	TCAGAAGCAGCCTTATGGCCGTCAACATACATATCACCATTATCGAACTCAACGCCCTGCATACGAAAAGACAGAATCTCTTCCAAGAGCTTGATGCGGT
	TATCCATCTGCTTATGGAAGCCAAGCATTGGGGATTGAGAAAGAGTAGAAATGCCACAAGCCTCAATAGCAGGTTTAAGAGCCTCGATACGCTCAAAGTC
	AAAATAATCAGCGTGACATTCAGAAGGGTAATAAGAACGAACCATAAAAAAGCCTCCAAGATTTGGAGGCATGAAAACATACAATTGGGAGGGTGTCAAT
	CCTGACGGTTATTTCCTAGACAAATTAGAGCCAATACCATCAGCTTTACCGTCTTTCCAGAAATTGTTCCAAGTATCGGCAACAGCTTTATCAATACCAT
	GAAAAATATCAACCACACCAGAAGCAGCATCAGTGACGACATTAGAAATATCCTTTGCAGTAGCGCCAATATGAGAAGAGCCATACCGCTGATTCTGCGT
	TTGCTGATGAACTAAGTCAACCTCAGCACTAACCTTGCGAGTCATTTCTTTGATTTGGTCATTGGTAAAATACTGACCAGCCGTTTGAGCTTGAGTAAGC
	ATTTGGCGCATAATCTCGGAAACCTGCTGTTGCTTGGAAAGATTGGTGTTTTCCATAATAGACGCAACGCGAGCAGTAGACTCCTTCTGTTGATAAGCAA
	GCATCTCATTTTGTGCATATACCTGGTCTTTCGTATTCTGGCGTGAAGTCGCCGACTGAATGCCAGCAATCTCTTTTTGAGTCTCATTTTGCATCTCGGC
	AATCTCTTTCTGATTGTCCAGTTGCATTTTAGTAAGCTCTTTTTGATTCTCAAATCCGGCGTCAACCATACCAGCAGAGGAAGCATCAGCACCAGCACGC
	TCCCAAGCATTAAGCTCAGGAAATGCAGCAGCAAGATAATCACGAGTATCCTTTCCTTTATCAGCGGCAGACTTGCCACCAAGTCCAACCAAATCAAGCA
	ACTTATCAGAAACGGCAGAAGTGCCAGCCTGCAACGTACCTTCAAGAAGTCCTTTACCAGCTTTAGCCATAGCACCAGAAACAAAACTAGGGACGGCCTC
	ATCAGGGTTAGGAACATTAGAGCCTTGAATGGCAGATTTAATACCAGCATCACCCATGCCTACAGTATTGTTATCGGTAGCAAGCACATCACCTTGAATG
	CCACCGGAGGCGGCTTTTTGACCGCCTCCAAACAATTTAGACATGGCGCCACCAGCAAGAGCAGAAGCAATACCGCCAGCAATAGCACCAAACATAAATC
	ACCTCACTTAAGTGGCTGGAGACAAATAATCTCTTTAATAACCTGATTCAGCGAAACCAATCCGCGGCATTTAGTAGCGGTAAAGTTAGACCAAACCATG
	AAACCAACATAAACATTATTGCCCGGCGTACGGGGAAGGACGTCAATAGTCACACAGTCCTTGACGGTATAATAACCACCATCATGGCGACCATCCAAAG
	GATAAACATCATAGGCAGTCGGGAGGGTAGTCGGAACCGAAGAAGACTCAAAGCGAACCAAACAGGCAAAAAATTTAGGGTCGGCATCAAAAGCAATATC
	AGCACCAACAGAAACAACCTGATTAGCGGCGTTGACAGATGTATCCATCTGAATGCAATGAAGAAAACCACCATTACCAGCATTAACCGTCAAACTATCA
	AAATATAACGTTGACGATGTAGCTTTAGGTGTCTGTAAAACAGGTGCCGAAGAAGCTGGAGTAACAGAAGTGAGAACCAGCTTATCAGAAAAAAAGTTTG
	AATTATGGCGAGAAATAAAAGTCTGAAACATGATTAAACTCCTAAGCAGAAAACCTACCGCGCTTCGCTTGGTCAACCCCTCAGCGGCAAAAATTAAAAT
	TTTTACCGCTTCGGCGTTATAACCTCACACTCAATCTTTTATCACGAAGTCATGATTGAATCGCGAGTGGTCGGCAGATTGCGATAAACGGTCACATTAA
	ATTTAACCTGACTATTCCACTGCAACAACTGAACGGACTGGAAACACTGGTCATAATCATGGTGGCGAATAAGTACGCGTTCTTGCAAATCACCAGAAGG
	CGGTTCCTGAATGAATGGGAAGCCTTCAAGAAGGTGATAAGCAGGAGAAACATACGAAGGCGCATAACGATACCACTGACCCTCAGCAATCTTAAACTTC
	TTAGACGAATCACCAGAACGGAAAACATCCTTCATAGAAATTTCACGCGGCGGCAAGTTGCCATACAAAACAGGGTCGCCAGCAATATCGGTATAAGTCA
	AAGCACCTTTAGCGTTAAGGTACTGAATCTCTTTAGTCGCAGTAGGCGGAAAACGAACAAGCGCAAGAGTAAACATAGTGCCATGCTCAGGAACAAAGAA
	ACGCGGCACAGAATGTTTATAGGTCTGTTGAACACGACCAGAAAACTGGCCTAACGACGTTTGGTCAGTTCCATCAACATCATAGCCAGATGCCCAGAGA
	TTAGAGCGCATGACAAGTAAAGGACGGTTGTCAGCGTCATAAGAGGTTTTACCTCCAAATGAAGAAATAACATCATGGTAACGCTGCATGAAGTAATCAC
	GTTCTTGGTCAGTATGCAAATTAGCATAAGCAGCTTGCAGACCCATAATGTCAATAGATGTGGTAGAAGTCGTCATTTGGCGAGAAAGCTCAGTCTCAGG
	AGGAAGCGGAGCAGTCCAAATGTTTTTGAGATGGCAGCAACGGAAACCATAACGAGCATCATCTTGATTAAGCTCATTAGGGTTAGCCTCGGTACGGTCA
	GGCATCCACGGCGCTTTAAAATAGTTGTTATAGATATTCAAATAACCCTGAAACAAATGCTTAGGGATTTTATTGGTATCAGGGTTAATCGTGCCAAGAA
	AAGCGGCATGGTCAATATAACCAGTAGTGTTAACAGTCGGGAGAGGAGTGGCATTAACACCATCCTTCATGAACTTAATCCACTGTTCACCATAAACGTG
	ACGATGAGGGACATAAAAAGTAAAAATGTCTACAGTAGAGTCAATAGCAAGGCCACGACGCAATGGAGAAAGACGGAGAGCGCCAACGGCGTCCATCTCG
	AAGGAGTCGCCAGCGATAACCGGAGTAGTTGAAATGGTAATAAGACGACCAATCTGACCAGCAAGGAAGCCAAGATGGGAAAGGTCATGCGGCATACGCT
	CGGCGCCAGTTTGAATATTAGACATAATTTATCCTCAAGTAAGGGGCCGAAGCCCCTGCAATTAAAATTGTTGACCACCTACATACCAAAGACGAGCGCC
	TTTACGCTTGCCTTTAGTACCTCGCAACGGCTGCGGACGACCAGGGCGAGCGCCAGAACGTTTTTTACCTTTAGACATTACATCACTCCTTCTGCACGTA
	ATTTTTGACGCACGTTTTCTTCTGCGTCAGTAAGAACGTCAGTGTTTCCTGCGCGTACACGCAAGGTAAACGCGAACAATTCAGCGGCTTTAACCGGACG
	CTCGACGCCATTAATAATGTTTTCCGTAAATTCAGCGCCTTCCATGATGAGACAGGCCGTTTGAATGTTGACGGGATGAACATAATAAGCAATGACGGCA
	GCAATAAACTCAACAGGAGCAGGAAAGCGAGGGTATCCTACAAAGTCCAGCGTACCATAAACGCAAGCCTCAACGCAGCGACGAGCACGAGAGCGGTCAG
	TAGCAATCCAAACTTTGTTACTCGTCAGAAAATCGAAATCATCTTCGGTTAAATCCAAAACGGCAGAAGCCTGAATGAGCTTAATAGAGGCCAAAGCGGT
	CTGGAAACGTACGGATTGTTCAGTAACTTGACTCATGATTTCTTACCTATTAGTGGTTGAACAGCATCGGACTCAGATAGTAATCCACGCTCTTTTAAAA
	TGTCAACAAGAGAATCTCTACCATGAACAAAATGTGACTCATATCTAAACCAGTCCTTGACGAACGTGCCAAGCATATTAAGCCACTTCTCCTCATCCAA
	CGCGTCAGTTTTTGACAGAATCGTTAGTTGATGGCGAAAGGTCGCAAAGTAAGAGCTTCTCGAGCTGCGCAAGGATAGGTCGAATTTTCTCATTTTCCGC
	CAGCAGTCCACTTCGATTTAATTCGTAAACAAGCAGTAGTAATTCCTGCTTTATCAAGATAATTTTTCGACTCATCAGAAATATCCGAAAGTGTTAACTT
	CTGCGTCATGGAAGCGATAAAACTCTGCAGGTTGGATACGCCAATCATTTTTATCGAAGCGCGCATAAATTTGAGCAGATTTGTCGTCACAGGTTGCGCC
	GCCAAAACGTCGGCTACAGTAACTTTTCCCAGCCTCAATCTCATCTCTCTTTTTGCGTTCTGCTTCAATATCTGGTTGAACGGCGTCGCGTCGTAACCCA
	GCTTGGTAAGTTGGATTAAGCACTCCGTGGACAGATTTGTCATTGTGAGCATTTTCATCCCGAAGTTGCGGCTCATTCTGATTCTGAACAGCTTCTTGGG
	AAGTAGCGACAGCTTGGTTTTTAGTGAGTTGTTCCATTCTTTAGCTCCTAGACCTTTAGCAGCAAGGTCCATATCTGACTTTTTGTTAACGTATTTAGCC
	ACATAGAAACCAACAGCCATATAACTGGTAGCTTTAAGCGGCTCACCTTTAGCATCAACAGGCCACAACCAACCAGAACGTGAAAAAGCGTCCTGCGTGT
	AGCGAACTGCGATGGGCATACTGTAACCATAAGGCCACGTATTTTGCAAGCTATTTAACTGGCGGCGATTGCGTACCCGACGACCAAAATTAGGGTCAAC
	GCTACCTGTAGGAAGTGTCCGCATAAAGTGCACCGCATGGAAATGAAGACGGCCATTAGCTGTACCATACTCAGGCACACAAAAATACTGATAGCAGTCG
	GCGTGTGAATCATTAGCCTTGCGACCCTCGGCAGCAAGAACCATACGACCAATATCACGAAAATAGTCACGCAAAGCATTGGGATTATCATAAAACGCCT
	CTAATCGGTCGTCAGCCAACGTGAGAGTGTCAAAAACGATAAACCAACCATCAGCATGAGCCTGTCGCATTGCATTCATCAAACGC");
	if(-s $vec_opt){
		print STDERR "Trying to read the contaminant sequences database from the fasta file: $vec_opt\n";
		open(SEQ_FAS, $vec_opt)or die "Unable to read the input database: $vec_opt\n";
		@seqs=();#empty default seqs
		my $seq_fas="";
		while(<SEQ_FAS>){
			s/\r//g;
			chomp;
			if(m/^>/){
				if($seq_fas=~m/\S{$kmer}/){#something usefull in there...
					push @seqs, $seq_fas;
				}
				$seq_fas="";
			}else{
				$seq_fas.=$_;
			}
		}
		if($seq_fas=~m/\S{$kmer}/){
			push @seqs, $seq_fas;
		}
		close(SEQ_FAS);
	}else{
		print STDERR "Filtering control reads using the default phiX sequence.\n";
	}
	
	#process circular contaminants
	foreach my $seq (@seqs){
		$seq=~s/\s//gs;
		$seq=~s/\t//gs;
		$seq=~s/\n//gs;
		$seq.=substr($seq,0,($kmer-1));
	}
	$vector_db = make_seq_words(\@seqs, $kmer);
	return($vector_db);
}