<commit_before>#!/bin/perl
###############################################################################
#   
#	USAGE: perl RNA-SeqAlign.pl [options]
#
#  Created by <NAME>
#
###############################################################################
use strict;
use warnings;
use File::Basename; # enable manipulating of the full path
use Cwd;
use lib '/homes/bioinfo/bioinfo_software/perl_modules/File-Slurp-9999.19/lib';
use File::Slurp;
# use List::Util qw(max);
# use List::Util qw(sum);
# use Bio::SeqIO;
# use Bio::Seq;
# use Bio::DB::Fasta;
use Getopt::Long;
use Pod::Usage;
###############################################################################
##############         Print informative message             ##################
###############################################################################
print "###########################################################\n";
print "#  RNA-SeqAlign.pl                                        #\n";
print "#                                                         #\n";
print "#  Created by <NAME> 3/15/14                    #\n";
print "#  github.com/i5K-KINBRE-script-share                     #\n";
print "#  perl RNA-SeqAlign.pl -help # for usage/options         #\n";
print "#  perl RNA-SeqAlign.pl -man # for more details           #\n";
print "###########################################################\n";
###############################################################################
##############                get arguments                  ##################
###############################################################################
my ($r_list,$project_name,$transcriptome,$min_len,$clean_read_file1,$clean_read_file2,@clean_r1,@clean_r2,$labels,$sams,$clean_read_singletons,$out_dir);
my $convert_header = 0;
my $mapq = 10;
my $man = 0;
my $help = 0;
GetOptions (
			  'help|?' => \$help, 
			  'man' => \$man,
			  'r|r_list:s' => \$r_list,
			  't|transcriptome:s' => \$transcriptome,
              'p|project_name:s' => \$project_name,
			  'c|convert_header' => \$convert_header,
			  'm|mapq:i' => \$mapq,
			  'l|read_length:i' => \$min_len

              )  
or pod2usage(2);
pod2usage(1) if $help;
pod2usage(-exitstatus => 0, -verbose => 2) if $man;
sub quote { qq!"$_[0]"! } ## interpolate slurped text
my $dirname = dirname(__FILE__); # github directories (all github directories must be in the same directory)
my $home = getcwd; # working directory (this is where output files will be printed)
mkdir "${home}/${project_name}_scripts";
mkdir "${home}/${project_name}_qsubs";
mkdir "${home}/${project_name}_prinseq";
###############################################################################
############## Create array of the sample names and read files    #############
###############################################################################
my @reads;
open (READ_LIST, '<', $r_list) or die "Can't open $r_list!\n";
while (<READ_LIST>)
{
    chomp;
    push @reads , [split];
}
#######################################################################
#########     Index the de novo transcriptome for Bowtie2    ##########
#######################################################################
close (SCRIPT);
open (SCRIPT, '>', "${home}/${project_name}_scripts/${project_name}_index.sh") or die "Can't open ${home}/${project_name}_scripts/${project_name}_index.sh!\n"; # create a shell script for each read-pair set
open (QSUBS_INDEX, '>', "${home}/${project_name}_qsubs/${project_name}_qsubs_index.sh") or die "Can't open ${home}/${project_name}_qsubs/${project_name}_qsubs_index.sh!\n";
print QSUBS_INDEX "#!/bin/bash\n";
$dirname =~ /(.*)\/RNA-Seq-annotation-and-comparison\/KSU_bioinfo_lab/;
my $git_dir = $1;
my $text_out = read_file("${dirname}/Bowtie2_index_template.txt"); ## read shell template with slurp
# remove sequences shorer than 200bp from the reference transcriptome\n";
my (${trans_filename}, ${trans_directories}, ${trans_suffix}) = fileparse($transcriptome,'\..*'); # break appart filenames
my $de_novo_fasta = "${trans_directories}${trans_filename}_gt_200.fasta";
my $index="${trans_directories}${trans_filename}_gt_200";
print SCRIPT eval quote($text_out);
print SCRIPT "\n";
print QSUBS_INDEX "qsub -l mem=10G,h_rt=10:00:00 ${home}/${project_name}_scripts/${project_name}_index.sh\n";
close (SCRIPT);
###############################################################################
############   Open qsub scripts for mapping and cleaning    ##################
###############################################################################
open (QSUBS_CLEAN, '>', "${home}/${project_name}_qsubs/${project_name}_qsubs_clean.sh") or die "Can't open ${home}/${project_name}_qsubs/${project_name}_qsubs_clean.sh!\n";
print QSUBS_CLEAN "#!/bin/bash\n";
open (QSUBS_MAP, '>', "${home}/${project_name}_qsubs/${project_name}_qsubs_map.sh") or die "Can't open ${home}/${project_name}_qsubs/${project_name}_qsubs_map.sh!\n";
print QSUBS_MAP "#!/bin/bash\n";
###############################################################################
##############     Write scripts for each sample             ##################
###############################################################################
for my $samples (@reads)
{
    my @r1 = split(',',$samples->[1]); # get list of forward reads
    my @r2 = split(',',$samples->[2]); # get list of reverse reads
    if (scalar(@r1) != scalar(@r2))
    {
        print "Error the number of forward and reverse read files does not match for sample $samples->[0]!\n"; # each forward file must have a corresponding reverse file
        exit;
    }
    #######################################################################
    ############ Convert headers of illumina paired-end data ##############
    #######################################################################
    for my $file (0..$#r1)
    {
        my (${filename}, ${directories}, ${suffix}) = fileparse($r1[$file],'\..*'); # break appart filenames
        my (${filename2}, ${directories2}, ${suffix2}) = fileparse($r2[$file],'\..*'); # break appart filenames
        $out_dir = ${directories};
        open (SCRIPT, '>', "${home}/${project_name}_scripts/${filename}_clean.sh") or die "Can't open ${home}/${project_name}_scripts/${filename}_clean.sh!\n"; # create a shell script for each read-pair set
        print SCRIPT "#!/bin/bash\n";
        if ($convert_header)
        {
            print SCRIPT "#######################################################################\n############ Convert headers of illumina paired-end data ##############\n#######################################################################\n";
                print SCRIPT "cat $r1[$file] | awk \'{if (NR % 4 == 1) {split(\$1, arr, \":\"); printf \"%s_%s:%s:%s:%s:%s#0/%s\\n\", arr[1], arr[3], arr[4], arr[5], arr[6], arr[7], substr(\$2, 1, 1), \$0} else if (NR % 4 == 3){print \"+\"} else {print \$0} }\' > ${directories}${filename}_header.fastq\n";
                $r1[$file] = "${directories}${filename}_header.fastq";
                print SCRIPT "cat $r2[$file] | awk \'{if (NR % 4 == 1) {split(\$1, arr, \":\"); printf \"%s_%s:%s:%s:%s:%s#0/%s\\n\", arr[1], arr[3], arr[4], arr[5], arr[6], arr[7], substr(\$2, 1, 1), \$0} else if (NR % 4 == 3){print \"+\"} else {print \$0} }\' > ${directories}${filename2}_header.fastq\n";
                $r2[$file] = "${directories}${filename2}_header.fastq";
            
        }
        #######################################################################
        ######### Clean reads for low quality without de-duplicating ##########
        #######################################################################
        print QSUBS_CLEAN "qsub -l h_rt=48:00:00,mem=10G ${home}/${project_name}_scripts/${filename}_clean.sh\n";
        $text_out = read_file("${dirname}/Prinseq_template.txt"); ## read shell template with slurp
        print SCRIPT eval quote($text_out);
        print SCRIPT "\n";
        if ($clean_read_file1)
        {
            $clean_read_file1 = "$clean_read_file1"." ${directories}${filename}_good_1.fastq";
            $clean_read_file2 = "$clean_read_file2"." ${directories}${filename}_good_2.fastq";
            $clean_read_singletons = "$clean_read_singletons". " ${directories}${filename}_good_1_singletons.fastq ${directories}${filename}_good_2_singletons.fastq";
        }
        else
        {
            $clean_read_file1 = " ${directories}${filename}_good_1.fastq";
            $clean_read_file2 = " ${directories}${filename}_good_2.fastq";
            $clean_read_singletons = " ${directories}${filename}_good_1_singletons.fastq ${directories}${filename}_good_2_singletons.fastq";
        }
    }
    #######################################################################
    ######### Map reads to de novo transcriptome using Bowtie2   ##########
    #######################################################################
    close (SCRIPT);
    open (SCRIPT, '>', "${home}/${project_name}_scripts/$samples->[0]_map.sh") or die "Can't open ${home}/${project_name}_scripts/$samples->[0]_map.sh!\n"; # create a shell script for each read-pair set
    $text_out = read_file("${dirname}/Bowtie2_map_template.txt"); ## read shell template with slurp
    print SCRIPT eval quote($text_out);
    print SCRIPT "\n";
    
    if ($labels)
    {
        $sams = "$sams".",${out_dir}$samples->[0]_200.sam";
        $labels  = "$labels".",$samples->[0]";
    }
    else
    {
        $sams = "${out_dir}$samples->[0]_200.sam";
        $labels  = "$samples->[0]";
    }
    print QSUBS_MAP "qsub -l h_rt=48:00:00,mem=2G -pe single 20 ${home}/${project_name}_scripts/$samples->[0]_map.sh\n";
    
}
#######################################################################
#########                Summarize read counts               ##########
#######################################################################
close (SCRIPT);
open (SCRIPT, '>', "${home}/${project_name}_scripts/${project_name}_count.sh") or die "Can't open ${home}/${project_name}_scripts/${project_name}_count.sh!\n"; # create a shell script
${dirname} =~ /(.*RNA-Seq-annotation-and-comparison\/KSU_bioinfo_lab)(\/RNA-SeqAlign)/;
my $shortdirname = $1;
$text_out = read_file("${dirname}/Count_reads_template.txt"); ## read shell template with slurp
print SCRIPT eval quote($text_out);
print SCRIPT "\n";
open (QSUBS_COUNT, '>', "${home}/${project_name}_qsubs/${project_name}_qsubs_count.sh") or die "Can't open ${home}/${project_name}_qsubs/${project_name}_qsubs_count.sh!\n";
print QSUBS_COUNT "#!/bin/bash\n";<commit_msg>debug path (relative to absolute)<commit_after>#!/bin/perl
###############################################################################
#   
#	USAGE: perl RNA-SeqAlign.pl [options]
#
#  Created by <NAME>
#
###############################################################################
use strict;
use warnings;
use File::Basename; # enable manipulating of the full path
use Cwd;
use lib '/homes/bioinfo/bioinfo_software/perl_modules/File-Slurp-9999.19/lib';
use File::Slurp;
# use List::Util qw(max);
# use List::Util qw(sum);
# use Bio::SeqIO;
# use Bio::Seq;
# use Bio::DB::Fasta;
use Getopt::Long;
use Pod::Usage;
###############################################################################
##############         Print informative message             ##################
###############################################################################
print "###########################################################\n";
print "#  RNA-SeqAlign.pl                                        #\n";
print "#                                                         #\n";
print "#  Created by <NAME> 3/15/14                    #\n";
print "#  github.com/i5K-KINBRE-script-share                     #\n";
print "#  perl RNA-SeqAlign.pl -help # for usage/options         #\n";
print "#  perl RNA-SeqAlign.pl -man # for more details           #\n";
print "###########################################################\n";
###############################################################################
##############                get arguments                  ##################
###############################################################################
my ($r_list,$project_name,$transcriptome,$min_len,$clean_read_file1,$clean_read_file2,@clean_r1,@clean_r2,$labels,$sams,$clean_read_singletons);
my $convert_header = 0;
my $mapq = 10;
my $man = 0;
my $help = 0;
GetOptions (
			  'help|?' => \$help, 
			  'man' => \$man,
			  'r|r_list:s' => \$r_list,
			  't|transcriptome:s' => \$transcriptome,
              'p|project_name:s' => \$project_name,
			  'c|convert_header' => \$convert_header,
			  'm|mapq:i' => \$mapq,
			  'l|read_length:i' => \$min_len

              )  
or pod2usage(2);
pod2usage(1) if $help;
pod2usage(-exitstatus => 0, -verbose => 2) if $man;
sub quote { qq!"$_[0]"! } ## interpolate slurped text
my $dirname = dirname(__FILE__); # github directories (all github directories must be in the same directory)
my $home = getcwd; # working directory (this is where output files will be printed)
mkdir "${home}/${project_name}_scripts";
mkdir "${home}/${project_name}_qsubs";
mkdir "${home}/${project_name}_prinseq";
###############################################################################
############## Create array of the sample names and read files    #############
###############################################################################
my @reads;
open (READ_LIST, '<', $r_list) or die "Can't open $r_list!\n";
while (<READ_LIST>)
{
    chomp;
    push @reads , [split];
}
#######################################################################
#########     Index the de novo transcriptome for Bowtie2    ##########
#######################################################################
close (SCRIPT);
open (SCRIPT, '>', "${home}/${project_name}_scripts/${project_name}_index.sh") or die "Can't open ${home}/${project_name}_scripts/${project_name}_index.sh!\n"; # create a shell script for each read-pair set
open (QSUBS_INDEX, '>', "${home}/${project_name}_qsubs/${project_name}_qsubs_index.sh") or die "Can't open ${home}/${project_name}_qsubs/${project_name}_qsubs_index.sh!\n";
print QSUBS_INDEX "#!/bin/bash\n";
$dirname =~ /(.*)\/RNA-Seq-annotation-and-comparison\/KSU_bioinfo_lab/;
my $git_dir = $1;
my $text_out = read_file("${dirname}/Bowtie2_index_template.txt"); ## read shell template with slurp
# remove sequences shorer than 200bp from the reference transcriptome\n";
my (${trans_filename}, ${trans_directories}, ${trans_suffix}) = fileparse($transcriptome,'\..*'); # break appart filenames
my $de_novo_fasta = "${trans_directories}${trans_filename}_gt_200.fasta";
my $index="${trans_directories}${trans_filename}_gt_200";
print SCRIPT eval quote($text_out);
print SCRIPT "\n";
print QSUBS_INDEX "qsub -l mem=10G,h_rt=10:00:00 ${home}/${project_name}_scripts/${project_name}_index.sh\n";
close (SCRIPT);
###############################################################################
############   Open qsub scripts for mapping and cleaning    ##################
###############################################################################
open (QSUBS_CLEAN, '>', "${home}/${project_name}_qsubs/${project_name}_qsubs_clean.sh") or die "Can't open ${home}/${project_name}_qsubs/${project_name}_qsubs_clean.sh!\n";
print QSUBS_CLEAN "#!/bin/bash\n";
open (QSUBS_MAP, '>', "${home}/${project_name}_qsubs/${project_name}_qsubs_map.sh") or die "Can't open ${home}/${project_name}_qsubs/${project_name}_qsubs_map.sh!\n";
print QSUBS_MAP "#!/bin/bash\n";
###############################################################################
##############     Write scripts for each sample             ##################
###############################################################################
for my $samples (@reads)
{
    my @r1 = split(',',$samples->[1]); # get list of forward reads
    my @r2 = split(',',$samples->[2]); # get list of reverse reads
    if (scalar(@r1) != scalar(@r2))
    {
        print "Error the number of forward and reverse read files does not match for sample $samples->[0]!\n"; # each forward file must have a corresponding reverse file
        exit;
    }
    #######################################################################
    ############ Convert headers of illumina paired-end data ##############
    #######################################################################
    for my $file (0..$#r1)
    {
        my (${filename}, ${directories}, ${suffix}) = fileparse($r1[$file],'\..*'); # break appart filenames
        my (${filename2}, ${directories2}, ${suffix2}) = fileparse($r2[$file],'\..*'); # break appart filenames
        open (SCRIPT, '>', "${home}/${project_name}_scripts/${filename}_clean.sh") or die "Can't open ${home}/${project_name}_scripts/${filename}_clean.sh!\n"; # create a shell script for each read-pair set
        print SCRIPT "#!/bin/bash\n";
        if ($convert_header)
        {
            print SCRIPT "#######################################################################\n############ Convert headers of illumina paired-end data ##############\n#######################################################################\n";
                print SCRIPT "cat $r1[$file] | awk \'{if (NR % 4 == 1) {split(\$1, arr, \":\"); printf \"%s_%s:%s:%s:%s:%s#0/%s\\n\", arr[1], arr[3], arr[4], arr[5], arr[6], arr[7], substr(\$2, 1, 1), \$0} else if (NR % 4 == 3){print \"+\"} else {print \$0} }\' > ${directories}${filename}_header.fastq\n";
                $r1[$file] = "${directories}${filename}_header.fastq";
                print SCRIPT "cat $r2[$file] | awk \'{if (NR % 4 == 1) {split(\$1, arr, \":\"); printf \"%s_%s:%s:%s:%s:%s#0/%s\\n\", arr[1], arr[3], arr[4], arr[5], arr[6], arr[7], substr(\$2, 1, 1), \$0} else if (NR % 4 == 3){print \"+\"} else {print \$0} }\' > ${directories}${filename2}_header.fastq\n";
                $r2[$file] = "${directories}${filename2}_header.fastq";
            
        }
        #######################################################################
        ######### Clean reads for low quality without de-duplicating ##########
        #######################################################################
        print QSUBS_CLEAN "qsub -l h_rt=48:00:00,mem=10G ${home}/${project_name}_scripts/${filename}_clean.sh\n";
        $text_out = read_file("${dirname}/Prinseq_template.txt"); ## read shell template with slurp
        print SCRIPT eval quote($text_out);
        print SCRIPT "\n";
        if ($clean_read_file1)
        {
            $clean_read_file1 = "$clean_read_file1"." ${directories}${filename}_good_1.fastq";
            $clean_read_file2 = "$clean_read_file2"." ${directories}${filename}_good_2.fastq";
            $clean_read_singletons = "$clean_read_singletons". " ${directories}${filename}_good_1_singletons.fastq ${directories}${filename}_good_2_singletons.fastq";
        }
        else
        {
            $clean_read_file1 = " ${directories}${filename}_good_1.fastq";
            $clean_read_file2 = " ${directories}${filename}_good_2.fastq";
            $clean_read_singletons = " ${directories}${filename}_good_1_singletons.fastq ${directories}${filename}_good_2_singletons.fastq";
        }
    }
    #######################################################################
    ######### Map reads to de novo transcriptome using Bowtie2   ##########
    #######################################################################
    close (SCRIPT);
    open (SCRIPT, '>', "${home}/${project_name}_scripts/$samples->[0]_map.sh") or die "Can't open ${home}/${project_name}_scripts/$samples->[0]_map.sh!\n"; # create a shell script for each read-pair set
    $text_out = read_file("${dirname}/Bowtie2_map_template.txt"); ## read shell template with slurp
    print SCRIPT eval quote($text_out);
    print SCRIPT "\n";
    
    if ($labels)
    {
        $sams = "$sams".",${home}/$samples->[0]_200.sam";
        $labels  = "$labels".",$samples->[0]";
    }
    else
    {
        $sams = "${home}/$samples->[0]_200.sam";
        $labels  = "$samples->[0]";
    }
    print QSUBS_MAP "qsub -l h_rt=48:00:00,mem=2G -pe single 20 ${home}/${project_name}_scripts/$samples->[0]_map.sh\n";
    
}
#######################################################################
#########                Summarize read counts               ##########
#######################################################################
close (SCRIPT);
open (SCRIPT, '>', "${home}/${project_name}_scripts/${project_name}_count.sh") or die "Can't open ${home}/${project_name}_scripts/${project_name}_count.sh!\n"; # create a shell script
${dirname} =~ /(.*RNA-Seq-annotation-and-comparison\/KSU_bioinfo_lab)(\/RNA-SeqAlign)/;
my $shortdirname = $1;
$text_out = read_file("${dirname}/Count_reads_template.txt"); ## read shell template with slurp
print SCRIPT eval quote($text_out);
print SCRIPT "\n";
open (QSUBS_COUNT, '>', "${home}/${project_name}_qsubs/${project_name}_qsubs_count.sh") or die "Can't open ${home}/${project_name}_qsubs/${project_name}_qsubs_count.sh!\n";
print QSUBS_COUNT "#!/bin/bash\n";
{'fraction_non_alphanumeric': 0.35255805172739285, 'fraction_numerical': 0.014253350953369832, 'mean_word_length': 4.543956043956044, 'pattern_counts': {'":': 4, '<': 11, '<?xml version=': 0, '>': 69, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3423632', 'n_tokens_mistral': 6644, 'n_tokens_neox': 6116, 'n_words': 1275}