convert-sam-for-rsem

   1 #!/usr/bin/perl
   2
   3 use Getopt::Long;
   4 use Pod::Usage;
   5 use strict;
   6
   7 my $standard_output = "&STDOUT";
   8
   9 my $out_file = $standard_output;
  10 my @tmp_dirs = ();
  11 my $help = 0;
  12
  13 GetOptions("o=s" => \$out_file,
  14            "T|temporary-directory=s" => \@tmp_dirs,
  15            "h|help" => \$help) or pd2usage(-exitval => 2, -verbose => 2);
  16
  17
  18 pod2usage(-verbose => 2) if ($help == 1);
  19 pod2usage(-msg => "Invalid number of arguments!", -exitval => 2, -verbose => 2) if (scalar(@ARGV) != 2);
  20
  21 my $command;
  22 my (@fields, @header) = ();
  23 my $M;
  24
  25 # Load fields
  26 my @lines = ();
  27 my $type;
  28
  29 open(INPUT, "$ARGV[0].ti");
  30 @lines = <INPUT>;
  31 chomp(@lines);
  32 close(INPUT);
  33
  34 @fields = ();
  35 ($M, $type) = split(/ /, $lines[0]);
  36 for (my $i = 0; $i < $M; $i++) {
  37     push(@fields, "SN:$lines[$i * 6 + 1]");
  38 }
  39
  40
  41 # Reorder header
  42 my $line;
  43
  44 open(INPUT, $ARGV[1]);
  45 @header = ();
  46 while (($line = <INPUT>) && substr($line, 0, 1) eq '@') {
  47     chomp($line);
  48     push(@header, $line);
  49 }
  50 close(INPUT);
  51
  52 my $n = scalar(@header);
  53 if ($n > 0) {
  54     my %hash = ();
  55     my @ktable = ();
  56
  57     my $tid = 0;
  58
  59     for (my $i = 0; $i < $n; $i++) {
  60         my @arr = split(/\t/, $header[$i]);
  61         if ($arr[0] ne "\@SQ") { push(@ktable, ""); next; }
  62         my $hasSN = 0;
  63         foreach my $key (@arr) {
  64             if (substr($key, 0, 3) eq "SN:") {
  65                 $hash{$key} = $i;
  66                 $hasSN = 1;
  67                 last;
  68             }
  69         }
  70         if (!$hasSN) { print STDERR "\"$header[$i]\" does not have a SN tag!\n"; exit(-1); }
  71         push(@ktable, $fields[$tid++]);
  72     }
  73
  74     if ($tid != $M) { print STDERR "Number of \@SQ lines is not correct!\n"; exit(-1); }
  75
  76     open(OUTPUT, ">$out_file");
  77     for (my $i = 0; $i < $n; $i++) {
  78         if ($ktable[$i] eq "") { print OUTPUT $header[$i]; }
  79         else { print OUTPUT $header[$hash{$ktable[$i]}]; }
  80         print OUTPUT "\n";
  81     }
  82     close(OUTPUT);
  83 }
  84
  85
  86 # extract alignment section
  87 $command = "grep ^[^@] $ARGV[1] > $ARGV[1].__temp";
  88 &runCommand($command);
  89
  90 # sort and output the alignment section
  91 $command = "sort -k 1,1 -s";
  92 if (scalar(@tmp_dirs) > 0) { $" = " -T "; $command .= " -T @tmp_dirs"; }
  93 $command .= " $ARGV[1].__temp";
  94 if ($out_file ne $standard_output) { $command .= " >> $out_file"; }
  95 &runCommand($command);
  96
  97 # delete temporary files
  98 $command = "rm -f $ARGV[1].__temp";
  99 &runCommand($command);
 100
 101 # finish
 102 print STDERR "Conversion is completed successfully!\n";
 103
 104 # command, {err_msg}
 105 sub runCommand {
 106     print STDERR $_[0]."\n";
 107     my $status = system($_[0]);
 108     if ($status != 0) {
 109         my $errmsg;
 110         if (scalar(@_) > 1) { $errmsg = $_[1]; }
 111         else { $errmsg = "\"$command\" failed! Plase check if you provide correct parameters/options for the script!"; }
 112         print STDERR $errmsg."\n";
 113         exit(-1);
 114     }
 115     print STDERR "\n";
 116 }
 117
 118 __END__
 119
 120 =head1 NAME
 121
 122 convert-sam-for-rsem
 123
 124 =head1 SYNOPSIS
 125
 126 =over
 127
 128  convert-sam-for-rsem [options] reference_name input_sam
 129
 130 =back
 131
 132 =head1 ARGUMENTS
 133
 134 =over
 135
 136 =item B<reference_name>
 137
 138 The name of the reference used. This should be the same name used by 'rsem-prepare-reference'.
 139
 140 =item B<input_sam>
 141
 142 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).
 143
 144 =back
 145
 146 =head1 OPTIONS
 147
 148 =over
 149
 150 =item B<-o> <file>
 151
 152 Output the converted SAM file into <file>. (Default: STDOUT)
 153
 154 =item B<-T/--temporary-directory> <directory>
 155
 156 '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".
 157
 158 =item B<-h/--help>
 159
 160 Show help information.
 161
 162 =back
 163
 164 =head1 DESCRIPTION
 165
 166 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.
 167
 168 Note: You do not need to run this script if Bowtie (not Bowtie 2) is used, or the order of @SQ lines is the same as 'reference_name.idx.fa' and the alignment lines of a same read group together and the mates of the same alignment are adjacent each other for paired-end reads.
 169
 170 Note: This program can only recognize SAM files. See ARGUMENTS section.
 171
 172 =head1 EXAMPLES
 173
 174 Suppose reference_name and input_sam are set to '/ref/mouse_125' and 'input.sam'.
 175
 176 1) Output to standard output and gzip the output to 'input_for_rsem.sam.gz':
 177
 178  convert-sam-for-rsem /ref/mouse_125 input.sam | gzip > input_for_rsem.sam.gz
 179
 180 2) Output to 'input_for_rsem.sam' directly:
 181
 182  convert-sam-for-rsem /ref/mouse_125 input.sam -o input_for_rsem.sam
 183
 184 =cut