Illumina_MySeq_Trim

#!/usr/bin/env perl
use strict;
use warnings;
#--INCLUDE PACKAGES-----------------------------------------------------------
use IO::String;
use File::Basename;
use File::Copy;
use Cwd;
#-----------------------------------------------------------------------------
#----SUBROUTINES--------------------------------------------------------------
#-----------------------------------------------------------------------------

sub get_file_data
    {
    my ($file_name) = @_;
    my @file_content;
    open (PROTEINFILE, $file_name);
    @file_content = <PROTEINFILE>;
    close PROTEINFILE;
    return @file_content;
    } 


    sub WriteArrayToFile
    {
    my ($filename, @in) = @_;
    my $a = join (@in, "\n");
    open (OUTFILE, ">$filename");
    foreach my $a (@in)
      {
      print OUTFILE $a;
      print OUTFILE "\n";
      }
    close (OUTFILE);
     }

sub parse_flags
{
  my ($flags) = @_;

  my @elements = split("-", $flags);
  my @r1; my @r2;
  foreach my $a (@elements)
      {
      $a=~ s/\s+$//;
      if ($a ne '')
         {
          my @tr = split (/ /, $a);
          if ($tr[0] && $tr[1]) { push (@r1, $tr[0]); push (@r2, $tr[1]);}
         }
      }
return (\@r1, \@r2);
} # end of parse_flags subroutine;

#-----------------------------------------------------------------------------
#-----------------------------------------------------------------------------
#-----------------------------------------------------------------------------
my @adapter_array;
my $system_string;

my $ForwardReads = $ARGV[0];
my $ReverseReads = $ARGV[1];
my $ResultDirectory = $ARGV[2];
if (not defined $ARGV[2] ) {
    $ResultDirectory = cwd();
}
my $on_off_switch = $ARGV[3];
if (not defined $ARGV[3] ) {
    $on_off_switch = "on";
#    $on_off_switch = "off";
}
my $runtime_flags = $ARGV[4];

#-----------------------------------------------------------------------------
# DEFAULT RUNTIME PARAMETERS
# we assume that data are MiSeq PE250 as default parameters
#-----------------------------------------------------------------------------


my $forward_read_length = 250; #default: 250 
my $reverse_read_length = 250; #default: 250


#-----------------------------------------------------------------------------
# RUNTIME FLAG MANAGEMENT
#
# -f_read_trim nt     trim 3'-end of forward reads so they are maximally nt nucleotides in length      
# -r_read_trim nt     trim 3'-end of reverse reads so they are maximally nt nucleotides in length   
# -transcriptome nt   indicate the data are RNAseq data with a read length of nt
# -nextera 0/1        to indicate if nextera transposons are present 1=yes 0=no
#
#-----------------------------------------------------------------------------

my $is_transcriptome =0; 
my $contains_nextera =0;

my @rf1 = @{(parse_flags($runtime_flags))[0]};
my @rf2 = @{(parse_flags($runtime_flags))[1]};

for (my $rf_count = 0; $rf_count <=scalar(@rf1); $rf_count++)
    {
    if ($rf1[$rf_count] eq "f_read_trim")  { $forward_read_length = $rf2[$rf_count]; }           
    if ($rf1[$rf_count] eq "r_read_trim")  { $reverse_read_length = $rf2[$rf_count]; }           
    if ($rf1[$rf_count] eq "nextera")  { $contains_nextera = 1; }           
    if ($rf1[$rf_count] eq "transcriptome")  { $forward_read_length = $rf2[$rf_count]; $reverse_read_length = $rf2[$rf_count]; $contains_nextera = 0; $is_transcriptome =1 }           

}

#-----------------------------------------------------------------------------
# test if reads zipped; if yes, unzip
#-----------------------------------------------------------------------------

my @testzip = split (/\./,$ForwardReads);
my $test_zip_length = scalar (@testzip);
if ($testzip[$test_zip_length-1] eq "gz") 
    {
    system ("gunzip ".$ForwardReads);
    $ForwardReads =~ s/.gz//g;
    printf $ForwardReads;
    }
@testzip = split (/\./,$ReverseReads);
$test_zip_length = scalar (@testzip);
if($testzip[$test_zip_length-1] eq "gz") 
    {
    system ("gunzip ".$ReverseReads);
    $ReverseReads =~ s/.gz//g;
    printf $ReverseReads;
    }

#-----------------------------------------------------------------------------
# Run fastQC
#-----------------------------------------------------------------------------

system("mkdir -p ".$ResultDirectory."/temp");
#system("fastqc ".$ForwardReads." -o ".$ResultDirectory."/temp");
#system("fastqc ".$ReverseReads." -o ".$ResultDirectory."/temp");
system("fastqc ".$ForwardReads." -o ".$ResultDirectory);  #not putting them into temp directory any more
system("fastqc ".$ReverseReads." -o ".$ResultDirectory);

#-----------------------------------------------------------------------------
# on/off switch
#-----------------------------------------------------------------------------

if ($on_off_switch eq "on")
{

#-----------------------------------------------------------------------------
# remove nextera transposons with trim_galore
#-----------------------------------------------------------------------------

my $f_b = $ForwardReads; $f_b =~ s/.fastq//;
my $r_b = $ReverseReads; $r_b =~ s/.fastq//;

if ($contains_nextera == 1)

{

   $system_string = "trim_galore --phred33 --paired -q 30 --nextera -o ".$ResultDirectory."/ ".$ForwardReads." ".$ReverseReads;
   system($system_string);
#   my $f_b = $ForwardReads; $f_b =~ s/.fastq//;
#   my $r_b = $ReverseReads; $r_b =~ s/.fastq//;

   $system_string = ("mv ".$ForwardReads." ".$f_b.".old");
   system ($system_string);
   printf $system_string."\n\n";

   $system_string =  ("mv ".$ReverseReads." ".$r_b.".old");
   system ($system_string);
   printf $system_string."\n\n";

   $system_string = ("mv ".$f_b."_val_1.fq ".$f_b.".fastq");
   system ($system_string);
   printf $system_string."\n\n";

   $system_string = ("mv ".$r_b."_val_2.fq ".$r_b.".fastq");
   system ($system_string);
   printf $system_string."\n\n";

}
#-------------get forward adapters--------------------------------------------
#-----------------------------------------------------------------------------

my $f = $ForwardReads."c"; $f =~ tr/.fastqc/_fastqc/;
my ($filename, $dirs, $suffix) = fileparse($f);
 
#system("unzip ".$ResultDirectory."/temp/".$filename.".zip -d ".$ResultDirectory."/temp/");

system("unzip ".$ResultDirectory."/".$filename.".zip -d ".$ResultDirectory."/temp/");
my @Fanalysis = get_file_data ($ResultDirectory."/temp/".$filename."/fastqc_data.txt");
($filename, $dirs, $suffix) = fileparse($ForwardReads);
my $Ffilename = $filename;
my $active =0;
foreach my $s (@Fanalysis)
  {
  my @ug = split (/\t/, $s);
  if ($ug[0] eq ">>Overrepresented sequences")
     {$active = 1;next;}
  if ($s eq "#Sequence	Count	Percentage	Possible Source\n")
     {next;} 
  if ($s eq ">>END_MODULE\n")
     {$active = 0;next;}
  if ($active == 1)
     {
     my @tt = split (/\t/, $s); 
     if (length $tt[0] >5)
        {
        push (@adapter_array, "-b ".$tt[0]);
        printf ("\n".$tt[0]);
        }
     }  
  }

#-----------------------------------------------------------------------------
#---------get reverse adapters------------------------------------------------
#-----------------------------------------------------------------------------


$f = $ReverseReads."c"; $f =~ tr/.fastqc/_fastqc/;
($filename, $dirs, $suffix) = fileparse($f);

#system("unzip ".$ResultDirectory."/temp/".$filename.".zip -d ".$ResultDirectory."/temp/");

system("unzip ".$ResultDirectory."/".$filename.".zip -d ".$ResultDirectory."/temp/");
my @Ranalysis = get_file_data ($ResultDirectory."/temp/".$filename."/fastqc_data.txt");
($filename, $dirs, $suffix) = fileparse($ReverseReads);
my $Rfilename = $filename;

$active =0;
foreach my $s (@Ranalysis)
  {
  my @ug = split (/\t/, $s);  
  if ($ug[0] eq ">>Overrepresented sequences")
     {$active = 1;next;}
  if ($s eq "#Sequence	Count	Percentage	Possible Source\n")
     {next;}
  if ($s eq ">>END_MODULE\n")
     {$active = 0;next;}
  if ($active == 1)
     {
     my @tt = split (/\t/, $s);
     if (length $tt[0] >5)
        {
        push (@adapter_array, "-b ".$tt[0]);
       printf ("\n".$tt[0]);
        }
     }
  }

push (@adapter_array, "-b GATCGGAAGAGCACACGTCTGAACTCCAGTCAC");
#push (@adapter_array, "-b AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCT");
push (@adapter_array, "-b AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT");
#push (@adapter_array, "-b AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATCATT");
push (@adapter_array, "-b NNNNNNN");
push (@adapter_array, "-e 0.15");
push (@adapter_array, "-q 30");
push (@adapter_array, "-n 2");


#-----------------------------------------------------------------------------
#--write adapter file---------------------------------------------------------
#-----------------------------------------------------------------------------

my $outfile_name = $ResultDirectory."/temp/".$Ffilename.".conf";

WriteArrayToFile ($outfile_name, @adapter_array);

#-----------------------------------------------------------------------------
#--trim adapter sequences-----------------------------------------------------
#-----------------------------------------------------------------------------
printf "\n\n";
my $confstring;
foreach my $ty (@adapter_array) { $confstring .= $ty." ";}
system ("cutadapt ".$confstring.$ForwardReads." > ".$ResultDirectory."/temp/".$Ffilename.".cutadapt");
system ("cutadapt ".$confstring.$ReverseReads." > ".$ResultDirectory."/temp/".$Rfilename.".cutadapt");


#-----------------------------------------------------------------------------
#-trim to Q30-----------------------------------------------------------------
#-----------------------------------------------------------------------------

if ($is_transcriptome == 0)
{
$confstring = "read_fastq -e base_33 -i ".$ResultDirectory."/temp/".$Ffilename.".cutadapt | trim_seq -m 30 -l 10 --trim=right | write_fastq -o ".$ResultDirectory."/temp/".$Ffilename.".cutadapt.q30 -x\n";
system ($confstring);
$confstring = "read_fastq -e base_33 -i ".$ResultDirectory."/temp/".$Rfilename.".cutadapt | trim_seq -m 30 -l 10 --trim=right | write_fastq -o ".$ResultDirectory."/temp/".$Rfilename.".cutadapt.q30 -x\n";
system ($confstring);
}

if ($is_transcriptome == 1)
{
$confstring = "read_fastq -e base_33 -i ".$ResultDirectory."/temp/".$Ffilename.".cutadapt | trim_seq -m 30 | write_fastq -o ".$ResultDirectory."/temp/".$Ffilename.".cutadapt.q30 -x\n";
system ($confstring);
$confstring = "read_fastq -e base_33 -i ".$ResultDirectory."/temp/".$Rfilename.".cutadapt | trim_seq -m 30 | write_fastq -o ".$ResultDirectory."/temp/".$Rfilename.".cutadapt.q30 -x\n";
system ($confstring);
}

#-----------------------------------------------------------------------------
#- remove poly-A and crop to defined length----------
#-----------------------------------------------------------------------------
chdir($ResultDirectory."/temp/");
$confstring = "homerTools trim -3 AAAAAAAA -len ".$forward_read_length." ".$ResultDirectory."/temp/".$Ffilename.".cutadapt.q30\n";
system($confstring);
$confstring = "homerTools trim -3 AAAAAAAA -len ".$reverse_read_length." ".$ResultDirectory."/temp/".$Rfilename.".cutadapt.q30\n";
system($confstring);

#-----------------------------------------------------------------------------
#-MINLEN: 50; paired only----------------------------------
#-----------------------------------------------------------------------------

#$confstring = "trimmomatic PE -phred33 ".$ResultDirectory."/temp/".$Ffilename.".cutadapt.q30.trimmed ".$ResultDirectory."/temp/".$Rfilename.".cutadapt.q30.trimmed ".$Ffilename.".trimmed.paired.fq.gz.fastq ".$Ffilename.".trimmed.unpaired.fq.gz ".$Rfilename.".trimmed.paired.fq.gz.fastq ".$Rfilename.".trimmed.unpaired.fq.gz MINLEN:50\n";

$confstring = "trimmomatic PE -phred33 ".$ResultDirectory."/temp/".$Ffilename.".cutadapt.q30.trimmed ".$ResultDirectory."/temp/".$Rfilename.".cutadapt.q30.trimmed ".$Ffilename.".trimmed.paired.fq.gz.fastq ".$Ffilename.".trimmed.unpaired.fq.gz ".$Rfilename.".trimmed.paired.fq.gz.fastq ".$Rfilename.".trimmed.unpaired.fq.gz MINLEN:50\n";

printf "\n\n".$confstring."\n\n";
system($confstring);

#-----------------------------------------------------------------------------
#-Run qc on paired files
#-----------------------------------------------------------------------------
system("fastqc ".$ResultDirectory."/temp/".$Ffilename.".trimmed.paired.fq.gz.fastq -o ".$ResultDirectory);
system("fastqc ".$ResultDirectory."/temp/".$Rfilename.".trimmed.paired.fq.gz.fastq -o ".$ResultDirectory);

#$confstring = "cp ".$ResultDirectory."/temp/".$Ffilename.".trimmed.paired.fq.gz.fastq ".$ResultDirectory."/".$Ffilename.".QCed.fastq";
#$confstring = "cp ".$ResultDirectory."/temp/".$Rfilename.".trimmed.paired.fq.gz.fastq ".$ResultDirectory."/".$Rfilename.".QCed.fastq";

$confstring = "cp ".$ResultDirectory."/temp/".$Ffilename.".trimmed.paired.fq.gz.fastq ".$ResultDirectory."/F.QCed.fastq";
system($confstring);
$confstring = "cp ".$ResultDirectory."/temp/".$Rfilename.".trimmed.paired.fq.gz.fastq ".$ResultDirectory."/R.QCed.fastq";
system($confstring);


#-----------------------------------------------------------------------------
#Convert to fasta
#-----------------------------------------------------------------------------

if ($is_transcriptome == 0)
{
$confstring = "read_fastq -i ".$ResultDirectory."/F.QCed.fastq -e base_33 | write_fasta -o ".$ResultDirectory."/F.QCed.fasta -x";
system ($confstring);
$confstring = "read_fastq -i ".$ResultDirectory."/R.QCed.fastq -e base_33 | write_fasta -o ".$ResultDirectory."/R.QCed.fasta -x";
system ($confstring);
}

#-----------------------------------------------------------------------------
#-----------------------------------------------------------------------------

#-----------------------------------------------------------------------------
#clearnup
#-----------------------------------------------------------------------------

system ("rm -rf ".$ResultDirectory."/temp/");
system ("rm ".$ResultDirectory."/*.unpaired.fq.gz");
system ("rm ".$ResultDirectory."/*.zip");

#system ("rm ".$ForwardReads);
#system ("rm ".$ReverseReads);

system ("rm ".$f_b.".old");
system ("rm ".$r_b.".old");


#-----------------------------------------------------------------------------
# on/off switch
#-----------------------------------------------------------------------------
} 
#end of on switch braket

if ($on_off_switch eq "off")
{

if ($ForwardReads ne $ResultDirectory."/F.QCed.fastq")
{
  $system_string = ("mv ".$ForwardReads." ".$ResultDirectory."/F.QCed.fastq");
  system ($system_string);
}
if ($ReverseReads ne $ResultDirectory."/R.QCed.fastq")
{
  $system_string =  ("mv ".$ReverseReads." ".$ResultDirectory."/R.QCed.fastq");
  system ($system_string);
}

$system_string = "read_fastq -i ".$ResultDirectory."/R.QCed.fastq -e base_33 | write_fasta -o ".$ResultDirectory."/R.QCed.fasta -x";
system ($system_string);


$system_string = "read_fastq -i ".$ResultDirectory."/F.QCed.fastq -e base_33 | write_fasta -o ".$ResultDirectory."/F.QCed.fasta -x";
system ($system_string);

}


#-----------------------------------------------------------------------------
# gzip files if transcriptome
#-----------------------------------------------------------------------------


if ($is_transcriptome == 1)
{
system ("gzip ".$ResultDirectory."/*.fastq")
}