2bRADExtraction.pl

#!/usr/bin/env perl
# Authors: Zheng Sun, Rongchao Zhang, Shi Huang
# Last update: 2020.06.03
use strict;
use warnings;
use Getopt::Long;
use FindBin qw($Bin);
use File::Basename qw(dirname basename);
no strict 'refs';

my $author="Zheng Sun, Rongchao Zhang, Shi Huang";
my $time="2020.06.03";

#set default parameters
my $ncount ||=0.08;
my $quality ||=30;
my $percent ||=80;
my $qbase ||=33;
my $format ||="fa";
my $gz ||="yes";
my $q_control ||="yes";
my $pear_cpu ||=1;
#set software path
my $pear ||="pear";

select STDOUT;$|=1;# Standard output for clearing cache

my (@input,$type,$site,$outdir,@outprefix);
my $help;
GetOptions(
	"i:s{1,2}"  => \@input,# single-end or double-end reads
	"t:i"  => \$type, #fa：reference genome data；fq：shotgun metagenomics data，single 2b-RAD tags，five concatenated 2b-RAD tags
	"s:i"  => \$site, # restriction sites
	"od:s" => \$outdir, # output directory
	"op:s{1,5}" => \@outprefix, # output prefix

	"gz:s" => \$gz, # do compression or not for the outputs

	"qc:s" => \$q_control, # QC: Whether QC is required
	"n:f"  => \$ncount, # QC: The # of ambiguity nucleotides N allowed
	"q:i"  => \$quality, # QC: Minimum per-base quality score (XX)
	"p:i"  => \$percent, # Minimum percentage of bases that must have per-base quality score over [-q]
	"b:i"  => \$qbase, # Quality values of bases
	"fm:s" => \$format, # The output file format: fa/fq

	#software
	"pe:s" => \$pear,#pear soft
	"pc:i" => \$pear_cpu,#cpu of pear soft

	"h|help:s" => \$help,
	);

sub usage{
	print STDERR "\e[;33;1m
Description
    It performs digital type-2B-restriction disgest of DNA data generated by a wide range of sequencing protocols by one of 16 restriction enzymes.For a given type 2b restriction enzyme, it can return a Fasta file including resulting 2b-RAD tags, and a statistical summary including raw number ofinput sequences, restriction enzyme used, number of restriction fragments produced, percentage of restriction fragments over the whole (meta)genomedata. The four application scenarios of this script are listed as below:
    1.To extract 2b-RAD tags from reference genome(s) data, run: perl EeTt.pl -i genome.fa(.gz) -t 1 -s 1 -od . -op sample
    2.To extract 2b-RAD tags from shotgun metagenomic sequencing data, run: perl EeTt.pl -i shotgun.R1.fq(.gz) (or shotgun.R1.fq.gz  shotgun.R2.fq.gz) -t 2 -s 1 -od . -op sample
    3.To extract single 2b-RAD tags from SE or PE sequencing data, run: perl EeTt.pl -i 2b-RADsingle.fq(.gz) (or 2b-RADsingle.R1.fq.gz from PE platform) -t 3 -s 1 -od . -op sample
    4.To split five concatenated 2b-RAD tags from PE sequencing data, run: perl EeTt.pl -i 2b-RAD.R1.fq(.gz) 2b-RAD.R2.fq(.gz) -t 4 -s 1 -od . -op sample1 sample2 sample3 sample4 sample5
Usage
    perl $0 -i <input_file> -t <type> -s <site> -od <outdir> -op <outprefix> [options]*
Required:
    -i  <str>    Input fasta/fastq filepath (.gz supported)
    -t  <int>    The type of input fasta/fastq file
                 [1] Reference genome data in the Fasta format
                 [2] Shotgun (meta)genome data in the Fastq format
                 [3] SE/PE sequencing data in the Fastq format
                 [4] PE sequencing data in the Fastq format
    -s  <int>    One of the type 2b restriction enzymes (sites).
                 [1]CspCI  [9]BplI
                 [2]AloI   [10]FalI
                 [3]BsaXI  [11]Bsp24I
                 [4]BaeI   [12]HaeIV
                 [5]BcgI   [13]CjePI
                 [6]CjeI   [14]Hin4I
                 [7]PpiI   [15]AlfI
                 [8]PsrI   [16]BslFI
    -od <str>    The output directory (automatically create if it does not exist)
    -op <str>    The output prefix (recommended: sample name(s))
Optional:
    -gz <str>    Whether the output file is compressed (yes or no) [$gz]
    -h|help      print this help
Optional (only applicable when -t equals 2, or 3, or 4, i.e. taking fastq data as input):
    -qc <str>    Whether quality control is required (yes or no) [$q_control]
    -n  <float>  Maximum percentage of ambiguity bases \"N\" [default: $ncount]
    -q  <int>    Minimum per-base quality score [default: $quality]
    -p  <int>    Minimum percentage of bases that must have per-base quality score over [-q] [default: $percent]
    -b  <int>    Phred quality score type [default: $qbase]
    -fm <str>    Output file format (fa or fq) [default: $format]
    -pe  <str>   Path of pear soft [$pear]
    -pc  <int>   Cpu of pear soft [$pear_cpu]
Author  $author $time\e[0m\n";
}

if(defined($help)){
	&usage;
	exit 0;
}

unless(@input && $type && $site && $outdir && @outprefix){# parameters checking
	&usage;
	print STDERR "Please check parameter -i -t -s -od -op\n";
	exit 1;
}



# checking input args
unless($gz eq "yes" || $gz eq "no"){
	&usage;
	print STDERR "Parameter -gz is wrong\n";
	exit 1;
}
unless($format eq "fa" || $format eq "fq"){
	&usage;
	print STDERR "Parameter -fm is wrong\n";
	exit 1;
}
unless($q_control eq "yes" || $q_control eq "no"){
	&usage;
	print STDERR "Parameter -qc is wrong\n";
	exit 1;
}

# Define the DNA sequences at a given restriction enzyme site
my ($enzyme,@site,@start,@end,$minpear,$maxpear);
if( 1 == $site ){#CspCI
    @site = (
            '[AGCT]{11}CAA[AGCT]{5}GTGG[AGCT]{10}',
            '[AGCT]{10}CCAC[AGCT]{5}TTG[AGCT]{11}',
            );
	$enzyme="CspCI";
	@start = (0,37,78,119,160);
	@end   = (41,82,123,164,205);
}elsif( 2 == $site ){#AloI
    @site = (
            '[AGCT]{7}GAAC[AGCT]{6}TCC[AGCT]{7}',
            '[AGCT]{7}GGA[AGCT]{6}GTTC[AGCT]{7}',
            );
	$enzyme="AloI";
	@start = (0,38,80,122,164);
	@end   = (42,84,126,168,210);
}elsif( 3 == $site ){#BsaXI
    @site = (
            '[AGCT]{9}AC[AGCT]{5}CTCC[AGCT]{7}',
            '[AGCT]{7}GGAG[AGCT]{5}GT[AGCT]{9}',
            );
	$enzyme="BsaXI";
	@start = (0,33,69,105,141);
	@end   = (35,71,107,143,180);
	$minpear ||=173;
	$maxpear ||=181;
}elsif( 4 == $site ){#BaeI
    @site = (
            '[AGCT]{10}AC[AGCT]{4}GTA[CT]C[AGCT]{7}',
            '[AGCT]{7}G[AG]TAC[AGCT]{4}GT[AGCT]{10}',
            );
	$enzyme="BaeI";
	@start = (0,38,79,120,161);
	@end   = (40,81,122,163,205);
	$minpear ||=198;
	$maxpear ||=206;
}elsif( 5 == $site ){#BcgI
    @site = (
            '[AGCT]{10}CGA[AGCT]{6}TGC[AGCT]{10}',
            '[AGCT]{10}GCA[AGCT]{6}TCG[AGCT]{10}',
            );
	$enzyme="BcgI";
	@start = (0,36,75,114,153);
	@end   = (38,77,116,155,195);
	$minpear ||=188;
	$maxpear ||=196;
}elsif( 6 == $site ){#CjeI
    @site = (
            '[AGCT]{8}CCA[AGCT]{6}GT[AGCT]{9}',
            '[AGCT]{9}AC[AGCT]{6}TGG[AGCT]{8}',
            );
	$enzyme="CjeI";
	@start = (0,40,83,126,169);
	@end   = (42,85,128,171,214);
}elsif( 7 == $site ){#PpiI
    @site = (
            '[AGCT]{7}GAAC[AGCT]{5}CTC[AGCT]{8}',
            '[AGCT]{8}GAG[AGCT]{5}GTTC[AGCT]{7}',
            );
	$enzyme="PpiI";
	@start = (0,37,77,117,157);
	@end   = (39,79,119,159,199);
}elsif( 8 == $site ){#PsrI
    @site = (
            '[AGCT]{7}GAAC[AGCT]{6}TAC[AGCT]{7}',
            '[AGCT]{7}GTA[AGCT]{6}GTTC[AGCT]{7}',
            );
	$enzyme="PsrI"; 
	@start = (0,37,77,117,157);
	@end   = (39,79,119,159,199);
}elsif( 9 == $site ){#BplI
    @site = (
            '[AGCT]{8}GAG[AGCT]{5}CTC[AGCT]{8}', #palindromes
            );
	$enzyme="BplI";
	@start = (0,37,77,117,157);
	@end   = (39,79,119,159,199);
}elsif( 10 == $site ){#FalI
    @site = (
            '[AGCT]{8}AAG[AGCT]{5}CTT[AGCT]{8}', #palindromes
            );
	$enzyme="FalI";
	@start = (0,37,77,117,157);
	@end   = (39,79,119,159,200);
	$minpear ||=193;
	$maxpear ||=201;
}elsif( 11 == $site ){#Bsp24I
    @site = (
            '[AGCT]{8}GAC[AGCT]{6}TGG[AGCT]{7}',
            '[AGCT]{7}CCA[AGCT]{6}GTC[AGCT]{8}',
            );
	$enzyme="Bsp24I";
	@start = (0,37,77,117,157);
	@end   = (39,79,119,159,200);
}elsif( 12 == $site ){#HaeIV
    @site = (
            '[AGCT]{7}GA[CT][AGCT]{5}[AG]TC[AGCT]{9}',
            '[AGCT]{9}GA[CT][AGCT]{5}[AG]TC[AGCT]{7}',
            );
	$enzyme="HaeIV";
	@start = (0,38,79,120,161);
	@end   = (40,81,122,163,204);
}elsif( 13 == $site ){#CjePI
    @site = (
            '[AGCT]{7}CCA[AGCT]{7}TC[AGCT]{8}',
            '[AGCT]{8}GA[AGCT]{7}TGG[AGCT]{7}',
            );
	$enzyme="CjePI";
	@start = (0,39,81,123,165);
	@end   = (41,83,125,167,209);
}elsif( 14 == $site ){#Hin4I
    @site = (
            '[AGCT]{8}GA[CT][AGCT]{5}[GAC]TC[AGCT]{8}',
            '[AGCT]{8}GA[CTG][AGCT]{5}[AG]TC[AGCT]{8}',
            );
	$enzyme="Hin4I";
	@start = (0,37,77,117,157);
	@end   = (39,79,119,159,199);
}elsif( 15 == $site ){#AlfI
    @site = (
            '[AGCT]{10}GCA[AGCT]{6}TGC[AGCT]{10}', #palindromes
            );
	$enzyme="AlfI";
	@start = (0,36,75,114,153);
	@end   = (38,77,116,155,194);
}elsif( 16 == $site ){#BslFI ??some question?? single enzyme
    @site = (
            '[AGCT]{6}GGGAC[AGCT]{14}',
            '[AGCT]{14}GTCCC[AGCT]{6}',
            );
	$enzyme="BslFI";
	@start = (0,34,72,110,148);
	@end   = (38,76,114,152,190);
}else{
	&usage;
	print STDERR "The parameter -s is wrong\n";
	exit 1;
}

&CheckDir($outdir);
my $raw_reads_num=0;# check # of raw reads 如果涉及到拼接，shotgun双端 需要单独记录拼接前的reads数目

if($#input==0 && $type==1 && $#outprefix==0){# reference genome(s)
	print STDOUT "COMMAND: perl $0 -i $input[0] -t 1 -s $site -od $outdir -op $outprefix[0] -gz $gz\n";
	print STDOUT "Electronic enzyme digestion of input genome(s) -- Start, ",`date`;
	&Electronic_enzyme;
	print STDOUT "Electronic enzyme digestion of input genome(s) -- End, ",`date`;
}elsif($#input==0 && $type==2 && $#outprefix==0){# single-end reads from shotgun metagenomics
	if($q_control eq "yes"){# need QC
		print STDOUT "COMMAND: perl $0 -i $input[0] -t 2 -s $site -od $outdir -op $outprefix[0] -gz $gz -qc $q_control -n $ncount -q $quality -p $percent -b $qbase -fm $format\n";
	}else{#no QC
		print STDOUT "COMMAND: perl $0 -i $input[0] -t 2 -s $site -od $outdir -op $outprefix[0] -gz $gz -qc $q_control -fm $format\n";
	}
	print STDOUT "Tags extraction from shotgun metagenomics data -- Start, ",`date`;
	&fastq;
	print STDOUT "Tags extraction from shotgun metagenomics data -- End, ",`date`;
}elsif($#input==1 && $type==2 && $#outprefix==0){# paired-end reads from shotgun metagenomics
#	unless(-e "$pear"){
#		&usage;
#		print STDERR "Can not find software $pear\n";
#		exit 1;
#	}
	if($q_control eq "yes"){# need QC
		print STDOUT "COMMAND: perl $0 -i $input[0] $input[1] -t 2 -s $site -od $outdir -op $outprefix[0] -gz $gz -qc $q_control -n $ncount -q $quality -p $percent -b $qbase -fm $format -pe $pear -pc $pear_cpu\n";
	}else{# no QC
		print STDOUT "COMMAND: perl $0 -i $input[0] $input[1] -t 2 -s $site -od $outdir -op $outprefix[0] -gz $gz -qc $q_control -fm $format -pe $pear -pc $pear_cpu\n";
	}
	print STDOUT "Tags extraction from shotgun metagenomics data -- Start, ",`date`;
	&pear; # Merge the paired-end reads using PEAR. PEAR is a very fast and accurate software tool to merge paired-end reads from next-generation sequencing experiments.
	$input[0]="$outdir/$outprefix[0].$enzyme.pear.fastq.gz";
	&fastq;
	&execute("rm -f $outdir/$outprefix[0].$enzyme.pear.fastq.gz");
	print STDOUT "Tags extraction from shotgun metagenomics data -- End, ",`date`;
}elsif($#input==0 && $type==3 && $#outprefix==0){# single iso-RAD tags
	if($q_control eq "yes"){# need QC
		print STDOUT "COMMAND: perl $0 -i $input[0] -t 3 -s $site -od $outdir -op $outprefix[0] -gz $gz -qc $q_control -n $ncount -q $quality -p $percent -b $qbase -fm $format\n";
	}else{#no QC
		print STDOUT "COMMAND: perl $0 -i $input[0] -t 3 -s $site -od $outdir -op $outprefix[0] -gz $gz -qc $q_control -fm $format\n";
	}
	print STDOUT "Data Split for single isoRAD tags from SE Platform  -- Start, ",`date`;
	&Single_Lable;
	print STDOUT "Data Split for Single isoRAD tags from SE Platform -- End, ",`date`;
}elsif($#input==1 && $type==4 && $#outprefix==4){# five concatenated isoRAD tags
#	unless(-e "$pear"){
#		&usage;
#		print STDERR "Can not find $pear\n";
#		exit 1;
#	}
	if($q_control eq "yes"){# need QC
		print STDOUT "COMMAND: perl $0 -i $input[0] $input[1] -t 4 -s $site -od $outdir -op ",join(" ",@outprefix[0..4])," -gz $gz -qc $q_control -n $ncount -q $quality -p $percent -b $qbase -fm $format -pe $pear -pc $pear_cpu\n";
	}else{#不需要质控
		print STDOUT "COMMAND: perl $0 -i $input[0] $input[1] -t 4 -s $site -od $outdir -op ",join(" ",@outprefix[0..4])," -gz $gz -qc $q_control -fm $format -pe $pear -pc $pear_cpu\n";
	}
	print STDOUT "Data split for five concatenated isoRAD tags from PE platform   -- Start, ",`date`;
	&Five_Lable;
	print STDOUT "Data split for five concatenated isoRAD tags from PE platform  -- End, ",`date`;
}else{
	&usage;
	print STDERR "Please check parameter -i -t -op\n";
	exit 1;
}


sub pear{
	my $r1=$input[0];
	my $r2=$input[1];
	if($r1=~/\.gz$/){# count the # of raw reads
		open R,"gzip -dc $r1|" or die "cannot open $r1\n";
	}else{
		open R,"$r1" or die "cannot open $r1\n";
	}
	while(<R>){#统计双端shotgun数据reads数
		$raw_reads_num++;
		<R>;
		<R>;
		<R>;
	}
	close R;
	my $outprefix="$outprefix[0].$enzyme";
	&execute("$pear -f $r1 -r $r2 -e -o $outdir/$outprefix -j $pear_cpu");#pear 拼接，当插入片段过短时，会无法拼接
	&execute("cat $outdir/$outprefix.assembled.fastq $outdir/$outprefix.unassembled.forward.fastq $outdir/$outprefix.unassembled.reverse.fastq | gzip > $outdir/$outprefix.pear.fastq.gz");#拼接后的reads和不能拼接的R1R2合并
	&execute("rm -f $outdir/$outprefix.assembled.fastq $outdir/$outprefix.unassembled.forward.fastq $outdir/$outprefix.unassembled.reverse.fastq");
	&execute("rm -f $outdir/$outprefix.discarded.fastq");
}
sub fastq{
	my $fastq=$input[0];
	my $outprefix=$outprefix[0];
	if($fastq=~/\.gz$/){
		open IN,"gzip -dc $fastq|" or die "cannot open $fastq\n";
	}else{
		open IN,"$fastq" or die "cannot open $fastq\n";
	}
	if($gz eq "yes"){
		open OU,"|gzip >$outdir/$outprefix.$enzyme.$format.gz" or die "cannot open $outdir/$outprefix.$enzyme.$format.gz\n";
	}else{
		open OU,">$outdir/$outprefix.$enzyme.$format" or die "cannot open $outdir/$outprefix.$enzyme.$format\n";
	}
	open STAT,">$outdir/$outprefix.$enzyme.stat.xls" or die "cannot open $outdir/$outprefix.$enzyme.stat.xls\n";
	if($raw_reads_num==0){# single-end reads from shotgun metagenomics
		print STAT "sample\tenzyme\tinput_reads_num\tenzyme_reads_num\tpercent\n";
	}else{#双端shotgun数据
		print STAT "sample\tenzyme\tinput_reads_num\tcombine_uncombineR1R2_reads_num\tenzyme_reads_num\tpercent\n";
	}
	my($input_reads_num,$enzyme_reads_num,$percent_sub);
	$enzyme_reads_num=0;
	while(<IN>){
		$input_reads_num++;# When inputting paired-end shotgun metagenomics reads，we only record the # of merged reads 当输入数据是双端shotgun数据时，此时记录的是拼接后的reads数以及未能拼接的R1R2之和
		my $line=$_ . <IN> . <IN> . <IN>;
		if($q_control eq "yes"){# QC
			next unless(&CheckN($line));
			next unless(&CheckQ($line));
		}
		my @tmp=split /\n/,$line;
		my %uniq;# deredundancy as some restriction sites will have palindromes
		for my $site(@site){
			while($tmp[1]=~/($site)/g){
				my $pos=pos($tmp[1]);
				my $seq=$1;
				my $len=length($seq);
				$pos=$pos-$len+1;
				my $qual=substr($tmp[3],$pos-1,$len);
				pos($tmp[1])=$pos;
				$uniq{$pos}{$len}="$tmp[0]-$pos\n$seq\n+\n$qual\n";
			}
		}
		for my $pos(sort {$a<=>$b} keys %uniq){
			for my $len(sort {$a<=>$b} keys %{$uniq{$pos}}){
				$enzyme_reads_num++;
				my @a=split /\n/,$uniq{$pos}{$len};
				if($format eq "fq"){
					print OU "$a[0]\n$a[1]\n$a[2]\n$a[3]\n";
				}elsif($format eq "fa"){
					$a[0]=~s/^@/>/;
					print OU "$a[0]\n$a[1]\n";
				}
			}
		}
		undef %uniq;
	}
	close IN;
	close OU;
	if($raw_reads_num==0){#single-end reads from shotgun metagenomics
		$percent_sub=sprintf "%.2f",$enzyme_reads_num/$input_reads_num*100;
		print STAT "$outprefix\t$enzyme\t$input_reads_num\t$enzyme_reads_num\t$percent_sub%\n";
	}else{#双端shotgun数据
		$percent_sub=sprintf "%.2f",$enzyme_reads_num/$raw_reads_num*100;
		print STAT "$outprefix\t$enzyme\t$raw_reads_num\t$input_reads_num\t$enzyme_reads_num\t$percent_sub%\n";
	}
	close STAT;
}



sub Five_Lable{# Data split for five concanated 2b-RAD tags
	my $r1=$input[0];
	my $r2=$input[1];
	my ($output,@fhandle);
	my $outprefix=$outprefix[0];# rename intermediate files using the first sample name
	my $input_reads_num;
	if($r1=~/\.gz$/){# record the # of raw reads
		open R,"gzip -dc $r1|" or die "cannot open $r1\n";
	}else{
		open R,"$r1" or die "cannot open $r1\n";
	}
	while(<R>){
		$input_reads_num++;
		<R>;
		<R>;
		<R>;
	}
	close R;
	&execute("$pear -f $r1 -r $r2 -e -n $minpear -m $maxpear -o $outdir/$outprefix -j $pear_cpu");# merge data
	open IN,"$outdir/$outprefix.assembled.fastq" or die "cannot open $outdir/$outprefix.assembled.fastq\n";
	for my $i(1..$#start+1){# open file handle
		my $fh="OU" . $i;
		my $j=$i-1;
		$output="$outdir/$outprefix[$j].$enzyme.$format";
		if($gz eq "yes"){
			open $fh,"|gzip > $output.gz" or die "cannot open $output.gz\n";
		}elsif($gz eq "no"){
			open $fh,"> $output" or die "cannot open $output\n";
		}
		push @fhandle,$fh;
	}
	my $stat_name=join("-",@outprefix[0..4]);
	open STAT,">$outdir/$stat_name.$enzyme.stat.xls" or die "cannot open $outdir/$stat_name.$enzyme.stat.xls\n";
	print STAT "sample\tenzyme\tinput_reads_num\tcombine_reads_num\tenzyme_reads_num\tqc_reads_num\tpercent\n";
	my($combine_reads_num,%enzyme_reads_num,%qc_reads_num);
	$combine_reads_num=0;
	for my $i(0..$#start){
		$qc_reads_num{$i}=0;
		$enzyme_reads_num{$i}=0;
	}
	while(<IN>){
		$combine_reads_num++;#拼接后reads统计
		my $line=$_ . <IN> . <IN> .<IN>;
		my @tmp=split /\n/,$line;
		for my $i(0..$#start){
			my $fh=$fhandle[$i];
			my $num=$i+1;
			my $id="$tmp[0]:$num";#1
			my $seq=substr($tmp[1],$start[$i],$end[$i]-$start[$i]+1);#2
			my $qual=substr($tmp[3],$start[$i],$end[$i]-$start[$i]+1);#4
			for my $j(0..$#site){
				if($seq=~s/^(\S*?)($site[$j])\S*/$2/){
					$enzyme_reads_num{$i}++;#含酶切位点reads
					my $begin=length($1);
					my $len=length($2);
					$qual=substr($qual,$begin,$len);
					my $sub_line="$id\n$seq\n$tmp[2]\n$qual";
					if($q_control eq "yes"){# check if QC needed
						next unless(&CheckN($sub_line));
						next unless(&CheckQ($sub_line));
					}
					$qc_reads_num{$i}++;#质控后reads
					if($format eq "fa"){
						$id=~s/^@/>/;
						print $fh "$id\n$seq\n";
					}elsif($format eq "fq"){
						print $fh "$sub_line\n";
					}
					last;
				}
			}
		}
	}
	close IN;
	for(@fhandle){
		close $_;
	}
	for my $i(0..$#start){
		my $percent_sub=sprintf "%.2f",$qc_reads_num{$i}/$input_reads_num*100;
		print STAT "$outprefix[$i]\t$enzyme\t$input_reads_num\t$combine_reads_num\t$enzyme_reads_num{$i}\t$qc_reads_num{$i}\t$percent_sub%\n";
	}
	close STAT;
	&execute("rm -f $outdir/$outprefix.assembled.fastq");
	&execute("rm -f $outdir/$outprefix.unassembled.forward.fastq $outdir/$outprefix.unassembled.reverse.fastq");
	&execute("rm -f $outdir/$outprefix.discarded.fastq");
}


sub Single_Lable{
	my $r=$input[0];
	my $outprefix=$outprefix[0];
	my $output;
	if($r=~/\.gz$/){
		open IN,"gzip -dc $r|" or die "cannot open $r\n";
	}else{
		open IN,"$r" or die "cannot open $r\n";
	}
	$output="$outdir/$outprefix.$enzyme.$format";
	if($gz eq "yes"){
		open OU,"|gzip > $output.gz" or die "cannot open $output.gz\n";
	}elsif($gz eq "no"){
		open OU,"> $output" or die "cannot open $output\n";
	}
	open STAT,">$outdir/$outprefix.$enzyme.stat.xls" or die "cannot open $outdir/$outprefix.$enzyme.stat.xls\n";
	print STAT "sample\tenzyme\tinput_reads_num\tenzyme_reads_num\tqc_reads_num\tpercent\n";
	my ($input_reads_num,$enzyme_reads_num,$qc_reads_num,$percent_sub);
	$qc_reads_num=0;
	while(<IN>){
		$input_reads_num++;#原始数据reads
		my $line=$_ . <IN> . <IN> . <IN>;
		my @tmp=split /\n/,$line;
		if(length($tmp[1])>50){# if input read length exceeds 50bp (such as those from PE150 platform), we will chop the sequence to the first 50-bp for the following analysis
			$tmp[1]=substr($tmp[1],0,50);
			$tmp[3]=substr($tmp[3],0,50);
		}
		for my $i(0..$#site){
			if($tmp[1]=~s/^(\S*?)($site[$i])\S*/$2/){#取的是第一个酶切序列，非贪婪匹配
				$enzyme_reads_num++;#有酶切位点的reads
				my $begin=length($1);
				my $len=length($2);
				$tmp[3]=substr($tmp[3],$begin,$len);# quality score
				my $sub_line=join("\n",@tmp[0..3]);
				if($q_control eq "yes"){# check if QC needed
					next unless(&CheckN($sub_line));
					next unless(&CheckQ($sub_line));
				}
				$qc_reads_num++;#质控后的reads
				if($format eq "fa"){
					$tmp[0]=~s/^@/>/;
					print OU "$tmp[0]\n$tmp[1]\n";
				}elsif($format eq "fq"){
					print OU "$sub_line\n";
				}
				last;
			}
		}
	}
	close IN;
	close OU;
	$percent_sub=sprintf "%.2f",$qc_reads_num/$input_reads_num*100;
	print STAT "$outprefix\t$enzyme\t$input_reads_num\t$enzyme_reads_num\t$qc_reads_num\t$percent_sub%\n";
	close STAT;
}

sub CheckN{
    my $line=shift;
    my @tmp=split /\n/,$line;
    my $length=length($tmp[1]);
    @tmp=split //,$tmp[1];
    my $count=0;
    for my $base(@tmp){
        if($base eq "N"){$count++;};
    }
    if($ncount>0 && $ncount<1){
        if($ncount>= $count/$length){
            return 1;
        }else{
            return 0;
        }
    }elsif($ncount==0 || $ncount>=1){
        if($ncount>=$count){
            return 1;
        }else{
            return 0;
        }
    }
}

sub CheckQ{
    my $line = shift;
    my @array = split /\n/,$line;
    @array = split //,$array[3];
    my $count = 0;
    foreach my $i( @array ){
        next unless( ord($i) >= $quality + $qbase );
        $count ++;
    }
    if( $count >= scalar(@array) * $percent / 100 ){
        return 1;
    }else{
        return 0;
    }
}

sub Electronic_enzyme{
	my $genome=$input[0];
	my $outprefix=$outprefix[0];
	my $cnt=0;
	$/=">";
	if($genome=~/\.gz$/){
		open IN,"gzip -dc $genome|" or die "cannot open $genome\n";
	}else{
		open IN,"$genome" or die "cannot open $genome\n";
	}
	if($gz eq "yes"){
		open OU,"|gzip > $outdir/$outprefix.$enzyme.fa.gz" or die "cannot open $outdir/$outprefix.$enzyme.fa.gz\n";
	}else{
		open OU,">$outdir/$outprefix.$enzyme.fa" or die "cannot open $outdir/$outprefix.$enzyme.fa\n";
	}
	open STAT,">$outdir/$outprefix.$enzyme.stat.xls" or die "cannot open $outdir/$outprefix.$enzyme.stat.xls\n";
	print STAT "sample\tenzyme\tinput_reads_num\tenzyme_reads_num\tpercent\n";
	my($input_reads_num,$enzyme_reads_num,$percent_sub);
	$enzyme_reads_num=0;
	<IN>;
	while(<IN>){
		chomp;
		$input_reads_num++;
		my @tmp=split /\n/;
		my $id=(split /\s+/,$tmp[0])[0];
		my $seq=join("",@tmp[1..$#tmp]);
		$seq=uc($seq); # convert the lowercase to uppercase bases 小写碱基转换为大写
		my %hash;my %hash_tag;
		for my $i(0..$#site){ # iterate all restriction sites 循环酶切位点
			while($seq=~/($site[$i])/g){ # digital digestion
				my $tag=$1;
				my $len=length($tag);
				my $pos=pos($seq);
				$pos=$pos-$len+1;
				my $pos_end=$pos+$len-1;
				pos($seq)=$pos; #调整位置
				$hash{$pos}{$pos_end}="$id-$pos-$pos_end";
				$hash_tag{$pos}{$pos_end}=$tag;
			}
		}
		# sort the 2b-RAD tags by the genome positions排序后输出
		for my $pos(sort {$a<=>$b} keys %hash){
			for my $pos_end(sort {$a<=>$b} keys %{$hash{$pos}}){
				$cnt++;
				$enzyme_reads_num++;
				print OU ">$hash{$pos}{$pos_end}-$cnt\n"; #染色体id-起始位置-终止位置-第n个标签
				print OU "$hash_tag{$pos}{$pos_end}\n";
			}
		}
		undef %hash;
		undef %hash_tag;
	}
	close IN;
	close OU;
	$percent_sub=sprintf "%.2f",$enzyme_reads_num/$input_reads_num*100;
	print STAT "$outprefix\t$enzyme\t$input_reads_num\t$enzyme_reads_num\t$percent_sub%\n";
	close STAT;
	$/="\n";
}


sub execute{
	my $cmd = shift;
	print "$cmd\n";
	my $exit_code=system($cmd);
	if($exit_code!=0){
		print STDERR "Command $cmd failed with an exit code of $exit_code.\n";
		exit($exit_code >> 8);
	}
}
sub CheckDir{
	my $file = shift;
	unless( -d $file ){
		if( -d dirname($file) && -w dirname($file) ){system("mkdir $file");}
		else{print STDERR "$file not exists and cannot be built\n";exit 1;}
		}
		return 1;
}