convert-sam-for-rsem

   1 #!/usr/bin/perl
   2
   3 use Getopt::Long;
   4 use Pod::Usage;
   5 use File::Basename;
   6 use File::Path 'rmtree';
   7 use strict;
   8
   9 my ($in_file, $out_file) = ();
  10
  11 my @tmp_dirs = ();
  12 my $help = 0;
  13
  14 GetOptions("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
  23 $in_file = $ARGV[0];
  24 $out_file = "$ARGV[1].bam";
  25
  26 my ($fn, $dir, $suf) = fileparse($in_file, qr/\.[^.]*/);
  27
  28 $suf = lc(substr($suf, 1));
  29 pod2usage(-msg => "Input file's suffix is neither sam nor bam!", -exitval => 2, -verbose => 2) if (($suf ne "sam") && ($suf ne "bam"));
  30 my $isSam = ($suf eq "sam");
  31
  32 ($fn, $dir, $suf) = fileparse($0);
  33
  34 my $temp_dir = "$out_file.temp";
  35 if (-d $temp_dir) { print "Warning: $temp_dir exists, convert-sam-for-rsem will write temporary files into this folder and delete it after it finishes!\n"; }
  36 else {
  37     if (!mkdir($temp_dir)) { print "Fail to create folder $temp_dir.\n"; exit(-1); }
  38 }
  39
  40 # Convert to SAM format if input is a BAM file
  41
  42 my $sam_file;
  43
  44 if (!$isSam) {
  45     $sam_file = "$temp_dir/input.sam";
  46     $command = $dir."sam/samtools view -h -o $sam_file $in_file";
  47     &runCommand($command);
  48 }
  49 else {
  50     $sam_file = $in_file;
  51 }
  52
  53 # Phase I, sort entries so that all entries of a same read groups together
  54
  55 my $tmp_sam = "$temp_dir/tmp.sam";
  56
  57 # grep header section
  58 $command = "grep ^@ $sam_file > $tmp_sam";
  59 &runCommand($command);
  60
  61 # sort alignment section
  62 $command = "grep ^[^@] $sam_file | sort -k 1,1 -s";
  63 if (scalar(@tmp_dirs) > 0) { $" = " -T "; $command .= " -T @tmp_dirs"; }
  64 $command .= " >> $tmp_sam";
  65 &runCommand($command);
  66
  67 # Phase II, parse the temporary SAM file to make paired-end alignments' two mates adjacent to each other
  68
  69 $command = $dir."rsem-scan-for-paired-end-reads $tmp_sam $out_file";
  70 &runCommand($command);
  71
  72 # delete temporary directory
  73 rmtree($temp_dir);
  74
  75 print STDERR "Conversion is completed. $out_file will be checked by 'rsem-sam-validator'.\n";
  76
  77 # Phase III, validate if the resulting bam file is correct
  78
  79 $command = $dir."rsem-sam-validator $out_file";
  80 &runCommand($command);
  81
  82
  83 # command, {err_msg}
  84 sub runCommand {
  85     print $_[0]."\n";
  86     my $status = system($_[0]);
  87     if ($status != 0) {
  88         my $errmsg;
  89         if (scalar(@_) > 1) { $errmsg = $_[1]; }
  90         else { $errmsg = "\"$command\" failed! Plase check if you provide correct parameters/options for the pipeline!"; }
  91         print $errmsg."\n";
  92         exit(-1);
  93     }
  94     print "\n";
  95 }
  96
  97 __END__
  98
  99 =head1 NAME
 100
 101 convert-sam-for-rsem
 102
 103 =head1 SYNOPSIS
 104
 105 =over
 106
 107  convert-sam-for-rsem [options] <input.sam/input.bam> output_file_name
 108
 109 =back
 110
 111 =head1 ARGUMENTS
 112
 113 =over
 114
 115 =item B<input.sam/input.bam>
 116
 117 The SAM or BAM file generated by user's aligner. We require this file contains the header section. If input is a SAM file, it must end with suffix 'sam' (case insensitive). If input is a BAM file, it must end with suffix 'bam' (case insensitive).
 118
 119 =item B<output_file_name>
 120
 121 The output name for the converted file. 'convert-sam-for-rsem' will output a BAM with the name 'output_file_name.bam'.
 122
 123 =back
 124
 125 =head1 OPTIONS
 126
 127 =over
 128
 129 =item B<-T/--temporary-directory> <directory>
 130
 131 '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".
 132
 133 =item B<-h/--help>
 134
 135 Show help information.
 136
 137 =back
 138
 139 =head1 DESCRIPTION
 140
 141 This program converts the SAM/BAM file generated by user's aligner into a BAM 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 and output header sections. This program will create a temporary directory called 'output_file_name.bam.temp' to store the intermediate files. The directory will be deleted automatically after the conversion. After the conversion, this program will call 'rsem-sam-validator' to validate the resulting BAM file.
 142
 143 Note: You do not need to run this script if `rsem-sam-validator' reports that your SAM/BAM file is valid.
 144
 145 Note: This program does not check the correctness of input file. You should make sure the input is a valid SAM/BAM format file.
 146
 147 =head1 EXAMPLES
 148
 149 Suppose input is set to 'input.sam' and output file name is "output"
 150
 151  convert-sam-for-rsem input.sam output
 152
 153 We will get a file called 'output.bam' as output.
 154
 155 =cut