#!/usr/bin/perl

# Comamnd usage:
#   $ perl quality_filtering.pl [Fastq_indir] [FastQC_outdir] [Clean_Fastq_outdir]

use strict;
use warnings;
use autodie;

use File::Copy;
use File::Path;
use File::Basename;

my ($fastq_indir, $rRNA_infile, $fastqc_outdir, $clean_fastq_outdir) = @ARGV;

mkpath([$fastqc_outdir, $clean_fastq_outdir]);

my $processor = qx(grep -c "processor" /proc/cpuinfo);
chomp($processor);

my $rRNA_build_file = dirname($rRNA_infile) . "/bowtie_build";
system("bowtie2-build --quiet $rRNA_infile $rRNA_build_file");

my @fastq_files = glob("$fastq_indir/*.fastq");

foreach my $fastq_file (@fastq_files){
	

	
	print "===============================================================\n";
	print "FastQC raw data quality checking...\n";
  system("fastqc --nogroup -t $processor --extract -o $fastqc_outdir $fastq_file");
  
  # . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
  # Check adapter contamination
  
  my $fastq_filesymbol = $1 if(basename($fastq_file) =~ /(.+)\.fastq/);
  
  my @adapter_seq;
  
  my $fastqc_file = "$fastqc_outdir/${fastq_filesymbol}_fastqc/fastqc_data.txt";
  open(FASTQC,$fastqc_file);
  my $sw = 0;
  while(<FASTQC>){ chomp;
		$sw = 0 if(/^>>END_MODULE/);
		if($sw == 1){
			my ($over_seq, $type) = (split /\t+/, $_)[0,3];
			push @adapter_seq, $over_seq if($type ne "No Hit");
		}	
		$sw = 1 if(/^#Sequence	Count	Percentage/);
	}
  close(FASTQC);
  
  rmtree(dirname($fastqc_file));
  
  # . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
  # Remove adapter sequence (If adapter sequence exists)
  
  my $fastq_adapter_file = "$clean_fastq_outdir/${fastq_filesymbol}_adapter.fastq";
  
  if(@adapter_seq >= 1){
		my $cutadapt_arg;
		$cutadapt_arg .= "-b $_ " foreach(@adapter_seq);
		print "$cutadapt_arg\n";
		print "===============================================================\n";
		print "Cutadapt adapter trimming...\n";
		system("cutadapt -m 20 $cutadapt_arg $fastq_file > $fastq_adapter_file");
	}else{
		copy($fastq_file, "$fastq_adapter_file");
	}
  
  # . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
  # Remove low quality read
  
  print "===============================================================\n";
  print "Prinseq quality filtering...\n";
  my $fastq_prinseq_file = "$clean_fastq_outdir/${fastq_filesymbol}_prinseq";
  system("prinseq-lite.pl -trim_qual_right 30 -min_qual_mean 30 -trim_tail_right 5 -min_len 20 -ns_max_p 5 -fastq $fastq_adapter_file -out_good $fastq_prinseq_file -out_bad null");
  
  # . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
  # Remove rRNA read
  
  print "===============================================================\n";
  print "Bowtie2 rRNA contamination filtering...\n";
  my $fastq_clean_file = "$clean_fastq_outdir/${fastq_filesymbol}_clean.fastq";
  my $fastq_rrna_file  = "$clean_fastq_outdir/${fastq_filesymbol}_rrna.fastq";
  system("bowtie2 --quiet -p $processor --un $fastq_clean_file -x $rRNA_build_file -q $fastq_prinseq_file.fastq > $fastq_rrna_file");
  
  # . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
  # Filtered fastq data quality check
  
  print "===============================================================\n";
	print "FastQC filtered data quality checking...\n";
  system("fastqc --nogroup -t $processor -o $fastqc_outdir $fastq_clean_file");
  
  unlink(glob("$fastqc_outdir/*.zip")); # Remove fastqc zip files
  unlink($fastq_adapter_file, "$fastq_prinseq_file.fastq", $fastq_rrna_file);
}

unlink(glob(dirname($rRNA_infile). "/*.bt2"));