convert-sam-for-rsem

   1 #!/usr/bin/perl
   2
   3 use Getopt::Long;
   4 use Pod::Usage;
   5 use strict;
   6
   7 my $out_file = "";
   8 my @tmp_dirs = ();
   9 my $help = 0;
  10
  11 GetOptions("o=s" => \$out_file,
  12            "T|temporary-directory=s" => \@tmp_dirs,
  13            "h|help" => \$help) or pd2usage(-exitval => 2, -verbose => 2);
  14
  15
  16 pod2usage(-verbose => 2) if ($help == 1);
  17 pod2usage(-msg => "Invalid number of arguments!", -exitval => 2, -verbose => 2) if (scalar(@ARGV) != 1);
  18
  19 my $command;
  20
  21 # grep header section
  22 $command = "grep ^@ $ARGV[0]";
  23 if ($out_file ne "") { $command .= " > $out_file"; }
  24 &runCommand($command);
  25
  26 # sort alignment section
  27 $command = "grep ^[^@] $ARGV[0] | sort -k 1,1 -s";
  28 if (scalar(@tmp_dirs) > 0) { $" = " -T "; $command .= " -T @tmp_dirs"; }
  29 if ($out_file ne "") { $command .= " >> $out_file"; }
  30 &runCommand($command);
  31
  32 # finish
  33 print STDERR "Conversion is completed successfully!\n";
  34
  35 # command, {err_msg}
  36 sub runCommand {
  37     print STDERR $_[0]."\n";
  38     my $status = system($_[0]);
  39     if ($status != 0) {
  40         my $errmsg;
  41         if (scalar(@_) > 1) { $errmsg = $_[1]; }
  42         else { $errmsg = "\"$command\" failed! Plase check if you provide correct parameters/options for the script!"; }
  43         print STDERR $errmsg."\n";
  44         exit(-1);
  45     }
  46     print STDERR "\n";
  47 }
  48
  49 __END__
  50
  51 =head1 NAME
  52
  53 convert-sam-for-rsem
  54
  55 =head1 SYNOPSIS
  56
  57 =over
  58
  59  convert-sam-for-rsem [options] input_sam
  60
  61 =back
  62
  63 =head1 ARGUMENTS
  64
  65 =over
  66
  67 =item B<input_sam>
  68
  69 The SAM file (*.sam) generated by user's aligner. If the aligner produces a BAM file, please use samtools to convert it to a SAM file (with header information).
  70
  71 =back
  72
  73 =head1 OPTIONS
  74
  75 =over
  76
  77 =item B<-o> <file>
  78
  79 Output the converted SAM file into <file>. (Default: STDOUT)
  80
  81 =item B<-T/--temporary-directory> <directory>
  82
  83 'convert-sam-for-rsem' will call 'sort' command and this is the '-T/--temporary-directory' option of 'sort' command. The following is the description from 'sort' : "use DIR for temporaries, not $TMPDIR or /tmp; multiple options specify multiple directories".
  84
  85 =item B<-h/--help>
  86
  87 Show help information.
  88
  89 =back
  90
  91 =head1 DESCRIPTION
  92
  93 This program converts the SAM file generated by user's aligner into a SAM file which RSEM can process. However, users should make sure their aligners use 'reference_name.idx.fa' generated by 'rsem-prepare-reference' as their references. In addition, their aligners should output header information and make two mates of the same alignment adjacent to each other for paired-end data. This program will output the converted file into standard output by default for the purpose of piping. By setting '-o' option, users can make the converted file written into disk.
  94
  95 Note: You do not need to run this script if Bowtie (not Bowtie 2) is used, or the alignment lines of a same read group together and the mates of the same alignment are adjacent each other for paired-end reads.
  96
  97 Note: This program can only recognize SAM files. See ARGUMENTS section.
  98
  99 =head1 EXAMPLES
 100
 101 Suppose input_sam is set to 'input.sam'.
 102
 103 1) Output to standard output and gzip the output to 'input_for_rsem.sam.gz':
 104
 105  convert-sam-for-rsem input.sam | gzip > input_for_rsem.sam.gz
 106
 107 2) Output to 'input_for_rsem.sam' directly:
 108
 109  convert-sam-for-rsem input.sam -o input_for_rsem.sam
 110
 111 =cut