#!/usr/bin/perl -w
#
# Converts bismark mapping output to SAM format.
# Usage: bismark_to_SAM.pl (-p) -c [chrom sizes file] -i [bismark mapping output] -o [SAM output]

use strict;
use warnings;
use Getopt::Std;
use Carp;

### Program variables ###

my $chrom_sizes_file; # Contains chromosome name used in mapping (column 1) and its length (column 2)
my $bismark_infile;
my $sam_outfile;
my $is_paired_end = 0;

### Main program ###

parse_cmd_line();
verify_input_file();
generate_SAM_header();
if($is_paired_end){
  convert_pe_bismark_to_SAM();
}
else{
  convert_se_bismark_to_SAM();
}

### Subroutines ###

sub parse_cmd_line{
  my %opts;
  getopts('c:i:o:p', \%opts);
  $chrom_sizes_file = $opts{'c'};
  warn "no chrom_sizes file indicated\n" unless defined $chrom_sizes_file;
  $bismark_infile = $opts{'i'};
  warn "no input file specified\n" unless defined $bismark_infile;
  $is_paired_end = 1 if exists $opts{'p'};
  if(!defined $chrom_sizes_file || !defined $bismark_infile){
    die usage();
  }
  $sam_outfile = $opts{'o'};
  $sam_outfile = $bismark_infile . ".sam" unless defined $sam_outfile;
}
  
sub usage{
  print <<EOF

Usage:
   bismark_to_SAM.pl (-p) -c [chrom sizes file] -i [bismark mapping output] -o [SAM output]

  -p                           -  paired-end bismark mapping results (default is single-end)
  -c [chrom sizes file]        -  file containing length of chromosomes/sequences used for Bowtie mapping
  -i [bismark mapping output]  -  file containing bismark mapping output
  -o [SAM output]              -  name for output file in SAM format (default: [input].sam)

EOF
    ;
  exit 1;
}

sub verify_input_file{
  open my $ifh, "<", $bismark_infile or die "Cannot open input file ($bismark_infile). $!";
  my $line = <$ifh>;
  $line = <$ifh> if $line =~ /^Bismark version/;
  chomp $line;
  my @NR = split "\t", $line;
  if($is_paired_end){
    die "Input file ($bismark_infile) is not a paired-end Bismark mapping result\n" if scalar @NR !=13;
  }
  else{
    die "Input file ($bismark_infile) is not a single-end Bismark mapping results\n" if scalar @NR !=10;
  }
  close $ifh or die $!;
}

sub generate_SAM_header{
  open my $ifh, "<", $chrom_sizes_file or die "Cannot open chrom_sizes file ($chrom_sizes_file). $!";
  open my $ofh, ">", $sam_outfile or die "Cannot open output file ($sam_outfile). $!";
  print {$ofh} "\@" . "HD\tVN:1.0\tSO:unsorted\n";          # @HD = header, VN = version, SO = sort order [technically sorted by chromosome and start coordinate]
  while(my $line = <$ifh>){
    chomp $line;
    my ($chr, $length, @rest) = split "\t", $line;
    die "Chromosome length is not numeric in file ($chrom_sizes_file) on line $.\n" if $length =~ /\D/;
    print {$ofh} "\@" . "SQ\tSN:$chr\tLN:$length\n";        # @SQ = sequence, SN = seq name, LN = length
  }
}

sub convert_se_bismark_to_SAM{
  open my $ifh, "<", $bismark_infile or die "Cannot open input file ($bismark_infile). $!";
  open my $ofh, "|-", "sort -k3,3 -k4,4n >>$sam_outfile" or die "Cannot sort or append to output file ($sam_outfile). $!";
  while(my $line = <$ifh>){
    next if $line =~ /^Bismark version/;
    chomp $line;
    my ($id, $strand, $chr, $start, $stop, $actual_seq, $ref_seq, @rest) = split "\t", $line;

    ### Input validations ###
    die "Strand information ($strand) is incorrect in file ($bismark_infile) on line $.\n" if $strand =~ /[^\+\-]/;
    die "Start position ($start) is not numeric in file ($bismark_infile) on line $.\n" if $start =~ /\D/;
    die "End position ($stop) is not numeric in file ($bismark_infile) on line $.\n" if $stop =~ /\D/;

    # Assumes bisulfite sequences has A,C,T,G and N only
    die "Bisulfite sequence ($actual_seq) contains invalid nucleotides in file ($bismark_infile) on line $.\n" if $actual_seq =~ /[^ACTGNactgn]/;
    
    # Allows all degenerate nucleotide sequences in reference genome
    die "Reference sequence ($ref_seq) contains invalid nucleotides in file ($bismark_infile) on line $.\n" if $ref_seq =~ /[^ACTGNRYMKSWBDHVactgnrymkswbdhv]/;

    my $flag;                                        # FLAG variable used for SAM format.
    ($strand eq "+") ? ($flag = 0) : ($flag = 16);   # Currently, this only handles single-end mapping, with 0 for "+" strand and 16 for "-" strand
#    $start ++;                                      # Uncomment if using 0-based mapping
    my $mapq = 255;                                  # Assume mapping quality is unavailable
    my $cigar = length($actual_seq) . "M";           # Assume no indel during mapping (only matches and mismatches)
    my $mrnm = "*";                                  # Paired-end variable
    my $mpos = 0;                                    # Paired-end variable
    my $isize = 0;                                   # Paired-end variable
    $actual_seq = revcomp($actual_seq) if $strand eq "-"; # Sequence represented on the forward genomic strand 
    $ref_seq = substr($ref_seq, 0, length($ref_seq) - 2); # Removes additional nucleotide at the end
    $ref_seq = revcomp($ref_seq) if $strand eq "-";       # Required for comparison with actual sequence
    my $qual = "I" x length($ref_seq);               # Assume highest quality, since quality information is unavailable
    my $hemming_dist = hemming_dist($actual_seq, $ref_seq);
    my $NM_tag = "NM:i:$hemming_dist";               # Optional tag: edit distance based on nucleotide differences 
    my $MD_tag = make_mismatch_string($actual_seq, $ref_seq);  # Optional tag: String to provide mismatch and indel information. Expect only mismatches

# SAM format: QNAME, FLAG, RNAME, 1-based START, MAPQ, CIGAR, MRNM, MPOS, ISIZE, SEQ, QUAL
    print {$ofh} join("\t", ($id, $flag, $chr, $start, $mapq, $cigar, $mrnm, $mpos, $isize, $actual_seq, $qual, $NM_tag, $MD_tag)), "\n";
  }
}

sub convert_pe_bismark_to_SAM{
  open my $ifh, "<", $bismark_infile or die "Cannot open input file ($bismark_infile). $!";
  open my $ofh, "|-", "sort -k3,3 -k4,4n >>$sam_outfile" or die "Cannot sort or append to output file ($sam_outfile). $!";
  while(my $line = <$ifh>){
    next if $line =~ /^Bismark version/;
    chomp $line;
    my ($id, $strand, $chr, $start, $stop, $actual_seq1, $ref_seq1, $methcall1, $actual_seq2, $ref_seq2, $methcall2, @rest) = split "\t", $line;

    ### Input validations ###
    die "Strand information ($strand) is incorrect in file ($bismark_infile) on line $.\n" if $strand =~ /[^\+\-]/;
    die "Start position ($start) is not numeric in file ($bismark_infile) on line $.\n" if $start =~ /\D/;
    die "End position ($stop) is not numeric in file ($bismark_infile) on line $.\n" if $stop =~ /\D/;

    # Assumes bisulfite sequences has A,C,T,G and N only
    die "Bisulfite sequence ($actual_seq1) contains invalid nucleotides in file ($bismark_infile) on line $.\n" if $actual_seq1 =~ /[^ACTGNactgn]/;
    die "Bisulfite sequence ($actual_seq2) contains invalid nucleotides in file ($bismark_infile) on line $.\n" if $actual_seq2 =~ /[^ACTGNactgn]/;
    
    # Allows all degenerate nucleotide sequences in reference genome
    die "Reference sequence ($ref_seq1) contains invalid nucleotides in file ($bismark_infile) on line $.\n" if $ref_seq1 =~ /[^ACTGNRYMKSWBDHVactgnrymkswbdhv]/;
    die "Reference sequence ($ref_seq2) contains invalid nucleotides in file ($bismark_infile) on line $.\n" if $ref_seq2 =~ /[^ACTGNRYMKSWBDHVactgnrymkswbdhv]/;

    # Generate output for Read 1
    my $flag1;                                       # FLAG variable used for SAM format.
    ($strand eq "+") ? ($flag1 = 99) : ($flag1 = 83); # 99 if read 1 is positive strand, 83 if read 1 is negative strand
#    $start ++;                                      # Uncomment if using 0-based mapping
    my $read1start;
    ($strand eq "+") ? ($read1start = $start) 
      : ($read1start = $stop - length($actual_seq2));
    my $mapq = 255;                                  # Assume mapping quality is unavailable
    my $cigar1 = length($actual_seq1) . "M";         # Assume no indel during mapping (only matches and mismatches)
    my $mrnm = $chr;                                 # Chromosome of mate
    my $mpos1;
    ($strand eq "+") ? ($mpos1 = $stop - length($actual_seq2))
      : ($mpos1 = $start);                           # Left-most position of mate
    my $isize1 = $mpos1 - $read1start;               # Insert size
    $actual_seq1 = revcomp($actual_seq1) if $strand eq "-"; # Sequence represented on the forward genomic strand 
    $ref_seq1 = substr($ref_seq1, 0, length($ref_seq1) - 2); # Removes additional 2 nucleotides at the end (not sure why it's like this)
    $ref_seq1 = revcomp($ref_seq1) if $strand eq "-";        # Required for comparison with actual sequence
    my $qual1 = "I" x length($ref_seq1);               # Assume highest quality, since quality information is unavailable
    my $hemming_dist = hemming_dist($actual_seq1, $ref_seq1);
    my $NM_tag1 = "NM:i:$hemming_dist";               # Optional tag: edit distance based on nucleotide differences 
    my $MD_tag1 = make_mismatch_string($actual_seq1, $ref_seq1);  # Optional tag: String to provide mismatch and indel information. Expect only mismatches

    # Generate output for Read 2
    my $flag2;                                       # FLAG variable used for SAM format.
    ($strand eq "+") ? ($flag2 = 147) : ($flag2 = 163); # 147 if read 1 is positive strand, 163 if read 1 is negative strand
    my $read2start;
    ($strand eq "-") ? ($read2start = $start) 
      : ($read2start = $stop - length($actual_seq2));
    my $cigar2 = length($actual_seq2) . "M";         # Assume no indel during mapping (only matches and mismatches)
    my $mpos2;
    ($strand eq "-") ?  ($mpos2 = $stop - length($actual_seq2))
      : ($mpos2 = $start);                           # Left-most position of mate
    my $isize2 = $mpos2 - $read2start;                    # Insert size
    $actual_seq2 = revcomp($actual_seq2) if $strand eq "+"; # Mate sequence represented on the forward genomic strand 
    $ref_seq2 = substr($ref_seq2, 0, length($ref_seq2) - 2); # Removes additional 2 nucleotides at the end (not sure why it's like this)
    $ref_seq2 = revcomp($ref_seq2) if $strand eq "+";        # Required for comparison with actual sequence
    my $qual2 = "I" x length($ref_seq2);              # Assume highest quality, since quality information is unavailable
    $hemming_dist = hemming_dist($actual_seq2, $ref_seq2);
    my $NM_tag2 = "NM:i:$hemming_dist";               # Optional tag: edit distance based on nucleotide differences 
    my $MD_tag2 = make_mismatch_string($actual_seq2, $ref_seq2);  # Optional tag: String to provide mismatch and indel information. Expect only mismatches

# SAM format: QNAME, FLAG, RNAME, 1-based START, MAPQ, CIGAR, MRNM, MPOS, ISIZE, SEQ, QUAL
    print {$ofh} join("\t", ($id, $flag1, $chr, $read1start, $mapq, $cigar1, $mrnm, $mpos1, $isize1, $actual_seq1, $qual1, $NM_tag1, $MD_tag1)), "\n";
    print {$ofh} join("\t", ($id, $flag2, $chr, $read2start, $mapq, $cigar2, $mrnm, $mpos2, $isize2, $actual_seq2, $qual2, $NM_tag2, $MD_tag2)), "\n";
  }
}

sub revcomp{
  my $seq = shift or croak "Missing seq to reverse complement";
  $seq = reverse $seq;
  $seq =~ tr/ACTGactg/TGACtgac/;
  return $seq;
}

sub hemming_dist{
  my $string1 = shift or croak "Missing string 1";
  my $string2 = shift or croak "Missing string 2";
  return my $hd = length($string1) - (($string1 ^ $string2) =~ tr/[\0]/[\0]/);
}

sub make_mismatch_string{
  my $actual_seq = shift or croak "Missing actual sequence";
  my $ref_seq = shift or croak "Missing reference sequence";
  my $MD_tag = "MD:Z:";
  my $tmp = ($actual_seq ^ $ref_seq);              # Bitwise comparison
  my $prev_mm_pos = 0;
  while($tmp =~ /[^\0]/g){                         # Where bitwise comparison showed a difference
    my $nuc_match = pos($tmp) - $prev_mm_pos - 1;  # Generate number of nucleotide that matches since last mismatch
    my $nuc_mm = substr($actual_seq, pos($tmp) - 1, 1) if pos($tmp) <= length($actual_seq);  # Obtain nucleotide that was different from reference
    $MD_tag .= "$nuc_match" if $nuc_match > 0;     # Ignore if mismatches are adjacent to each other
    $MD_tag .= "$nuc_mm" if defined $nuc_mm;       # Ignore if there is no mismatch (prevents uninitialized string concatenation)
    $prev_mm_pos = pos($tmp);                      # Position of last mismatch
  }
  my $end_matches = length($actual_seq) - $prev_mm_pos;  # Provides number of matches from last mismatch till end of sequence
  $MD_tag .= "$end_matches" if $end_matches > 0;   # Ignore if mismatch is at the end of sequence
  return $MD_tag;
}


