Mercurial > repos > rliterman > csp2

annotate CSP2/CSP2_env/env-d9b9114564458d9d-741b3de822f2aaca6c6caa4325c4afce/opt/mummer-3.23/dnadiff @ 69:33d812a61356

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload commit 2e9511a184a1ca667c7be0c6321a36dc4e3d116d
author jpayne
date Tue, 18 Mar 2025 17:55:14 -0400
parents
children
rev   line source
jpayne@69 1 #!/usr/bin/env perl
jpayne@69 2
jpayne@69 3 #-------------------------------------------------------------------------------
jpayne@69 4 # Programmer: Adam M Phillippy, University of Maryland
jpayne@69 5 # File: dnadiff
jpayne@69 6 # Date: 11 / 29 / 06
jpayne@69 7 #
jpayne@69 8 # Try 'dnadiff -h' for more information.
jpayne@69 9 #
jpayne@69 10 #-------------------------------------------------------------------------------
jpayne@69 11
jpayne@69 12 use lib "/mnt/c/Users/crash/Documents/BobLiterman/CSP2_Galaxy/CSP2/CSP2_env/env-d9b9114564458d9d-741b3de822f2aaca6c6caa4325c4afce/opt/mummer-3.23/scripts";
jpayne@69 13 use Foundation;
jpayne@69 14 use File::Spec::Functions;
jpayne@69 15 use strict;
jpayne@69 16
jpayne@69 17 my $BIN_DIR = "/mnt/c/Users/crash/Documents/BobLiterman/CSP2_Galaxy/CSP2/CSP2_env/env-d9b9114564458d9d-741b3de822f2aaca6c6caa4325c4afce/opt/mummer-3.23";
jpayne@69 18 my $SCRIPT_DIR = "/mnt/c/Users/crash/Documents/BobLiterman/CSP2_Galaxy/CSP2/CSP2_env/env-d9b9114564458d9d-741b3de822f2aaca6c6caa4325c4afce/opt/mummer-3.23/scripts";
jpayne@69 19
jpayne@69 20 my $VERSION_INFO = q~
jpayne@69 21 DNAdiff version 1.3
jpayne@69 22 ~;
jpayne@69 23
jpayne@69 24 my $HELP_INFO = q~
jpayne@69 25 USAGE: dnadiff [options] <reference> <query>
jpayne@69 26 or dnadiff [options] -d <delta file>
jpayne@69 27
jpayne@69 28 DESCRIPTION:
jpayne@69 29 Run comparative analysis of two sequence sets using nucmer and its
jpayne@69 30 associated utilities with recommended parameters. See MUMmer
jpayne@69 31 documentation for a more detailed description of the
jpayne@69 32 output. Produces the following output files:
jpayne@69 33
jpayne@69 34 .report - Summary of alignments, differences and SNPs
jpayne@69 35 .delta - Standard nucmer alignment output
jpayne@69 36 .1delta - 1-to-1 alignment from delta-filter -1
jpayne@69 37 .mdelta - M-to-M alignment from delta-filter -m
jpayne@69 38 .1coords - 1-to-1 coordinates from show-coords -THrcl .1delta
jpayne@69 39 .mcoords - M-to-M coordinates from show-coords -THrcl .mdelta
jpayne@69 40 .snps - SNPs from show-snps -rlTHC .1delta
jpayne@69 41 .rdiff - Classified ref breakpoints from show-diff -rH .mdelta
jpayne@69 42 .qdiff - Classified qry breakpoints from show-diff -qH .mdelta
jpayne@69 43 .unref - Unaligned reference IDs and lengths (if applicable)
jpayne@69 44 .unqry - Unaligned query IDs and lengths (if applicable)
jpayne@69 45
jpayne@69 46 MANDATORY:
jpayne@69 47 reference Set the input reference multi-FASTA filename
jpayne@69 48 query Set the input query multi-FASTA filename
jpayne@69 49 or
jpayne@69 50 delta file Unfiltered .delta alignment file from nucmer
jpayne@69 51
jpayne@69 52 OPTIONS:
jpayne@69 53 -d|delta Provide precomputed delta file for analysis
jpayne@69 54 -h
jpayne@69 55 --help Display help information and exit
jpayne@69 56 -p|prefix Set the prefix of the output files (default "out")
jpayne@69 57 -V
jpayne@69 58 --version Display the version information and exit
jpayne@69 59 ~;
jpayne@69 60
jpayne@69 61
jpayne@69 62 my $USAGE_INFO = q~
jpayne@69 63 USAGE: dnadiff [options] <reference> <query>
jpayne@69 64 or dnadiff [options] -d <delta file>
jpayne@69 65 ~;
jpayne@69 66
jpayne@69 67
jpayne@69 68 my @DEPEND_INFO =
jpayne@69 69 (
jpayne@69 70 "$BIN_DIR/delta-filter",
jpayne@69 71 "$BIN_DIR/show-diff",
jpayne@69 72 "$BIN_DIR/show-snps",
jpayne@69 73 "$BIN_DIR/show-coords",
jpayne@69 74 "$BIN_DIR/nucmer",
jpayne@69 75 "$SCRIPT_DIR/Foundation.pm"
jpayne@69 76 );
jpayne@69 77
jpayne@69 78 my $DELTA_FILTER = "$BIN_DIR/delta-filter";
jpayne@69 79 my $SHOW_DIFF = "$BIN_DIR/show-diff";
jpayne@69 80 my $SHOW_SNPS = "$BIN_DIR/show-snps";
jpayne@69 81 my $SHOW_COORDS = "$BIN_DIR/show-coords";
jpayne@69 82 my $NUCMER = "$BIN_DIR/nucmer";
jpayne@69 83
jpayne@69 84 my $SNPBuff = 20; # required buffer around "good" snps
jpayne@69 85 my $OPT_Prefix = "out"; # prefix for all output files
jpayne@69 86 my $OPT_RefFile; # reference file
jpayne@69 87 my $OPT_QryFile; # query file
jpayne@69 88 my $OPT_DeltaFile; # unfiltered alignment file
jpayne@69 89 my $OPT_ReportFile = ".report"; # report file
jpayne@69 90 my $OPT_DeltaFile1 = ".1delta"; # 1-to-1 delta alignment
jpayne@69 91 my $OPT_DeltaFileM = ".mdelta"; # M-to-M delta alignment
jpayne@69 92 my $OPT_CoordsFile1 = ".1coords"; # 1-to-1 alignment coords
jpayne@69 93 my $OPT_CoordsFileM = ".mcoords"; # M-to-M alignment coords
jpayne@69 94 my $OPT_SnpsFile = ".snps"; # snps output file
jpayne@69 95 my $OPT_DiffRFile = ".rdiff"; # diffile for R
jpayne@69 96 my $OPT_DiffQFile = ".qdiff"; # diffile for Q
jpayne@69 97 my $OPT_UnRefFile = ".unref"; # unaligned ref IDs and lengths
jpayne@69 98 my $OPT_UnQryFile = ".unqry"; # unaligned qry IDs and lengths
jpayne@69 99
jpayne@69 100 my $TIGR; # TIGR Foundation object
jpayne@69 101
jpayne@69 102
jpayne@69 103 sub RunAlignment();
jpayne@69 104 sub RunFilter();
jpayne@69 105 sub RunCoords();
jpayne@69 106 sub RunSNPs();
jpayne@69 107 sub RunDiff();
jpayne@69 108 sub MakeReport();
jpayne@69 109
jpayne@69 110 sub FastaSizes($$);
jpayne@69 111
jpayne@69 112 sub FileOpen($$);
jpayne@69 113 sub FileClose($$);
jpayne@69 114
jpayne@69 115 sub GetOpt();
jpayne@69 116
jpayne@69 117
jpayne@69 118 #--------------------------------------------------------------------- main ----
jpayne@69 119 main:
jpayne@69 120 {
jpayne@69 121 GetOpt();
jpayne@69 122
jpayne@69 123 RunAlignment() unless defined($OPT_DeltaFile);
jpayne@69 124 RunFilter();
jpayne@69 125 RunCoords();
jpayne@69 126 RunSNPs();
jpayne@69 127 RunDiff();
jpayne@69 128 MakeReport();
jpayne@69 129
jpayne@69 130 exit(0);
jpayne@69 131 }
jpayne@69 132
jpayne@69 133
jpayne@69 134 #------------------------------------------------------------- RunAlignment ----
jpayne@69 135 # Run nucmer
jpayne@69 136 sub RunAlignment()
jpayne@69 137 {
jpayne@69 138 print STDERR "Building alignments\n";
jpayne@69 139 my $cmd = "$NUCMER --maxmatch -p $OPT_Prefix $OPT_RefFile $OPT_QryFile";
jpayne@69 140 my $err = "ERROR: Failed to run nucmer, aborting.\n";
jpayne@69 141
jpayne@69 142 system($cmd) == 0 or die $err;
jpayne@69 143 $OPT_DeltaFile = $OPT_Prefix . ".delta";
jpayne@69 144 }
jpayne@69 145
jpayne@69 146
jpayne@69 147 #---------------------------------------------------------------- RunFilter ----
jpayne@69 148 # Run delta-filter
jpayne@69 149 sub RunFilter()
jpayne@69 150 {
jpayne@69 151 print STDERR "Filtering alignments\n";
jpayne@69 152 my $cmd1 = "$DELTA_FILTER -1 $OPT_DeltaFile > $OPT_DeltaFile1";
jpayne@69 153 my $cmd2 = "$DELTA_FILTER -m $OPT_DeltaFile > $OPT_DeltaFileM";
jpayne@69 154 my $err = "ERROR: Failed to run delta-filter, aborting.\n";
jpayne@69 155
jpayne@69 156 system($cmd1) == 0 or die $err;
jpayne@69 157 system($cmd2) == 0 or die $err;
jpayne@69 158 }
jpayne@69 159
jpayne@69 160
jpayne@69 161 #------------------------------------------------------------------ RunSNPs ----
jpayne@69 162 # Run show-snps
jpayne@69 163 sub RunSNPs()
jpayne@69 164 {
jpayne@69 165 print STDERR "Analyzing SNPs\n";
jpayne@69 166 my $cmd = "$SHOW_SNPS -rlTHC $OPT_DeltaFile1 > $OPT_SnpsFile";
jpayne@69 167 my $err = "ERROR: Failed to run show-snps, aborting.\n";
jpayne@69 168
jpayne@69 169 system($cmd) == 0 or die $err;
jpayne@69 170 }
jpayne@69 171
jpayne@69 172
jpayne@69 173 #---------------------------------------------------------------- RunCoords ----
jpayne@69 174 # Run show-coords
jpayne@69 175 sub RunCoords()
jpayne@69 176 {
jpayne@69 177 print STDERR "Extracting alignment coordinates\n";
jpayne@69 178 my $cmd1 = "$SHOW_COORDS -rclTH $OPT_DeltaFile1 > $OPT_CoordsFile1";
jpayne@69 179 my $cmd2 = "$SHOW_COORDS -rclTH $OPT_DeltaFileM > $OPT_CoordsFileM";
jpayne@69 180 my $err = "ERROR: Failed to run show-coords, aborting.\n";
jpayne@69 181
jpayne@69 182 system($cmd1) == 0 or die $err;
jpayne@69 183 system($cmd2) == 0 or die $err;
jpayne@69 184 }
jpayne@69 185
jpayne@69 186
jpayne@69 187 #------------------------------------------------------------------ RunDiff ----
jpayne@69 188 # Run show-diff
jpayne@69 189 sub RunDiff()
jpayne@69 190 {
jpayne@69 191 print STDERR "Extracting alignment breakpoints\n";
jpayne@69 192 my $cmd1 = "$SHOW_DIFF -rH $OPT_DeltaFileM > $OPT_DiffRFile";
jpayne@69 193 my $cmd2 = "$SHOW_DIFF -qH $OPT_DeltaFileM > $OPT_DiffQFile";
jpayne@69 194 my $err = "ERROR: Failed to run show-diff, aborting.\n";
jpayne@69 195
jpayne@69 196 system($cmd1) == 0 or die $err;
jpayne@69 197 system($cmd2) == 0 or die $err;
jpayne@69 198 }
jpayne@69 199
jpayne@69 200
jpayne@69 201 #--------------------------------------------------------------- MakeReport ----
jpayne@69 202 # Output alignment report
jpayne@69 203 sub MakeReport()
jpayne@69 204 {
jpayne@69 205 print STDERR "Generating report file\n";
jpayne@69 206
jpayne@69 207 my ($fhi, $fho); # filehandle-in and filehandle-out
jpayne@69 208 my (%refs, %qrys) = ((),()); # R and Q ID->length
jpayne@69 209 my ($rqnAligns1, $rqnAlignsM) = (0,0); # alignment counter
jpayne@69 210 my ($rSumLen1, $qSumLen1) = (0,0); # alignment length sum
jpayne@69 211 my ($rSumLenM, $qSumLenM) = (0,0); # alignment length sum
jpayne@69 212 my ($rqSumLen1, $rqSumLenM) = (0,0); # combined alignment length sum
jpayne@69 213 my ($rqSumIdy1, $rqSumIdyM) = (0,0); # weighted alignment identity sum
jpayne@69 214 my ($qnIns, $rnIns) = (0,0); # insertion count
jpayne@69 215 my ($qSumIns, $rSumIns) = (0,0); # insertion length sum
jpayne@69 216 my ($qnTIns, $rnTIns) = (0,0); # tandem insertion count
jpayne@69 217 my ($qSumTIns, $rSumTIns) = (0,0); # tandem insertion length sum
jpayne@69 218 my ($qnInv, $rnInv) = (0,0); # inversion count
jpayne@69 219 my ($qnRel, $rnRel) = (0,0); # relocation count
jpayne@69 220 my ($qnTrn, $rnTrn) = (0,0); # translocation count
jpayne@69 221 my ($rnSeqs, $qnSeqs) = (0,0); # sequence count
jpayne@69 222 my ($rnASeqs, $qnASeqs) = (0,0); # aligned sequence count
jpayne@69 223 my ($rnBases, $qnBases) = (0,0); # bases count
jpayne@69 224 my ($rnABases, $qnABases) = (0,0); # aligned bases count
jpayne@69 225 my ($rnBrk, $qnBrk) = (0,0); # breakpoint count
jpayne@69 226 my ($rqnSNPs, $rqnIndels) = (0,0); # snp and indel counts
jpayne@69 227 my ($rqnGSNPs, $rqnGIndels) = (0,0); # good snp and indel counts
jpayne@69 228 my %rqSNPs = # SNP hash
jpayne@69 229 ( "."=>{"A"=>0,"C"=>0,"G"=>0,"T"=>0},
jpayne@69 230 "A"=>{"."=>0,"C"=>0,"G"=>0,"T"=>0},
jpayne@69 231 "C"=>{"."=>0,"A"=>0,"G"=>0,"T"=>0},
jpayne@69 232 "G"=>{"."=>0,"A"=>0,"C"=>0,"T"=>0},
jpayne@69 233 "T"=>{"."=>0,"A"=>0,"C"=>0,"G"=>0} );
jpayne@69 234 my %rqGSNPs = # good SNP hash
jpayne@69 235 ( "."=>{"A"=>0,"C"=>0,"G"=>0,"T"=>0},
jpayne@69 236 "A"=>{"."=>0,"C"=>0,"G"=>0,"T"=>0},
jpayne@69 237 "C"=>{"."=>0,"A"=>0,"G"=>0,"T"=>0},
jpayne@69 238 "G"=>{"."=>0,"A"=>0,"C"=>0,"T"=>0},
jpayne@69 239 "T"=>{"."=>0,"A"=>0,"C"=>0,"G"=>0} );
jpayne@69 240
jpayne@69 241 my $header; # delta header
jpayne@69 242
jpayne@69 243 #-- Get delta header
jpayne@69 244 $fhi = FileOpen("<", $OPT_DeltaFile);
jpayne@69 245 $header .= <$fhi>;
jpayne@69 246 $header .= <$fhi>;
jpayne@69 247 $header .= "\n";
jpayne@69 248 FileClose($fhi, $OPT_DeltaFile);
jpayne@69 249
jpayne@69 250 #-- Collect all reference and query IDs and lengths
jpayne@69 251 FastaSizes($OPT_RefFile, \%refs);
jpayne@69 252 FastaSizes($OPT_QryFile, \%qrys);
jpayne@69 253
jpayne@69 254 #-- Count ref and qry seqs and lengths
jpayne@69 255 foreach ( values(%refs) ) {
jpayne@69 256 $rnSeqs++;
jpayne@69 257 $rnBases += $_;
jpayne@69 258 }
jpayne@69 259 foreach ( values(%qrys) ) {
jpayne@69 260 $qnSeqs++;
jpayne@69 261 $qnBases += $_;
jpayne@69 262 }
jpayne@69 263
jpayne@69 264 #-- Count aligned seqs, aligned bases, and breakpoints for each R and Q
jpayne@69 265 $fhi = FileOpen("<", $OPT_CoordsFileM);
jpayne@69 266 while (<$fhi>) {
jpayne@69 267 chomp;
jpayne@69 268 my @A = split "\t";
jpayne@69 269 scalar(@A) == 13
jpayne@69 270 or die "ERROR: Unrecognized format $OPT_CoordsFileM, aborting.\n";
jpayne@69 271
jpayne@69 272 #-- Add to M-to-M alignment counts
jpayne@69 273 $rqnAlignsM++;
jpayne@69 274 $rSumLenM += $A[4];
jpayne@69 275 $qSumLenM += $A[5];
jpayne@69 276 $rqSumIdyM += ($A[6] / 100.0) * ($A[4] + $A[5]);
jpayne@69 277 $rqSumLenM += ($A[4] + $A[5]);
jpayne@69 278
jpayne@69 279 #-- If new ID, add to sequence and base count
jpayne@69 280 if ( $refs{$A[11]} > 0 ) {
jpayne@69 281 $rnASeqs++;
jpayne@69 282 $rnABases += $refs{$A[11]};
jpayne@69 283 $refs{$A[11]} *= -1; # If ref has alignment, length will be -neg
jpayne@69 284 }
jpayne@69 285 if ( $qrys{$A[12]} > 0 ) {
jpayne@69 286 $qnASeqs++;
jpayne@69 287 $qnABases += $qrys{$A[12]};
jpayne@69 288 $qrys{$A[12]} *= -1; # If qry has alignment, length will be -neg
jpayne@69 289 }
jpayne@69 290
jpayne@69 291 #-- Add to breakpoint counts
jpayne@69 292 my ($lo, $hi);
jpayne@69 293 if ( $A[0] < $A[1] ) { $lo = $A[0]; $hi = $A[1]; }
jpayne@69 294 else { $lo = $A[1]; $hi = $A[0]; }
jpayne@69 295 $rnBrk++ if ( $lo != 1 );
jpayne@69 296 $rnBrk++ if ( $hi != $A[7] );
jpayne@69 297
jpayne@69 298 if ( $A[2] < $A[3] ) { $lo = $A[2]; $hi = $A[3]; }
jpayne@69 299 else { $lo = $A[3]; $hi = $A[2]; }
jpayne@69 300 $qnBrk++ if ( $lo != 1 );
jpayne@69 301 $qnBrk++ if ( $hi != $A[8] );
jpayne@69 302 }
jpayne@69 303 FileClose($fhi, $OPT_CoordsFileM);
jpayne@69 304
jpayne@69 305 #-- Calculate average %idy, length, etc.
jpayne@69 306 $fhi = FileOpen("<", $OPT_CoordsFile1);
jpayne@69 307 while (<$fhi>) {
jpayne@69 308 chomp;
jpayne@69 309 my @A = split "\t";
jpayne@69 310 scalar(@A) == 13
jpayne@69 311 or die "ERROR: Unrecognized format $OPT_CoordsFile1, aborting.\n";
jpayne@69 312
jpayne@69 313 #-- Add to 1-to-1 alignment counts
jpayne@69 314 $rqnAligns1++;
jpayne@69 315 $rSumLen1 += $A[4];
jpayne@69 316 $qSumLen1 += $A[5];
jpayne@69 317 $rqSumIdy1 += ($A[6] / 100.0) * ($A[4] + $A[5]);
jpayne@69 318 $rqSumLen1 += ($A[4] + $A[5]);
jpayne@69 319 }
jpayne@69 320 FileClose($fhi, $OPT_CoordsFile1);
jpayne@69 321
jpayne@69 322 #-- If you are reading this, you need to get out more...
jpayne@69 323
jpayne@69 324 #-- Count reference diff features and indels
jpayne@69 325 $fhi = FileOpen("<", $OPT_DiffRFile);
jpayne@69 326 while (<$fhi>) {
jpayne@69 327 chomp;
jpayne@69 328 my @A = split "\t";
jpayne@69 329 defined($A[4])
jpayne@69 330 or die "ERROR: Unrecognized format $OPT_DiffRFile, aborting.\n";
jpayne@69 331 my $gap = $A[4];
jpayne@69 332 my $ins = $gap;
jpayne@69 333
jpayne@69 334 #-- Add to tandem insertion counts
jpayne@69 335 if ( $A[1] eq "GAP" ) {
jpayne@69 336 scalar(@A) == 7
jpayne@69 337 or die "ERROR: Unrecognized format $OPT_DiffRFile, aborting.\n";
jpayne@69 338 $ins = $A[6] if ( $A[6] > $gap );
jpayne@69 339 if ( $A[4] <= 0 && $A[5] <= 0 && $A[6] > 0 ) {
jpayne@69 340 $rnTIns++;
jpayne@69 341 $rSumTIns += $A[6];
jpayne@69 342 }
jpayne@69 343 }
jpayne@69 344
jpayne@69 345 #-- Remove unaligned sequence from count
jpayne@69 346 if ( $A[1] ne "DUP" ) {
jpayne@69 347 $rnABases -= $gap if ( $gap > 0 );
jpayne@69 348 }
jpayne@69 349
jpayne@69 350 #-- Add to insertion count
jpayne@69 351 if ( $ins > 0 ) {
jpayne@69 352 $rnIns++;
jpayne@69 353 $rSumIns += $ins;
jpayne@69 354 }
jpayne@69 355
jpayne@69 356 #-- Add to rearrangement counts
jpayne@69 357 $rnInv++ if ( $A[1] eq "INV" );
jpayne@69 358 $rnRel++ if ( $A[1] eq "JMP" );
jpayne@69 359 $rnTrn++ if ( $A[1] eq "SEQ" );
jpayne@69 360 }
jpayne@69 361 FileClose($fhi, $OPT_DiffRFile);
jpayne@69 362
jpayne@69 363 #-- Count query diff features and indels
jpayne@69 364 $fhi = FileOpen("<", $OPT_DiffQFile);
jpayne@69 365 while (<$fhi>) {
jpayne@69 366 chomp;
jpayne@69 367 my @A = split "\t";
jpayne@69 368 defined($A[4])
jpayne@69 369 or die "ERROR: Unrecognized format $OPT_DiffRFile, aborting.\n";
jpayne@69 370 my $gap = $A[4];
jpayne@69 371 my $ins = $gap;
jpayne@69 372
jpayne@69 373 #-- Add to tandem insertion counts
jpayne@69 374 if ( $A[1] eq "GAP" ) {
jpayne@69 375 scalar(@A) == 7
jpayne@69 376 or die "ERROR: Unrecognized format $OPT_DiffRFile, aborting.\n";
jpayne@69 377 $ins = $A[6] if ( $A[6] > $gap );
jpayne@69 378 if ( $A[4] <= 0 && $A[5] <= 0 && $A[6] > 0 ) {
jpayne@69 379 $qnTIns++;
jpayne@69 380 $qSumTIns += $A[6];
jpayne@69 381 }
jpayne@69 382 }
jpayne@69 383
jpayne@69 384 #-- Remove unaligned sequence from count
jpayne@69 385 if ( $A[1] ne "DUP" ) {
jpayne@69 386 $qnABases -= $gap if ( $gap > 0 );
jpayne@69 387 }
jpayne@69 388
jpayne@69 389 #-- Add to insertion count
jpayne@69 390 if ( $ins > 0 ) {
jpayne@69 391 $qnIns++;
jpayne@69 392 $qSumIns += $ins;
jpayne@69 393 }
jpayne@69 394
jpayne@69 395 #-- Add to rearrangement counts
jpayne@69 396 $qnInv++ if ( $A[1] eq "INV" );
jpayne@69 397 $qnRel++ if ( $A[1] eq "JMP" );
jpayne@69 398 $qnTrn++ if ( $A[1] eq "SEQ" );
jpayne@69 399 }
jpayne@69 400 FileClose($fhi, $OPT_DiffQFile);
jpayne@69 401
jpayne@69 402 #-- Count SNPs
jpayne@69 403 $fhi = FileOpen("<", $OPT_SnpsFile);
jpayne@69 404 while(<$fhi>) {
jpayne@69 405 chomp;
jpayne@69 406 my @A = split "\t";
jpayne@69 407 scalar(@A) == 12
jpayne@69 408 or die "ERROR: Unrecognized format $OPT_SnpsFile, aborting\n";
jpayne@69 409
jpayne@69 410 my $r = uc($A[1]);
jpayne@69 411 my $q = uc($A[2]);
jpayne@69 412
jpayne@69 413 #-- Plain SNPs
jpayne@69 414 $rqSNPs{$r}{$q}++;
jpayne@69 415 if ( !exists($rqSNPs{$q}{$r}) ) { $rqSNPs{$q}{$r} = 0; }
jpayne@69 416 if ( $r eq '.' || $q eq '.' ) { $rqnIndels++; }
jpayne@69 417 else { $rqnSNPs++; }
jpayne@69 418
jpayne@69 419 #-- Good SNPs with sufficient match buffer
jpayne@69 420 if ( $A[4] >= $SNPBuff ) {
jpayne@69 421 $rqGSNPs{$r}{$q}++;
jpayne@69 422 if ( !exists($rqGSNPs{$q}{$r}) ) { $rqGSNPs{$q}{$r} = 0; }
jpayne@69 423 if ( $r eq '.' || $q eq '.' ) { $rqnGIndels++; }
jpayne@69 424 else { $rqnGSNPs++; }
jpayne@69 425 }
jpayne@69 426 }
jpayne@69 427 FileClose($fhi, $OPT_SnpsFile);
jpayne@69 428
jpayne@69 429
jpayne@69 430 #-- Output report
jpayne@69 431 $fho = FileOpen(">", $OPT_ReportFile);
jpayne@69 432
jpayne@69 433 print $fho $header;
jpayne@69 434 printf $fho "%-15s %20s %20s\n", "", "[REF]", "[QRY]";
jpayne@69 435
jpayne@69 436 print $fho "[Sequences]\n";
jpayne@69 437
jpayne@69 438 printf $fho "%-15s %20d %20d\n",
jpayne@69 439 "TotalSeqs", $rnSeqs, $qnSeqs;
jpayne@69 440 printf $fho "%-15s %20s %20s\n",
jpayne@69 441 "AlignedSeqs",
jpayne@69 442 ( sprintf "%10d(%.2f%%)",
jpayne@69 443 $rnASeqs, ($rnSeqs ? $rnASeqs / $rnSeqs * 100.0 : 0) ),
jpayne@69 444 ( sprintf "%10d(%.2f%%)",
jpayne@69 445 $qnASeqs, ($rnSeqs ? $qnASeqs / $qnSeqs * 100.0 : 0) );
jpayne@69 446 printf $fho "%-15s %20s %20s\n",
jpayne@69 447 "UnalignedSeqs",
jpayne@69 448 ( sprintf "%10d(%.2f%%)",
jpayne@69 449 $rnSeqs - $rnASeqs,
jpayne@69 450 ($rnSeqs ? ($rnSeqs - $rnASeqs) / $rnSeqs * 100.0 : 0) ),
jpayne@69 451 ( sprintf "%10d(%.2f%%)",
jpayne@69 452 $qnSeqs - $qnASeqs,
jpayne@69 453 ($qnSeqs ? ($qnSeqs - $qnASeqs) / $qnSeqs * 100.0 : 0) );
jpayne@69 454
jpayne@69 455 print $fho "\n[Bases]\n";
jpayne@69 456
jpayne@69 457 printf $fho "%-15s %20d %20d\n",
jpayne@69 458 "TotalBases", $rnBases, $qnBases;
jpayne@69 459 printf $fho "%-15s %20s %20s\n",
jpayne@69 460 "AlignedBases",
jpayne@69 461 ( sprintf "%10d(%.2f%%)",
jpayne@69 462 $rnABases, ($rnBases ? $rnABases / $rnBases * 100.0 : 0) ),
jpayne@69 463 ( sprintf "%10d(%.2f%%)",
jpayne@69 464 $qnABases, ($qnBases ? $qnABases / $qnBases * 100.0 : 0) );
jpayne@69 465 printf $fho "%-15s %20s %20s\n",
jpayne@69 466 "UnalignedBases",
jpayne@69 467 ( sprintf "%10d(%.2f%%)",
jpayne@69 468 $rnBases - $rnABases,
jpayne@69 469 ($rnBases ? ($rnBases - $rnABases) / $rnBases * 100.0 : 0) ),
jpayne@69 470 ( sprintf "%10d(%.2f%%)",
jpayne@69 471 $qnBases - $qnABases,
jpayne@69 472 ($qnBases ? ($qnBases - $qnABases) / $qnBases * 100.0 : 0) );
jpayne@69 473
jpayne@69 474 print $fho "\n[Alignments]\n";
jpayne@69 475
jpayne@69 476 printf $fho "%-15s %20d %20d\n",
jpayne@69 477 "1-to-1", $rqnAligns1, $rqnAligns1;
jpayne@69 478 printf $fho "%-15s %20d %20d\n",
jpayne@69 479 "TotalLength", $rSumLen1, $qSumLen1;
jpayne@69 480 printf $fho "%-15s %20.2f %20.2f\n",
jpayne@69 481 "AvgLength",
jpayne@69 482 ($rqnAligns1 ? $rSumLen1 / $rqnAligns1 : 0),
jpayne@69 483 ($rqnAligns1 ? $qSumLen1 / $rqnAligns1 : 0);
jpayne@69 484 printf $fho "%-15s %20.2f %20.2f\n",
jpayne@69 485 "AvgIdentity",
jpayne@69 486 ($rqSumLen1 ? $rqSumIdy1 / $rqSumLen1 * 100.0 : 0),
jpayne@69 487 ($rqSumLen1 ? $rqSumIdy1 / $rqSumLen1 * 100.0 : 0);
jpayne@69 488
jpayne@69 489 print $fho "\n";
jpayne@69 490
jpayne@69 491 printf $fho "%-15s %20d %20d\n",
jpayne@69 492 "M-to-M", $rqnAlignsM, $rqnAlignsM;
jpayne@69 493 printf $fho "%-15s %20d %20d\n",
jpayne@69 494 "TotalLength", $rSumLenM, $qSumLenM;
jpayne@69 495 printf $fho "%-15s %20.2f %20.2f\n",
jpayne@69 496 "AvgLength",
jpayne@69 497 ($rqnAlignsM ? $rSumLenM / $rqnAlignsM : 0),
jpayne@69 498 ($rqnAlignsM ? $qSumLenM / $rqnAlignsM : 0);
jpayne@69 499 printf $fho "%-15s %20.2f %20.2f\n",
jpayne@69 500 "AvgIdentity",
jpayne@69 501 ($rqSumLenM ? $rqSumIdyM / $rqSumLenM * 100.0 : 0),
jpayne@69 502 ($rqSumLenM ? $rqSumIdyM / $rqSumLenM * 100.0 : 0);
jpayne@69 503
jpayne@69 504 print $fho "\n[Feature Estimates]\n";
jpayne@69 505
jpayne@69 506 printf $fho "%-15s %20d %20d\n",
jpayne@69 507 "Breakpoints", $rnBrk, $qnBrk;
jpayne@69 508 printf $fho "%-15s %20d %20d\n",
jpayne@69 509 "Relocations", $rnRel, $qnRel;
jpayne@69 510 printf $fho "%-15s %20d %20d\n",
jpayne@69 511 "Translocations", $rnTrn, $qnTrn;
jpayne@69 512 printf $fho "%-15s %20d %20d\n",
jpayne@69 513 "Inversions", $rnInv, $qnInv;
jpayne@69 514
jpayne@69 515 print $fho "\n";
jpayne@69 516
jpayne@69 517 printf $fho "%-15s %20d %20d\n",
jpayne@69 518 "Insertions", $rnIns, $qnIns;
jpayne@69 519 printf $fho "%-15s %20d %20d\n",
jpayne@69 520 "InsertionSum", $rSumIns, $qSumIns;
jpayne@69 521 printf $fho "%-15s %20.2f %20.2f\n",
jpayne@69 522 "InsertionAvg",
jpayne@69 523 ($rnIns ? $rSumIns / $rnIns : 0),
jpayne@69 524 ($qnIns ? $qSumIns / $qnIns : 0);
jpayne@69 525
jpayne@69 526 print $fho "\n";
jpayne@69 527
jpayne@69 528 printf $fho "%-15s %20d %20d\n",
jpayne@69 529 "TandemIns", $rnTIns, $qnTIns;
jpayne@69 530 printf $fho "%-15s %20d %20d\n",
jpayne@69 531 "TandemInsSum", $rSumTIns, $qSumTIns;
jpayne@69 532 printf $fho "%-15s %20.2f %20.2f\n",
jpayne@69 533 "TandemInsAvg",
jpayne@69 534 ($rnTIns ? $rSumTIns / $rnTIns : 0),
jpayne@69 535 ($qnTIns ? $qSumTIns / $qnTIns : 0);
jpayne@69 536
jpayne@69 537 print $fho "\n[SNPs]\n";
jpayne@69 538
jpayne@69 539 printf $fho "%-15s %20d %20d\n",
jpayne@69 540 "TotalSNPs", $rqnSNPs, $rqnSNPs;
jpayne@69 541 foreach my $r (keys %rqSNPs) {
jpayne@69 542 foreach my $q (keys %{$rqSNPs{$r}}) {
jpayne@69 543 if ( $r ne "." && $q ne "." ) {
jpayne@69 544 printf $fho "%-15s %20s %20s\n",
jpayne@69 545 "$r$q",
jpayne@69 546 ( sprintf "%10d(%.2f%%)",
jpayne@69 547 $rqSNPs{$r}{$q},
jpayne@69 548 ($rqnSNPs ? $rqSNPs{$r}{$q} / $rqnSNPs * 100.0 : 0) ),
jpayne@69 549 ( sprintf "%10d(%.2f%%)",
jpayne@69 550 $rqSNPs{$q}{$r},
jpayne@69 551 ($rqnSNPs ? $rqSNPs{$q}{$r} / $rqnSNPs * 100.0 : 0) );
jpayne@69 552 }
jpayne@69 553 }
jpayne@69 554 }
jpayne@69 555
jpayne@69 556 print $fho "\n";
jpayne@69 557
jpayne@69 558 printf $fho "%-15s %20d %20d\n",
jpayne@69 559 "TotalGSNPs", $rqnGSNPs, $rqnGSNPs;
jpayne@69 560 foreach my $r (keys %rqGSNPs) {
jpayne@69 561 foreach my $q (keys %{$rqGSNPs{$r}}) {
jpayne@69 562 if ( $r ne "." && $q ne "." ) {
jpayne@69 563 printf $fho "%-15s %20s %20s\n",
jpayne@69 564 "$r$q",
jpayne@69 565 ( sprintf "%10d(%.2f%%)",
jpayne@69 566 $rqGSNPs{$r}{$q},
jpayne@69 567 ($rqnGSNPs ? $rqGSNPs{$r}{$q} / $rqnGSNPs * 100.0 : 0) ),
jpayne@69 568 ( sprintf "%10d(%.2f%%)",
jpayne@69 569 $rqGSNPs{$q}{$r},
jpayne@69 570 ($rqnGSNPs ? $rqGSNPs{$q}{$r} / $rqnGSNPs * 100.0 : 0) );
jpayne@69 571 }
jpayne@69 572 }
jpayne@69 573 }
jpayne@69 574
jpayne@69 575 print $fho "\n";
jpayne@69 576
jpayne@69 577 printf $fho "%-15s %20d %20d\n",
jpayne@69 578 "TotalIndels", $rqnIndels, $rqnIndels;
jpayne@69 579 foreach my $r (keys %rqSNPs) {
jpayne@69 580 foreach my $q (keys %{$rqSNPs{$r}}) {
jpayne@69 581 if ( $q eq "." ) {
jpayne@69 582 printf $fho "%-15s %20s %20s\n",
jpayne@69 583 "$r$q",
jpayne@69 584 ( sprintf "%10d(%.2f%%)",
jpayne@69 585 $rqSNPs{$r}{$q},
jpayne@69 586 ($rqnIndels ? $rqSNPs{$r}{$q} / $rqnIndels * 100.0 : 0) ),
jpayne@69 587 ( sprintf "%10d(%.2f%%)",
jpayne@69 588 $rqSNPs{$q}{$r},
jpayne@69 589 ($rqnIndels ? $rqSNPs{$q}{$r} / $rqnIndels * 100.0 : 0) );
jpayne@69 590 }
jpayne@69 591 }
jpayne@69 592 }
jpayne@69 593 foreach my $r (keys %rqSNPs) {
jpayne@69 594 foreach my $q (keys %{$rqSNPs{$r}}) {
jpayne@69 595 if ( $r eq "." ) {
jpayne@69 596 printf $fho "%-15s %20s %20s\n",
jpayne@69 597 "$r$q",
jpayne@69 598 ( sprintf "%10d(%.2f%%)",
jpayne@69 599 $rqSNPs{$r}{$q},
jpayne@69 600 ($rqnIndels ? $rqSNPs{$r}{$q} / $rqnIndels * 100.0 : 0) ),
jpayne@69 601 ( sprintf "%10d(%.2f%%)",
jpayne@69 602 $rqSNPs{$q}{$r},
jpayne@69 603 ($rqnIndels ? $rqSNPs{$q}{$r} / $rqnIndels * 100.0 : 0) );
jpayne@69 604 }
jpayne@69 605 }
jpayne@69 606 }
jpayne@69 607
jpayne@69 608 print $fho "\n";
jpayne@69 609
jpayne@69 610 printf $fho "%-15s %20d %20d\n",
jpayne@69 611 "TotalGIndels", $rqnGIndels, $rqnGIndels;
jpayne@69 612 foreach my $r (keys %rqGSNPs) {
jpayne@69 613 foreach my $q (keys %{$rqGSNPs{$r}}) {
jpayne@69 614 if ( $q eq "." ) {
jpayne@69 615 printf $fho "%-15s %20s %20s\n",
jpayne@69 616 "$r$q",
jpayne@69 617 ( sprintf "%10d(%.2f%%)",
jpayne@69 618 $rqGSNPs{$r}{$q},
jpayne@69 619 ($rqnGIndels ? $rqGSNPs{$r}{$q} / $rqnGIndels * 100.0 : 0) ),
jpayne@69 620 ( sprintf "%10d(%.2f%%)",
jpayne@69 621 $rqGSNPs{$q}{$r},
jpayne@69 622 ($rqnGIndels ? $rqGSNPs{$q}{$r} / $rqnGIndels * 100.0 : 0) );
jpayne@69 623 }
jpayne@69 624 }
jpayne@69 625 }
jpayne@69 626 foreach my $r (keys %rqGSNPs) {
jpayne@69 627 foreach my $q (keys %{$rqGSNPs{$r}}) {
jpayne@69 628 if ( $r eq "." ) {
jpayne@69 629 printf $fho "%-15s %20s %20s\n",
jpayne@69 630 "$r$q",
jpayne@69 631 ( sprintf "%10d(%.2f%%)",
jpayne@69 632 $rqGSNPs{$r}{$q},
jpayne@69 633 ($rqnGIndels ? $rqGSNPs{$r}{$q} / $rqnGIndels * 100.0 : 0) ),
jpayne@69 634 ( sprintf "%10d(%.2f%%)",
jpayne@69 635 $rqGSNPs{$q}{$r},
jpayne@69 636 ($rqnGIndels ? $rqGSNPs{$q}{$r} / $rqnGIndels * 100.0 : 0) );
jpayne@69 637 }
jpayne@69 638 }
jpayne@69 639 }
jpayne@69 640
jpayne@69 641 FileClose($fho, $OPT_ReportFile);
jpayne@69 642
jpayne@69 643
jpayne@69 644 #-- Output unaligned reference and query IDs, if applicable
jpayne@69 645 if ( $rnSeqs != $rnASeqs ) {
jpayne@69 646 $fho = FileOpen(">", $OPT_UnRefFile);
jpayne@69 647 while ( my ($key, $val) = each(%refs) ) {
jpayne@69 648 print $fho "$key\tUNI\t1\t$val\t$val\n" unless $val < 0;
jpayne@69 649 }
jpayne@69 650 FileClose($fho, $OPT_UnRefFile);
jpayne@69 651 }
jpayne@69 652 if ( $qnSeqs != $qnASeqs ) {
jpayne@69 653 $fho = FileOpen(">", $OPT_UnQryFile);
jpayne@69 654 while ( my ($key, $val) = each(%qrys) ) {
jpayne@69 655 print $fho "$key\tUNI\t1\t$val\t$val\n" unless $val < 0;
jpayne@69 656 }
jpayne@69 657 FileClose($fho, $OPT_UnQryFile);
jpayne@69 658 }
jpayne@69 659 }
jpayne@69 660
jpayne@69 661
jpayne@69 662 #--------------------------------------------------------------- FastaSizes ----
jpayne@69 663 # Compute lengths for a multi-fasta file and store in hash reference
jpayne@69 664 sub FastaSizes($$)
jpayne@69 665 {
jpayne@69 666
jpayne@69 667 my $file = shift;
jpayne@69 668 my $href = shift;
jpayne@69 669 my ($tag, $len);
jpayne@69 670
jpayne@69 671 my $fhi = FileOpen("<", $file);
jpayne@69 672 while (<$fhi>) {
jpayne@69 673 chomp;
jpayne@69 674
jpayne@69 675 if ( /^>/ ) {
jpayne@69 676 $href->{$tag} = $len if defined($tag);
jpayne@69 677 ($tag) = /^>(\S+)/;
jpayne@69 678 $len = 0;
jpayne@69 679 } else {
jpayne@69 680 if ( /\s/ ) {
jpayne@69 681 die "ERROR: Whitespace found in FastA $file, aborting.\n";
jpayne@69 682 }
jpayne@69 683 $len += length;
jpayne@69 684 }
jpayne@69 685 }
jpayne@69 686 $href->{$tag} = $len if defined($tag);
jpayne@69 687 FileClose($fhi, $file);
jpayne@69 688 }
jpayne@69 689
jpayne@69 690
jpayne@69 691 #----------------------------------------------------------------- FileOpen ----
jpayne@69 692 # Open file, return filehandle, or die
jpayne@69 693 sub FileOpen($$)
jpayne@69 694 {
jpayne@69 695 my ($mode, $name) = @_;
jpayne@69 696 my $fhi;
jpayne@69 697 open($fhi, $mode, $name)
jpayne@69 698 or die "ERROR: Could not open $name, aborting. $!\n";
jpayne@69 699 return $fhi;
jpayne@69 700 }
jpayne@69 701
jpayne@69 702
jpayne@69 703 #---------------------------------------------------------------- FileClose ----
jpayne@69 704 # Close file, or die
jpayne@69 705 sub FileClose($$)
jpayne@69 706 {
jpayne@69 707 my ($fho, $name) = @_;
jpayne@69 708 close($fho) or die "ERROR: Could not close $name, aborting. $!\n"
jpayne@69 709 }
jpayne@69 710
jpayne@69 711
jpayne@69 712 #------------------------------------------------------------------- GetOpt ----
jpayne@69 713 # Get command options and check file permissions
jpayne@69 714 sub GetOpt()
jpayne@69 715 {
jpayne@69 716 #-- Initialize TIGR::Foundation
jpayne@69 717 $TIGR = new TIGR::Foundation;
jpayne@69 718 if ( !defined($TIGR) ) {
jpayne@69 719 print STDERR "ERROR: TIGR::Foundation could not be initialized";
jpayne@69 720 exit(1);
jpayne@69 721 }
jpayne@69 722
jpayne@69 723 #-- Set help and usage information
jpayne@69 724 $TIGR->setHelpInfo($HELP_INFO);
jpayne@69 725 $TIGR->setUsageInfo($USAGE_INFO);
jpayne@69 726 $TIGR->setVersionInfo($VERSION_INFO);
jpayne@69 727 $TIGR->addDependInfo(@DEPEND_INFO);
jpayne@69 728
jpayne@69 729 #-- Get options
jpayne@69 730 my $err = !$TIGR->TIGR_GetOptions
jpayne@69 731 (
jpayne@69 732 "d|delta=s" => \$OPT_DeltaFile,
jpayne@69 733 "p|prefix=s" => \$OPT_Prefix,
jpayne@69 734 );
jpayne@69 735
jpayne@69 736 #-- Check if the parsing was successful
jpayne@69 737 if ( $err
jpayne@69 738 || (defined($OPT_DeltaFile) && scalar(@ARGV) != 0)
jpayne@69 739 || (!defined($OPT_DeltaFile) && scalar(@ARGV) != 2) ) {
jpayne@69 740 $TIGR->printUsageInfo();
jpayne@69 741 print STDERR "Try '$0 -h' for more information.\n";
jpayne@69 742 exit(1);
jpayne@69 743 }
jpayne@69 744
jpayne@69 745 my @errs;
jpayne@69 746
jpayne@69 747 $TIGR->isExecutableFile($DELTA_FILTER)
jpayne@69 748 or push(@errs, $DELTA_FILTER);
jpayne@69 749
jpayne@69 750 $TIGR->isExecutableFile($SHOW_DIFF)
jpayne@69 751 or push(@errs, $SHOW_DIFF);
jpayne@69 752
jpayne@69 753 $TIGR->isExecutableFile($SHOW_SNPS)
jpayne@69 754 or push(@errs, $SHOW_SNPS);
jpayne@69 755
jpayne@69 756 $TIGR->isExecutableFile($SHOW_COORDS)
jpayne@69 757 or push(@errs, $SHOW_COORDS);
jpayne@69 758
jpayne@69 759 $TIGR->isExecutableFile($NUCMER)
jpayne@69 760 or push(@errs, $NUCMER);
jpayne@69 761
jpayne@69 762 if ( defined($OPT_DeltaFile) ) {
jpayne@69 763 $TIGR->isReadableFile($OPT_DeltaFile)
jpayne@69 764 or push(@errs, $OPT_DeltaFile);
jpayne@69 765
jpayne@69 766 my $fhi = FileOpen("<", $OPT_DeltaFile);
jpayne@69 767 $_ = <$fhi>;
jpayne@69 768 FileClose($fhi, $OPT_DeltaFile);
jpayne@69 769
jpayne@69 770 ($OPT_RefFile, $OPT_QryFile) = /^(.+) (.+)$/;
jpayne@69 771 }
jpayne@69 772 else {
jpayne@69 773 $OPT_RefFile = File::Spec->rel2abs($ARGV[0]);
jpayne@69 774 $OPT_QryFile = File::Spec->rel2abs($ARGV[1]);
jpayne@69 775 }
jpayne@69 776
jpayne@69 777 $TIGR->isReadableFile($OPT_RefFile)
jpayne@69 778 or push(@errs, $OPT_RefFile);
jpayne@69 779
jpayne@69 780 $TIGR->isReadableFile($OPT_QryFile)
jpayne@69 781 or push(@errs, $OPT_QryFile);
jpayne@69 782
jpayne@69 783 $OPT_ReportFile = $OPT_Prefix . $OPT_ReportFile;
jpayne@69 784 $TIGR->isCreatableFile("$OPT_ReportFile")
jpayne@69 785 or $TIGR->isWritableFile("$OPT_ReportFile")
jpayne@69 786 or push(@errs, "$OPT_ReportFile");
jpayne@69 787
jpayne@69 788 $OPT_DeltaFile1 = $OPT_Prefix . $OPT_DeltaFile1;
jpayne@69 789 $TIGR->isCreatableFile("$OPT_DeltaFile1")
jpayne@69 790 or $TIGR->isWritableFile("$OPT_DeltaFile1")
jpayne@69 791 or push(@errs, "$OPT_DeltaFile1");
jpayne@69 792
jpayne@69 793 $OPT_DeltaFileM = $OPT_Prefix . $OPT_DeltaFileM;
jpayne@69 794 $TIGR->isCreatableFile("$OPT_DeltaFileM")
jpayne@69 795 or $TIGR->isWritableFile("$OPT_DeltaFileM")
jpayne@69 796 or push(@errs, "$OPT_DeltaFileM");
jpayne@69 797
jpayne@69 798 $OPT_CoordsFile1 = $OPT_Prefix . $OPT_CoordsFile1;
jpayne@69 799 $TIGR->isCreatableFile("$OPT_CoordsFile1")
jpayne@69 800 or $TIGR->isWritableFile("$OPT_CoordsFile1")
jpayne@69 801 or push(@errs, "$OPT_CoordsFile1");
jpayne@69 802
jpayne@69 803 $OPT_CoordsFileM = $OPT_Prefix . $OPT_CoordsFileM;
jpayne@69 804 $TIGR->isCreatableFile("$OPT_CoordsFileM")
jpayne@69 805 or $TIGR->isWritableFile("$OPT_CoordsFileM")
jpayne@69 806 or push(@errs, "$OPT_CoordsFileM");
jpayne@69 807
jpayne@69 808 $OPT_SnpsFile = $OPT_Prefix . $OPT_SnpsFile;
jpayne@69 809 $TIGR->isCreatableFile("$OPT_SnpsFile")
jpayne@69 810 or $TIGR->isWritableFile("$OPT_SnpsFile")
jpayne@69 811 or push(@errs, "$OPT_SnpsFile");
jpayne@69 812
jpayne@69 813 $OPT_DiffRFile = $OPT_Prefix . $OPT_DiffRFile;
jpayne@69 814 $TIGR->isCreatableFile("$OPT_DiffRFile")
jpayne@69 815 or $TIGR->isWritableFile("$OPT_DiffRFile")
jpayne@69 816 or push(@errs, "$OPT_DiffRFile");
jpayne@69 817
jpayne@69 818 $OPT_DiffQFile = $OPT_Prefix . $OPT_DiffQFile;
jpayne@69 819 $TIGR->isCreatableFile("$OPT_DiffQFile")
jpayne@69 820 or $TIGR->isWritableFile("$OPT_DiffQFile")
jpayne@69 821 or push(@errs, "$OPT_DiffQFile");
jpayne@69 822
jpayne@69 823 $OPT_UnRefFile = $OPT_Prefix . $OPT_UnRefFile;
jpayne@69 824 $TIGR->isCreatableFile("$OPT_UnRefFile")
jpayne@69 825 or $TIGR->isWritableFile("$OPT_UnRefFile")
jpayne@69 826 or push(@errs, "$OPT_UnRefFile");
jpayne@69 827
jpayne@69 828 $OPT_UnQryFile = $OPT_Prefix . $OPT_UnQryFile;
jpayne@69 829 $TIGR->isCreatableFile("$OPT_UnQryFile")
jpayne@69 830 or $TIGR->isWritableFile("$OPT_UnQryFile")
jpayne@69 831 or push(@errs, "$OPT_UnQryFile");
jpayne@69 832
jpayne@69 833 if ( scalar(@errs) ) {
jpayne@69 834 print STDERR "ERROR: The following critical files could not be used\n";
jpayne@69 835 while ( scalar(@errs) ) { print(STDERR pop(@errs),"\n"); }
jpayne@69 836 print STDERR "Check your paths and file permissions and try again\n";
jpayne@69 837 exit(1);
jpayne@69 838 }
jpayne@69 839 }