#!/usr/local/bin/perl

## This script works if your are working remotely with CBSU machines.

# Adjust the two variables below:
$INDEX="/workdir/nisha/reference_V7/cassavaV7_chrAndScaffoldsCombined_numeric.fa"; # provide genome bwa index prefix (with path)
$threads=30;  # number of threads

$cmd = 'ls *.fastq.gz >fastqList.txt';  ## this has to run through fastq files directory
system($cmd);

open (FIL, "fastqList.txt");
while (<FIL>) {
    chomp;
    if($_ =~ /_R1.fastq.gz/) {
        @tmp = split (/_/, $_); #Files are named in the SM_ID1_ID2_otherinfo_1.fq.gz or SM_ID1_ID2_otherinfo_2.fq.gz
        
        $forwardFile = $tmp[0]."_".$tmp[1]."_".$tmp[2]."_".$tmp[3]."_".$tmp[4]."_R1.fastq.gz";
        $reverseFile = $tmp[0]."_".$tmp[1]."_".$tmp[2]."_".$tmp[3]."_".$tmp[4]."_R2.fastq.gz";
        
        #Sample name
        $sample = "$tmp[4]";
        $sam = $sample.".sam";
        $t = localtime;

        #Set Read Group
        $sm = "$tmp[4]";
        $rg = "$tmp[2]";
        $lb = "lib1";
        $pl = "ILLUMINA";
        
        print "\n\nProcessing fastq files for the sample -  $sample \n";
        print "Reading\n", "\t$forwardFile\n", "\t$reverseFile\n";
	
	# Aligning with BWA-MEM
        $unSortedBAM = "$sample"."_unsorted.bam";
        $cmd = "bwa mem -t $threads -R '\@RG\\tID:$rg\\tSM:$sm\\tLB:$lb\\tPL:$pl' -w 3 -K 10000000 -c 10000";
        $cmd .= " $INDEX $forwardFile $reverseFile | samtools view -Sb1 - > $unSortedBAM";
        print("\t\t$cmd\n");
        system($cmd);
        
        print "\tFinished alignement of $sample at $t\n";

        print "\tSorting BAM file started at $t\n";
        
        # Generating sorted BAM file
$bamFile = $sample.".bam";
        $cmd = "samtools sort -@ $threads -m 4G $unSortedBAM -o $bamFile";
        print("\t\t$cmd\n");
        system($cmd);
        
        print "\tFinished at $t\n";
# removing fastq and unsortedBAM files
        #$cmd = "rm $sam";
        #print("\t$cmd\n");
        #system($cmd);
        
        $cmd = "rm $unSortedBAM";
        print("\t$cmd\n");
        system($cmd);

	$cmd = "rm $forwardFile";
	print("\t$cmd\n");
	system($cmd);

	$cmd = "rm $reverseFile";
	print("\t$cmd\n");
	system($cmd);
        
        # indexing BAM file
        $cmd = "samtools index $bamFile";
        print("\t\t$cmd\n");
        system($cmd);
print "Finished BAM file generation of the sample - $sample\n";
        
               }
}