#!/usr/bin/perl
use strict;

### Script to map the reads of P.falciparum IT chr5 against the malaria reference genome.

print "\nSawadee krap!\n\nThis script will automatically\n - map sequence reads to the malaria reference genome, \n";
print " - call SNPs, and \n - prepare files for viewing in Artemis.\n\nThe whole procedure will take several minutes, please be patient...\n";

print "\n*** The following command will index the reference sequence ***\n";
my $command = "bwa index Pf3D7_05.fasta";
print "$command\n";
system($command);

print "\n*** The following commands will align forward and reverse reads ***\n";
$command = "bwa aln Pf3D7_05.fasta IT.Chr5_1.fastq > F.sai";
print "$command\n";
system($command);
$command = "bwa aln Pf3D7_05.fasta IT.Chr5_2.fastq > R.sai";
print "\n$command\n";
system($command);

print "\n*** The following command will join the SAI files and generate a SAM file ***\n";
$command = "bwa sampe -a 700 Pf3D7_05.fasta F.sai R.sai IT.Chr5_1.fastq  IT.Chr5_2.fastq > BWA.sam";
print "$command\n";
system($command);

print "\n*** The following command will convert SAM to BAM file format ***\n";
$command = "samtools view -b -S BWA.sam > BWA.bam";
print "$command\n";
system($command);

print "\n*** The following command will sort the mapped reads ***\n";
$command = "samtools sort BWA.bam BWA.IT.Chr5";
print "$command\n";
system($command);

print "\n*** The following command will index the BAM file ***\n";
$command = "samtools index BWA.IT.Chr5.bam";
print "$command\n";
system($command);


print "\n*** The following command will do the SNP calling and generate a BCF file ***\n";
$command = "samtools mpileup -ugf Pf3D7_05.fasta BWA.IT.Chr5.bam  |  bcftools view -bcvg - > IT.Chr5.bcf";
print "$command\n";
system($command);

print "\n*** The following command will index the BCF file ***\n";
$command = "bcftools index IT.Chr5.bcf";
print "$command\n";
system($command);

print "\n*** The following command will start Artemis for you! ***\n";
$command = "art -Dbam=BWA.IT.Chr5.bam,IT.Chr5.bcf Pf3D7_05.embl &";
print "$command\n";
system($command);

print "\nThank you for having used this script, the script is DONE now.\nHave a nice day!\n\n";

