misc/export2sam.pl

   1 #!/usr/bin/env perl\r
   2 #\r
   3 #\r
   4 # export2sam.pl converts GERALD export files to SAM format.\r
   5 #\r
   6 #\r
   7 #\r
   8 ########## License:\r
   9 #\r
  10 # The MIT License\r
  11 #\r
  12 # Original SAMtools work copyright (c) 2008-2009 Genome Research Ltd.\r
  13 # Modified SAMtools work copyright (c) 2010 Illumina, Inc.\r
  14 #\r
  15 # Permission is hereby granted, free of charge, to any person obtaining a copy\r
  16 # of this software and associated documentation files (the "Software"), to deal\r
  17 # in the Software without restriction, including without limitation the rights\r
  18 # to use, copy, modify, merge, publish, distribute, sublicense, and/or sell\r
  19 # copies of the Software, and to permit persons to whom the Software is\r
  20 # furnished to do so, subject to the following conditions:\r
  21 #\r
  22 # The above copyright notice and this permission notice shall be included in\r
  23 # all copies or substantial portions of the Software.\r
  24 # \r
  25 # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR\r
  26 # IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,\r
  27 # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE\r
  28 # AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER\r
  29 # LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,\r
  30 # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN\r
  31 # THE SOFTWARE.\r
  32 #\r
  33 #\r
  34 #\r
  35 #\r
  36 ########## ChangeLog:\r
  37 #\r
  38 # Version: 2.3.1 (18MAR2011)\r
  39 #\r
  40 #   - Restore file '-' as stdin input.\r
  41 #\r
  42 # Version: 2.3.0 (24JAN2011)\r
  43 #\r
  44 #   - Add support for export reserved chromosome name "CONTROL",\r
  45 #       which is translated to optional field "XC:Z:CONTROL".\r
  46 #   - Check for ".gz" file extension on export files and open\r
  47 #       these as gzip pipes when the extension is found.\r
  48 #\r
  49 # Version: 2.2.0 (16NOV2010)\r
  50 #\r
  51 #   - Remove any leading zeros in export fields: RUNNO,LANE,TILE,X,Y\r
  52 #   - For export records with reserved chromosome name identifiers\r
  53 #       "QC" and "RM", add the optional field "XC:Z:QC" or "XC:Z:RM"\r
  54 #       to the SAM record, so that these cases can be distinguished\r
  55 #       from other unmatched reads.\r
  56 #\r
  57 # Version: 2.1.0 (21SEP2010)\r
  58 #\r
  59 #   - Additional export record error checking.\r
  60 #   - Convert export records with chromomsome value of "RM" to unmapped \r
  61 #       SAM records.\r
  62 #\r
  63 # Version: 2.0.0 (15FEB2010)\r
  64 #\r
  65 #   Script updated by Illumina in conjunction with CASAVA 1.7.0\r
  66 #   release.\r
  67 #\r
  68 #   Major changes are as follows:\r
  69 #   - The CIGAR string has been updated to include all gaps from\r
  70 #       ELANDv2 alignments.\r
  71 #   - The ELAND single read alignment score is always stored in the\r
  72 #       optional "SM" field and the ELAND paired read alignment score\r
  73 #       is stored in the optional "AS" field when it exists.\r
  74 #   - The MAPQ value is set to the higher of the two alignment scores,\r
  75 #       but no greater than 254, i.e. min(254,max(SM,AS))\r
  76 #   - The SAM "proper pair" bit (0x0002) is now set for read pairs\r
  77 #       meeting ELAND's expected orientation and insert size criteria.\r
  78 #   - The default quality score translation is set for export files\r
  79 #       which contain Phread+64 quality values. An option,\r
  80 #       "--qlogodds", has been added to translate quality values from\r
  81 #       the Solexa+64 format used in export files prior to Pipeline\r
  82 #       1.3\r
  83 #   - The export match descriptor is now reverse-complemented when\r
  84 #       necessary such that it always corresponds to the forward\r
  85 #       strand of the reference, to be consistent with other\r
  86 #       information in the SAM record. It is now written to the\r
  87 #       optional 'XD' field (rather than 'MD') to acknowledge its\r
  88 #       minor differences from the samtools match descriptor (see\r
  89 #       additional detail below).\r
  90 #   - An option, "--nofilter", has been added to include reads which\r
  91 #       have failed primary analysis quality filtration. Such reads\r
  92 #       will have the corresponding SAM flag bit (0x0200) set.\r
  93 #   - Labels in the export 'contig' field are preserved by setting\r
  94 #       RNAME to "$export_chromosome/$export_contig" when the contig\r
  95 #       label exists.\r
  96 #\r
  97 #\r
  98 # Contact: lh3\r
  99 # Version: 0.1.2 (03JAN2009)\r
 100 #\r
 101 #\r
 102 #\r
 103 ########## Known Conversion Limitations:\r
 104 #\r
 105 # - Export records for reads that map to a position < 1 (allowed\r
 106 #     in export format), are converted to unmapped reads in the SAM\r
 107 #     record.\r
 108 # - Export records contain the reserved chromosome names: "NM",\r
 109 #     "QC","RM" and "CONTROL". "NM" indicates that the aligner could\r
 110 #     not map the read to the reference sequence set. "QC" means that\r
 111 #     the aligner did not attempt to map the read due to some\r
 112 #     technical limitation. "RM" means that the read mapped to a set\r
 113 #     of 'contaminant' sequences specified in GERALD's RNA-seq\r
 114 #     workflow. "CONTROL" means that the read is a control. All of\r
 115 #     these alignment types are collapsed to the single unmapped\r
 116 #     alignment state in the SAM record, but the optional SAM "XC"\r
 117 #     field is used to record the original reserved chromosome name of\r
 118 #     the read for all but the "NM" case.\r
 119 # - The export match descriptor is slightly different than the\r
 120 #     samtools match descriptor. For this reason it is stored in the\r
 121 #     optional SAM field 'XD' (and not 'MD'). Note that the export\r
 122 #     match descriptor differs from the samtools version in two\r
 123 #     respects: (1) indels are explicitly closed with the '$'\r
 124 #     character and (2) insertions must be enumerated in the match\r
 125 #     descriptor. For example a 35-base read with a two-base insertion\r
 126 #     is described as: 20^2$14\r
 127 #\r
 128 #\r
 129 #\r
 130 \r
 131 my $version = "2.3.1";\r
 132 \r
 133 use strict;\r
 134 use warnings;\r
 135 \r
 136 use Getopt::Long;\r
 137 use File::Spec;\r
 138 use List::Util qw(min max);\r
 139 \r
 140 \r
 141 use constant {\r
 142   EXPORT_MACHINE => 0,\r
 143   EXPORT_RUNNO => 1,\r
 144   EXPORT_LANE => 2,\r
 145   EXPORT_TILE => 3,\r
 146   EXPORT_X => 4,\r
 147   EXPORT_Y => 5,\r
 148   EXPORT_INDEX => 6,\r
 149   EXPORT_READNO => 7,\r
 150   EXPORT_READ => 8,\r
 151   EXPORT_QUAL => 9,\r
 152   EXPORT_CHROM => 10,\r
 153   EXPORT_CONTIG => 11,\r
 154   EXPORT_POS => 12,\r
 155   EXPORT_STRAND => 13, \r
 156   EXPORT_MD => 14,\r
 157   EXPORT_SEMAP => 15,\r
 158   EXPORT_PEMAP => 16,\r
 159   EXPORT_PASSFILT => 21,\r
 160   EXPORT_SIZE => 22,\r
 161 };\r
 162 \r
 163 \r
 164 use constant {\r
 165   SAM_QNAME => 0,\r
 166   SAM_FLAG => 1,\r
 167   SAM_RNAME => 2,\r
 168   SAM_POS => 3,\r
 169   SAM_MAPQ => 4,\r
 170   SAM_CIGAR => 5,\r
 171   SAM_MRNM => 6,\r
 172   SAM_MPOS => 7,\r
 173   SAM_ISIZE => 8,\r
 174   SAM_SEQ => 9,\r
 175   SAM_QUAL => 10,\r
 176 };\r
 177 \r
 178 \r
 179 # function prototypes for Richard's code\r
 180 sub match_desc_to_cigar($);\r
 181 sub match_desc_frag_length($);\r
 182 sub reverse_compl_match_descriptor($);\r
 183 sub write_header($;$;$);\r
 184 \r
 185 \r
 186 &export2sam;\r
 187 exit;\r
 188 \r
 189 \r
 190 \r
 191 \r
 192 sub export2sam {\r
 193 \r
 194   my $cmdline = $0 . " " . join(" ",@ARGV);\r
 195   my $arg_count = scalar @ARGV;\r
 196   my $progname = (File::Spec->splitpath($0))[2];\r
 197 \r
 198   my $is_logodds_qvals = 0; # if true, assume files contain logodds (i.e. "solexa") quality values\r
 199   my $is_nofilter = 0;\r
 200   my $read1file;\r
 201   my $read2file;\r
 202   my $print_version = 0;\r
 203   my $help = 0;\r
 204 \r
 205   my $result = GetOptions( "qlogodds" => \$is_logodds_qvals, \r
 206                            "nofilter" => \$is_nofilter,\r
 207                            "read1=s"  => \$read1file,\r
 208                            "read2=s"  => \$read2file,\r
 209                            "version"  => \$print_version,\r
 210                            "help"     => \$help );\r
 211 \r
 212   my $usage = <<END;\r
 213 \r
 214 $progname converts GERALD export files to SAM format.\r
 215 \r
 216 Usage: $progname --read1=FILENAME [ options ] | --version | --help\r
 217 \r
 218   --read1=FILENAME  read1 export file or '-' for stdin (mandatory)\r
 219                       (file may be gzipped with ".gz" extension)\r
 220   --read2=FILENAME  read2 export file or '-' for stdin\r
 221                       (file may be gzipped with ".gz" extension)\r
 222   --nofilter        include reads that failed the basecaller\r
 223                       purity filter\r
 224   --qlogodds        assume export file(s) use logodds quality values\r
 225                       as reported by OLB (Pipeline) prior to v1.3\r
 226                       (default: phred quality values)\r
 227 \r
 228 END\r
 229 \r
 230   my $version_msg = <<END;\r
 231 \r
 232 $progname version: $version\r
 233 \r
 234 END\r
 235 \r
 236   if((not $result) or $help or ($arg_count==0)) {\r
 237     die($usage);\r
 238   }\r
 239 \r
 240   if(@ARGV) {\r
 241     print STDERR "\nERROR: Unrecognized arguments: " . join(" ",@ARGV) . "\n\n";\r
 242     die($usage);\r
 243   }\r
 244 \r
 245   if($print_version) {\r
 246     die($version_msg);\r
 247   }\r
 248 \r
 249   if(not defined($read1file)) {\r
 250     print STDERR "\nERROR: read1 export file must be specified\n\n";\r
 251     die($usage);\r
 252   }\r
 253 \r
 254   unless((-f $read1file) or ($read1file eq '-')) {\r
 255     die("\nERROR: Can't find read1 export file: '$read1file'\n\n");\r
 256   }\r
 257 \r
 258   if (defined $read2file) {\r
 259     unless((-f $read2file) or ($read2file eq '-')) {\r
 260       die("\nERROR: Can't find read2 export file: '$read2file'\n\n");\r
 261     }\r
 262     if($read1file eq $read2file) {\r
 263       die("\nERROR: read1 and read2 export filenames are the same: '$read1file'\n\n");\r
 264     }\r
 265   }\r
 266 \r
 267   my ($fh1, $fh2, $is_paired);\r
 268 \r
 269   my $read1cmd="$read1file";\r
 270   $read1cmd = "gzip -dc $read1file |" if($read1file =~ /\.gz$/);\r
 271   open($fh1, $read1cmd)\r
 272       or die("\nERROR: Can't open read1 process: '$read1cmd'\n\n");\r
 273   $is_paired = defined $read2file;\r
 274   if ($is_paired) {\r
 275     my $read2cmd="$read2file";\r
 276     $read2cmd = "gzip -dc $read2file |" if($read2file =~ /\.gz$/);\r
 277     open($fh2, $read2cmd)\r
 278         or die("\nERROR: Can't open read2 process: '$read2cmd'\n\n");\r
 279   }\r
 280   # quality value conversion table\r
 281   my @conv_table;\r
 282   if($is_logodds_qvals){ # convert from solexa+64 quality values (pipeline pre-v1.3):\r
 283     for (-64..64) {\r
 284       $conv_table[$_+64] = int(33 + 10*log(1+10**($_/10.0))/log(10)+.499);\r
 285     }\r
 286   } else {               # convert from phred+64 quality values (pipeline v1.3+):\r
 287     for (-64..-1) {\r
 288       $conv_table[$_+64] = undef;\r
 289     }\r
 290     for (0..64) {\r
 291       $conv_table[$_+64] = int(33 + $_);\r
 292     }\r
 293   }\r
 294   # write the header\r
 295   print write_header( $progname, $version, $cmdline );\r
 296   # core loop\r
 297   my $export_line_count = 0;\r
 298   while (<$fh1>) {\r
 299     $export_line_count++;\r
 300     my (@s1, @s2);\r
 301     &export2sam_aux($_, $export_line_count, \@s1, \@conv_table, $is_paired, 1, $is_nofilter);\r
 302     if ($is_paired) {\r
 303       my $read2line = <$fh2>;\r
 304       if(not $read2line){\r
 305         die("\nERROR: read1 and read2 export files do not contain the same number of reads.\n  Extra reads observed in read1 file at line no: $export_line_count.\n\n");\r
 306       }\r
 307       &export2sam_aux($read2line, $export_line_count, \@s2, \@conv_table, $is_paired, 2, $is_nofilter);\r
 308 \r
 309       if (@s1 && @s2) { # then set mate coordinate\r
 310         if($s1[SAM_QNAME] ne $s2[SAM_QNAME]){\r
 311           die("\nERROR: Non-paired reads in export files on line: $export_line_count.\n  Read1: $_  Read2: $read2line\n");\r
 312         }\r
 313 \r
 314         my $isize = 0;\r
 315         if ($s1[SAM_RNAME] ne '*' && $s1[SAM_RNAME] eq $s2[SAM_RNAME]) { # then calculate $isize\r
 316           my $x1 = ($s1[SAM_FLAG] & 0x10)? $s1[SAM_POS] + length($s1[SAM_SEQ]) : $s1[SAM_POS];\r
 317           my $x2 = ($s2[SAM_FLAG] & 0x10)? $s2[SAM_POS] + length($s2[SAM_SEQ]) : $s2[SAM_POS];\r
 318           $isize = $x2 - $x1;\r
 319         }\r
 320 \r
 321         foreach ([\@s1,\@s2,$isize],[\@s2,\@s1,-$isize]){ \r
 322           my ($sa,$sb,$is) = @{$_};\r
 323           if ($sb->[SAM_RNAME] ne '*') {\r
 324             $sa->[SAM_MRNM] = ($sb->[SAM_RNAME] eq $sa->[SAM_RNAME]) ? "=" : $sb->[SAM_RNAME];\r
 325             $sa->[SAM_MPOS] = $sb->[SAM_POS];\r
 326             $sa->[SAM_ISIZE] = $is;\r
 327             $sa->[SAM_FLAG] |= 0x20 if ($sb->[SAM_FLAG] & 0x10);\r
 328           } else {\r
 329             $sa->[SAM_FLAG] |= 0x8;\r
 330           }\r
 331         } \r
 332       }\r
 333     }\r
 334     print join("\t", @s1), "\n" if (@s1);\r
 335     print join("\t", @s2), "\n" if (@s2 && $is_paired);\r
 336   }\r
 337   close($fh1);\r
 338   if($is_paired) {\r
 339     while(my $read2line = <$fh2>){\r
 340       $export_line_count++;\r
 341       die("\nERROR: read1 and read2 export files do not contain the same number of reads.\n  Extra reads observed in read2 file at line no: $export_line_count.\n\n");\r
 342     }\r
 343     close($fh2);\r
 344   }\r
 345 }\r
 346 \r
 347 sub export2sam_aux {\r
 348   my ($line, $line_no, $s, $ct, $is_paired, $read_no, $is_nofilter) = @_;\r
 349   chomp($line);\r
 350   my @t = split("\t", $line);\r
 351   if(scalar(@t) < EXPORT_SIZE) {\r
 352     my $msg="\nERROR: Unexpected number of fields in export record on line $line_no of read$read_no export file. Found " . scalar(@t) . " fields but expected " . EXPORT_SIZE . ".\n";\r
 353     $msg.="\t...erroneous export record:\n" . $line . "\n\n";\r
 354     die($msg);\r
 355   }\r
 356   @$s = ();\r
 357   my $isPassFilt = ($t[EXPORT_PASSFILT] eq 'Y');\r
 358   return if(not ($isPassFilt or $is_nofilter));\r
 359   # read name\r
 360   my $samQnamePrefix = $t[EXPORT_MACHINE] . (($t[EXPORT_RUNNO] ne "") ? "_" .  int($t[EXPORT_RUNNO]) : "");\r
 361   $s->[SAM_QNAME] = join(':', $samQnamePrefix, int($t[EXPORT_LANE]), int($t[EXPORT_TILE]),\r
 362                          int($t[EXPORT_X]), int($t[EXPORT_Y]));\r
 363   # initial flag (will be updated later)\r
 364   $s->[SAM_FLAG] = 0;\r
 365   if($is_paired) {\r
 366     if($t[EXPORT_READNO] != $read_no){\r
 367       die("\nERROR: read$read_no export file contains record with read number: " .$t[EXPORT_READNO] . " on line: $line_no\n\n");\r
 368     }\r
 369     $s->[SAM_FLAG] |= 1 | 1<<(5 + $read_no);\r
 370   }\r
 371   $s->[SAM_FLAG] |= 0x200 if (not $isPassFilt);\r
 372 \r
 373   # read & quality\r
 374   my $is_export_rev = ($t[EXPORT_STRAND] eq 'R');\r
 375   if ($is_export_rev) { # then reverse the sequence and quality\r
 376     $s->[SAM_SEQ] = reverse($t[EXPORT_READ]);\r
 377     $s->[SAM_SEQ] =~ tr/ACGTacgt/TGCAtgca/;\r
 378     $s->[SAM_QUAL] = reverse($t[EXPORT_QUAL]);\r
 379   } else {\r
 380     $s->[SAM_SEQ] = $t[EXPORT_READ];\r
 381     $s->[SAM_QUAL] = $t[EXPORT_QUAL];\r
 382   }\r
 383   my @convqual = ();\r
 384   foreach (unpack('C*', $s->[SAM_QUAL])){\r
 385     my $val=$ct->[$_];\r
 386     if(not defined $val){\r
 387       my $msg="\nERROR: can't interpret export quality value: " . $_ . " in read$read_no export file, line: $line_no\n";\r
 388       if( $_ < 64 ) { $msg .= "  Use --qlogodds flag to translate logodds (solexa) quality values.\n"; }\r
 389       die($msg . "\n");\r
 390     }\r
 391     push @convqual,$val;\r
 392   }\r
 393 \r
 394   $s->[SAM_QUAL] = pack('C*',@convqual); # change coding\r
 395 \r
 396 \r
 397   # coor\r
 398   my $has_coor = 0;\r
 399   $s->[SAM_RNAME] = "*";\r
 400   if (($t[EXPORT_CHROM] eq 'NM') or\r
 401       ($t[EXPORT_CHROM] eq 'QC') or\r
 402       ($t[EXPORT_CHROM] eq 'RM') or\r
 403       ($t[EXPORT_CHROM] eq 'CONTROL')) {\r
 404     $s->[SAM_FLAG] |= 0x4; # unmapped\r
 405     push(@$s,"XC:Z:".$t[EXPORT_CHROM]) if($t[EXPORT_CHROM] ne 'NM');\r
 406   } elsif ($t[EXPORT_CHROM] =~ /(\d+):(\d+):(\d+)/) {\r
 407     $s->[SAM_FLAG] |= 0x4; # TODO: should I set BAM_FUNMAP in this case?\r
 408     push(@$s, "H0:i:$1", "H1:i:$2", "H2:i:$3")\r
 409   } elsif ($t[EXPORT_POS] < 1) {\r
 410     $s->[SAM_FLAG] |= 0x4; # unmapped\r
 411   } else {\r
 412     $s->[SAM_RNAME] = $t[EXPORT_CHROM];\r
 413     $s->[SAM_RNAME] .= "/" . $t[EXPORT_CONTIG] if($t[EXPORT_CONTIG] ne '');\r
 414     $has_coor = 1;\r
 415   }\r
 416   $s->[SAM_POS] = $has_coor? $t[EXPORT_POS] : 0;\r
 417 \r
 418 #  print STDERR "t[14] = " . $t[14] . "\n";\r
 419   my $matchDesc = '';\r
 420   $s->[SAM_CIGAR] = "*";\r
 421   if($has_coor){\r
 422     $matchDesc = ($is_export_rev) ? reverse_compl_match_descriptor($t[EXPORT_MD]) : $t[EXPORT_MD];\r
 423 \r
 424     if($matchDesc =~ /\^/){\r
 425       # construct CIGAR string using Richard's function\r
 426       $s->[SAM_CIGAR] = match_desc_to_cigar($matchDesc); # indel processing\r
 427     } else {\r
 428       $s->[SAM_CIGAR] = length($s->[SAM_SEQ]) . "M";\r
 429     }\r
 430   }\r
 431 \r
 432 #  print STDERR "cigar_string = $cigar_string\n";\r
 433 \r
 434   $s->[SAM_FLAG] |= 0x10 if ($has_coor && $is_export_rev);\r
 435   if($has_coor){\r
 436     my $semap = ($t[EXPORT_SEMAP] ne '') ? $t[EXPORT_SEMAP] : 0;\r
 437     my $pemap = 0;\r
 438     if($is_paired) {\r
 439       $pemap = ($t[EXPORT_PEMAP] ne '') ? $t[EXPORT_PEMAP] : 0;\r
 440 \r
 441       # set `proper pair' bit if non-blank, non-zero PE alignment score:\r
 442       $s->[SAM_FLAG] |= 0x02 if ($pemap > 0);\r
 443     }\r
 444     $s->[SAM_MAPQ] = min(254,max($semap,$pemap));\r
 445   } else {\r
 446     $s->[SAM_MAPQ] = 0;\r
 447   }\r
 448   # mate coordinate\r
 449   $s->[SAM_MRNM] = '*';\r
 450   $s->[SAM_MPOS] = 0;\r
 451   $s->[SAM_ISIZE] = 0;\r
 452   # aux\r
 453   push(@$s, "BC:Z:$t[EXPORT_INDEX]") if ($t[EXPORT_INDEX]);\r
 454   if($has_coor){\r
 455     # The export match descriptor differs slightly from the samtools match descriptor.\r
 456     # In order for the converted SAM files to be as compliant as possible,\r
 457     # we put the export match descriptor in optional field 'XD' rather than 'MD':\r
 458     push(@$s, "XD:Z:$matchDesc");\r
 459     push(@$s, "SM:i:$t[EXPORT_SEMAP]") if ($t[EXPORT_SEMAP] ne '');\r
 460     push(@$s, "AS:i:$t[EXPORT_PEMAP]") if ($is_paired and ($t[EXPORT_PEMAP] ne ''));\r
 461   }\r
 462 }\r
 463 \r
 464 \r
 465 \r
 466 #\r
 467 # the following code is taken from Richard Shaw's sorted2sam.pl file\r
 468 #\r
 469 sub reverse_compl_match_descriptor($)\r
 470 {\r
 471 #    print "\nREVERSING THE MATCH DESCRIPTOR!\n";\r
 472     my ($match_desc) = @_;\r
 473     my $rev_compl_match_desc = reverse($match_desc);\r
 474     $rev_compl_match_desc =~ tr/ACGT\^\$/TGCA\$\^/;\r
 475 \r
 476     # Unreverse the digits of numbers.\r
 477     $rev_compl_match_desc = join('',\r
 478                                  map {($_ =~ /\d+/)\r
 479                                       ? join('', reverse(split('', $_)))\r
 480                                       : $_} split(/(\d+)/,\r
 481                                                   $rev_compl_match_desc));\r
 482 \r
 483     return $rev_compl_match_desc;\r
 484 }\r
 485 \r
 486 \r
 487 \r
 488 sub match_desc_to_cigar($)\r
 489 {\r
 490     my ($match_desc) = @_;\r
 491 \r
 492     my @match_desc_parts = split(/(\^.*?\$)/, $match_desc);\r
 493     my $cigar_str = '';\r
 494     my $cigar_del_ch = 'D';\r
 495     my $cigar_ins_ch = 'I';\r
 496     my $cigar_match_ch = 'M';\r
 497 \r
 498     foreach my $match_desc_part (@match_desc_parts) {\r
 499         next if (!$match_desc_part);\r
 500 \r
 501         if ($match_desc_part =~ /^\^([ACGTN]+)\$$/) {\r
 502             # Deletion\r
 503             $cigar_str .= (length($1) . $cigar_del_ch);\r
 504         } elsif ($match_desc_part =~ /^\^(\d+)\$$/) {\r
 505             # Insertion\r
 506             $cigar_str .= ($1 . $cigar_ins_ch);\r
 507         } else {\r
 508             $cigar_str .= (match_desc_frag_length($match_desc_part)\r
 509                            . $cigar_match_ch);\r
 510         }\r
 511     }\r
 512 \r
 513     return $cigar_str;\r
 514 }\r
 515 \r
 516 \r
 517 #------------------------------------------------------------------------------\r
 518 \r
 519 sub match_desc_frag_length($)\r
 520                            {\r
 521     my ($match_desc_str) = @_;\r
 522     my $len = 0;\r
 523 \r
 524     my @match_desc_fields = split(/([ACGTN]+)/, $match_desc_str);\r
 525 \r
 526     foreach my $match_desc_field (@match_desc_fields) {\r
 527         next if ($match_desc_field eq '');\r
 528 \r
 529         $len += (($match_desc_field =~ /(\d+)/)\r
 530                  ? $1 : length($match_desc_field));\r
 531     }\r
 532 \r
 533     return $len;\r
 534 }\r
 535 \r
 536 \r
 537 # argument holds the command line\r
 538 sub write_header($;$;$) \r
 539 {\r
 540         my ($progname,$version,$cl) = @_;\r
 541         my $complete_header = "";\r
 542         $complete_header .= "\@PG\tID:$progname\tVN:$version\tCL:$cl\n";\r
 543 \r
 544         return $complete_header;\r
 545 }\r