<commit_before>#!/usr/bin/perl -w

use File::Which;
use Cwd qw(abs_path);
use Getopt::Long;

# INITIALIZE PATH AND FLAGS--TSW
my $binPath = abs_path($0);
$0 =~ s/^.*\///;
$binPath =~ s/\/$0$//;

my $bowtie = "bowtie"; # by default;
my $species = "Hsa"; # by default;
my $dbDir = "$binPath/../$species"; #default
my $pairedEnd = 0; # no by default
my $runExprFlag = 0; # no by default
my $onlyExprFlag = 0; # no by default
my $trim = "once"; # 1 by default
my $cores = 1; #default
my $readLength; 

my $legacyFlag = 0;
my $verboseFlag = 0;

GetOptions("bowtieProg=s" => \$bowtie,
			  "sp=s" => \$species,
			  "db=s" => \$dbDir,
			  "c" => \$cores,
			  "pe" => \$pairedEnd,
			  "expr" => \$runExprFlag,
			  "exprONLY" => \$onlyExprFlag,
			  "trim=s" => \$trim,
			  "help" => \$helpFlag,
			  "legacy" => \$legacyFlag,
			  "verbose" => \$verboseFlag,
			  "readLen=i" => \$readLength);


sub sysErrMsg {
  my $sysCommand = shift;
  not system($sysCommand) or die "[vastdb align error]: $sysCommand Failed in $0!";
}

sub verbPrint {
  my $verbMsg = shift;
  if($verboseFlag) {
    chomp($verbMsg);
    print STDERR "[vastdb align]: $verbMsg\n";
  }
}

# Set up output directory structure
unless($legacyFlag) {
  mkdir("align_out") unless (-e "align_out");
  mkdir("expr_out") unless (-e "expr_out");
  mkdir("align_out/$species") unless (-e "align_out/$species");
  mkdir("expr_out/$species") unless (-e "align_out/$species");
}

# Use pigz if installed  --KH
my $zip = which('pigz');
if ($zip eq '') {
    $zip = which('gzip');
    if ($zip eq '') {
        die "Error: need gzip or pigz\n";
    }
} else {
    verbPrint "Found pigz...";
    $zip .= " -fp $cores ";
}

# Command line flags here
if($pairedEnd and !defined($ARGV[0]) and !defined($ARGV[1])) { $helpFlag = 1; }


if (!$ARGV[0] or $helpFlag){
    print "\nCommand: \nvastdb align fastq_file_1 [fastq_file_2] -sp Mmu/Hsa [-expr/-exprONLY -trim once/twice -pe -c 1 -bowtieProg bowtie]\n\n";
    
    print ">> fastq_file can be compressed or uncompressed. It MUST have this name format: Sample-readlength.fq\n";
    print ">>>> Only underscores (\"_\") are allowed in the Sample name. No dashes, dots, spaces, etc.\n";
    print ">>>> If running with the paired end option (\"\-PE\"), they must be called: Sample_1-length.fq Sample_2-length.fq\n";
    print ">> The process can be started with \"genome substracted\" samples if a Sample-lenght-e.fq is used\n";
    print ">> Species implemented: Human (hg19) or Mouse (mm9)\n";
    print ">> For expression analyses: -expr (PSIs plus cRPKM calculations) OR -exprONLY (only cRPKMs) OR none\n";
    print ">> For trimming, it can be trimmed once (at 3') or twice (in an overlapping manner).\n";
    print ">>>> If nothing is provided, the program will decide based on the length of the reads (Default is twice if length>50)\n";
    print ">> -c N: Number of threds used when running bowtie (default=1).\n";
    die ">> Recommended to allow at least 15GB of RAM (~10GB are needed for mapping to the genome). For large files (~1 lane), >25GB\n\n";
}

die "Needs species\n" if !$species;



### Getting sample name and length:
$file=$ARGV[0];
my $zipped = ($file =~ /\.gz$/) ? 1 : 0;
if ($file=~/\-e\.f/){
    $genome_sub=1;
    ($root,$length)=$file=~/(.+?)\-(.+?)\-e\.(fastq|fq)(\.gz)?/;
    $fq=$&;
    die "Only for 50nt or 36nt if genome substracted\n" if $length!=36 && $length!=50;
}
else {
    if ($pairedEnd){
	($root,$length)=$file=~/(.+?)\_1\-(.+?)\.(fastq|fq)(\.gz)?/;
	$fq1=$&;
	$fq2=$ARGV[1];
    #$fq2=~s/\.gz//;
	$fq="$root-$length.fq";
    }
    else {
	($root,$length)=$file=~/(.+?)\-(.+?)\.(fastq|fq)(\.gz)?/;
	$fq=$&;
    }
}
###

#length options:
if ($length>=50){
    $difLE=$length-50;
    $le=50;
}
elsif ($length>=36){
    $difLE=$length-36;
    $le=36;
}
else {
    die "Minimum reads length: 50nt (Human) and 36nt (Mouse)\n";
}
die "Reads <50nt not available for Human\n" if $le==36 && $species eq "Hsa";
#####

#sysErrMsg "gunzip $file" if $file=~/\.gz/;
#sysErrMsg "gunzip $file2" if $file2=~/\.gz/ && $pairedEnd;

if (!$genome_sub){
#### Expression analysis (it maps only the first $le nucleotides of the read)
 if ($runExprFlag || $onlyExprFlag){
     verbPrint "Mapping reads against mRNA sequences\n";

     my $cmd = "$bowtie -p $cores -m 1 -v 2 -3 $difLE $dbDir/EXPRESSION/mRNA - expr_out/$species"."mRNA-$le-$root.out";

     if (!$pairedEnd) {
         $cmd = "$fq | $cmd";
     } else {
         $cmd = "$fq1 | $cmd";
     }

     if ($zipped) {
         $cmd = "gzip -dc $cmd";
     } else {
         $cmd = "cat $cmd";
     }

     sysErrMsg $cmd;
     verbPrint "Calculating cRPKMs\n";
     sysErrMsg "$binPath/expr_RPKM.pl $dbDir/EXPRESSION/$species"."mRNA-$le-$root.out expr_out/$species"."_mRNA-$le.eff"; 
 }
 if ($onlyExprFlag){
     #print "Compressing raw fastq files\n";
     #sysErrMsg "gzip $fq" if !$pairedEnd;
     #sysErrMsg "gzip $fq1 $fq2" if $pairedEnd;
     die "Expression analysis done\n";
 }
###

#### Merge PE
 if ($pairedEnd){
     verbPrint "Concatenating paired end reads\n";
     sysErrMsg "cat $fq1 $fq2 > $fq";
     #sysErrMsg "gzip $fq1 $fq2";
 }
 
#### Trimming
 if ($difLE>=10){
     if ($trim eq "twice" || !$trim){
	 if ($length > ($le*2)+10){
	     $half_length=sprintf("%.0f",$length/2);
	     verbPrint "Trimming and splitting fastq sequences from $length to $half_length nt\n";
	     sysErrMsg "$binPath/Trim-twiceOVER.pl $fq $half_length";
	     verbPrint "Trimming and splitting fastq sequences to $le nt sequences\n";
	     sysErrMsg "$binPath/Trim-twiceOVER.pl $root-$length-$half_length.fq $le";
	     sysErrMsg "rm $root-$length-$half_length.fq";
	     sysErrMsg "mv $root-$length-$half_length-$le.fq $root-$le.fq";
	 }
	 else {
	     verbPrint "Trimming and splitting fastq sequences to $le nt sequences\n";
	     sysErrMsg "$binPath/Trim-twiceOVER.pl $fq $le";
	     sysErrMsg "mv $root-$length-$le.fq $root-$le.fq";
	 }
     }
     elsif ($trim eq "once"){
	 verbPrint "Trimming fastq sequences to $le nt sequences\n";
	 sysErrMsg "$binPath/Trim-once.pl $fq $le";
	 sysErrMsg "mv $root-$length-$le.fq $root-$le.fq";
     }
 }
 verbPrint "Compressing raw fastq file\n" unless $length==50 || $length==36;
 #sysErrMsg "gzip $fq" unless $length==50 || $length==36;
####
 
#### Get effective reads (i.e. genome substraction).
 verbPrint "Doing genome substraction\n";
 sysErrMsg "$bowtie -p $cores -m 1 -v 2 --un $root-$le-e.fq --max /dev/null $dbDir/FILES/gDNA $root-$le.fq /dev/null";
 verbPrint "Compressing trimmed reads\n";
 sysErrMsg "$zip $root-$le.fq";<commit_msg>fixed fq1 and fq2 filename parsing, can now handle fastq files from any directory<commit_after>#!/usr/bin/perl -w

use File::Which;
use Cwd qw(abs_path);
use Getopt::Long;

# INITIALIZE PATH AND FLAGS--TSW
my $binPath = abs_path($0);
$0 =~ s/^.*\///;
$binPath =~ s/\/$0$//;

my $bowtie = "bowtie"; # by default;
my $species = "Hsa"; # by default;
my $dbDir = "$binPath/../$species"; #default
my $pairedEnd = 0; # no by default
my $runExprFlag = 0; # no by default
my $onlyExprFlag = 0; # no by default
my $trim = "once"; # 1 by default
my $cores = 1; #default
my $readLength; 

my $legacyFlag = 0;
my $verboseFlag = 0;

GetOptions("bowtieProg=s" => \$bowtie,
			  "sp=s" => \$species,
			  "db=s" => \$dbDir,
			  "c" => \$cores,
			  "pe" => \$pairedEnd,
			  "expr" => \$runExprFlag,
			  "exprONLY" => \$onlyExprFlag,
			  "trim=s" => \$trim,
			  "help" => \$helpFlag,
			  "legacy" => \$legacyFlag,
			  "verbose" => \$verboseFlag,
			  "readLen=i" => \$readLength);


sub sysErrMsg {
  my $sysCommand = shift;
  not system($sysCommand) or die "[vastdb align error]: $sysCommand Failed in $0!";
}

sub verbPrint {
  my $verbMsg = shift;
  if($verboseFlag) {
    chomp($verbMsg);
    print STDERR "[vastdb align]: $verbMsg\n";
  }
}

# Set up output directory structure
unless($legacyFlag) {
  mkdir("align_out") unless (-e "align_out");
  mkdir("expr_out") unless (-e "expr_out");
  mkdir("align_out/$species") unless (-e "align_out/$species");
  mkdir("expr_out/$species") unless (-e "align_out/$species");
}

# Use pigz if installed  --KH
my $zip = which('pigz');
if ($zip eq '') {
    $zip = which('gzip');
    if ($zip eq '') {
        die "Error: need gzip or pigz\n";
    } 
} else {
    $zip .= " -fp $cores ";
}
verbPrint "Found $zip...";

# Command line flags here
if($pairedEnd and !defined($ARGV[0]) and !defined($ARGV[1])) { $helpFlag = 1; }

if (!defined($ARGV[0]) or $helpFlag){
    print "\nCommand: \nvastdb align fastq_file_1 [fastq_file_2] -sp Mmu/Hsa [-expr/-exprONLY -trim once/twice -pe -c 1 -bowtieProg bowtie]\n\n";
    
    print ">> fastq_file can be compressed or uncompressed. It MUST have this name format: Sample-readlength.fq\n";
    print ">>>> Only underscores (\"_\") are allowed in the Sample name. No dashes, dots, spaces, etc.\n";
    print ">>>> If running with the paired end option (\"\-PE\"), they must be called: Sample_1-length.fq Sample_2-length.fq\n";
    print ">> The process can be started with \"genome substracted\" samples if a Sample-lenght-e.fq is used\n";
    print ">> Species implemented: Human (hg19) or Mouse (mm9)\n";
    print ">> For expression analyses: -expr (PSIs plus cRPKM calculations) OR -exprONLY (only cRPKMs) OR none\n";
    print ">> For trimming, it can be trimmed once (at 3') or twice (in an overlapping manner).\n";
    print ">>>> If nothing is provided, the program will decide based on the length of the reads (Default is twice if length>50)\n";
    print ">> -c N: Number of threds used when running bowtie (default=1).\n";
    die ">> Recommended to allow at least 15GB of RAM (~10GB are needed for mapping to the genome). For large files (~1 lane), >25GB\n\n";
}

die "Needs species\n" if !$species;



### Getting sample name and length:
my $fq1 = $ARGV[0];
my $fq2;
my $fileName1 = $fq1;
my $fileName2;
my $zipped = ($fq1 =~ /\.gz$/) ? 1 : 0;

my($root, $length);

$fileName1 =~ s/^.*\///g;

if ($fileName1 =~ /\-e\.f/){
    $genome_sub=1;
    ($root,$length)=$fileName1=~/(.+?)\-(.+?)\-e\.(fastq|fq)(\.gz)?/;
    $fq=$&;
    die "Only for 50nt or 36nt if genome substracted\n" if $length!=36 && $length!=50;
} else {
    # allow readlength to be given by -readLen x --TSW
    if($readLength) {
         $length = $readLength;
         $fileName1 =~ /(.+?)\.(fastq|fq)(\.gz)?/;
         $root = $1;
    } else { # default behavior by --MI
         ($root,$length)=$fileName1=~/(.+?)\_{0,1}1{0,1}\-(.+?)\.(fastq|fq)(\.gz)?/;
    }
    if ($pairedEnd){
		$fq2 = $ARGV[1];
		$fileName2 = $fq2;
      $fileName2 =~ s/^.*\///g;
    }
}

#verbPrint "$fileName1\n$fq1\n;"; die "";
###

#length options:
if ($length>=50){
    $difLE=$length-50;
    $le=50;
} elsif ($length>=36){
    $difLE=$length-36;
    $le=36;
} else {
    die "Minimum reads length: 50nt (Human) and 36nt (Mouse)\n";
}
die "Reads <50nt not available for Human\n" if $le==36 && $species eq "Hsa";
#####

#sysErrMsg "gunzip $file" if $file=~/\.gz/;
#sysErrMsg "gunzip $file2" if $file2=~/\.gz/ && $pairedEnd;

if (!$genome_sub){
#### Expression analysis (it maps only the first $le nucleotides of the read)
 if ($runExprFlag || $onlyExprFlag){
     verbPrint "Mapping reads against mRNA sequences\n";

     my $cmd = "$bowtie -p $cores -m 1 -v 2 -3 $difLE $dbDir/EXPRESSION/mRNA - expr_out/$species"."_mRNA-$le-$root.out";

     if ($pairedEnd) {
         $cmd = "$fq1 $fq2 | $cmd";  # altered this to cat both for/rev reads into bowtie.
     } else {
         $cmd = "$fq1 | $cmd";
     }

     if ($zipped) {
			$cmd = "zcat $cmd";
	 # Is zcat a better choice?  --TSW
    #     $cmd = "gzip -dc $cmd";
     } else {
         $cmd = "cat $cmd";
     }

     sysErrMsg $cmd;
     verbPrint "Calculating cRPKMs\n";
     sysErrMsg "$binPath/expr_RPKM.pl $dbDir/EXPRESSION/$species"."mRNA-$le-$root.out expr_out/$species"."_mRNA-$le.eff"; 
 }
 if ($onlyExprFlag){
     #print "Compressing raw fastq files\n";
     #sysErrMsg "gzip $fq" if !$pairedEnd;
     #sysErrMsg "gzip $fq1 $fq2" if $pairedEnd;
     print STDERR "Expression analysis done\n";
     exit(1);
 }
###

#### Merge PE
 if ($pairedEnd){
     verbPrint "Concatenating paired end reads\n";
     sysErrMsg "cat $fq1 $fq2 > $fq";  # away with this as well?
     #sysErrMsg "gzip $fq1 $fq2";
 } else {
   $fq = $fq1;
 }
 
#### Trimming
 if ($difLE >= 10){
   if ($trim eq "twice" || !$trim){
	  if ($length > ($le*2)+10){
	     $half_length = sprintf("%.0f", $length / 2);
	     verbPrint "Trimming and splitting fastq sequences from $length to $half_length nt\n";
	     sysErrMsg "$binPath/Trim-twiceOVER.pl $fq $half_length > $root-$length-$half_length.fq";
	     verbPrint "Trimming and splitting fastq sequences to $le nt sequences\n";
	     sysErrMsg "$binPath/Trim-twiceOVER.pl $root-$length-$half_length.fq $le > $root-$le.fq";
	     sysErrMsg "rm $root-$length-$half_length.fq";
	     #sysErrMsg "mv $root-$length-$half_length-$le.fq $root-$le.fq";  #piping to stdout removes need for this --TSW
	  } else {
	     verbPrint "Trimming and splitting fastq sequences to $le nt sequences\n";
	     sysErrMsg "$binPath/Trim-twiceOVER.pl $fq $le > $root-$le.fq";
	  #   sysErrMsg "mv $root-$length-$le.fq $root-$le.fq"; #piping to stdout removes need for this --TSW
	  }
   } elsif ($trim eq "once"){
	 verbPrint "Trimming fastq sequences to $le nt sequences\n";
	 sysErrMsg "$binPath/Trim-once.pl $fq $le > $root-$le.fq";
	# sysErrMsg "mv $root-$length-$le.fq $root-$le.fq"; #piping to stdout removes need for this --TSW
   }
 }
 verbPrint "Compressing raw fastq file\n" unless $length==50 || $length==36;
 #sysErrMsg "gzip $fq" unless $length==50 || $length==36;
####
 
#### Get effective reads (i.e. genome substraction).
 verbPrint "Doing genome substraction\n";
 sysErrMsg "$bowtie -p $cores -m 1 -v 2 --un $root-$le-e.fq --max /dev/null $dbDir/FILES/gDNA $root-$le.fq /dev/null";
 verbPrint "Compressing trimmed reads\n";
 sysErrMsg "$zip $root-$le.fq";
{'fraction_non_alphanumeric': 0.1895499892341922, 'fraction_numerical': 0.016148711691667265, 'mean_word_length': 3.1781109445277362, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 97, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, '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': '296645', 'n_tokens_mistral': 5685, 'n_tokens_neox': 5206, 'n_words': 1472}