rsem-calculate-expression

   1 #!/usr/bin/perl
   2
   3 use Getopt::Long;
   4 use Pod::Usage;
   5 use File::Basename;
   6 use Switch;
   7 use strict;
   8
   9 #const
  10 my $BURNIN = 200;
  11 my $CHAINLEN = 1000;
  12 my $SAMPLEGAP = 1;
  13 my $CONFIDENCE = 0.95;
  14 my $NSPC = 50;
  15
  16 my $NMB = 1024; # default
  17
  18 my $status = 0;
  19
  20 my $read_type = 1; # default, single end with qual
  21
  22 my $bowtie_path = "";
  23 my $C = 2;
  24 my $E = 99999999;
  25 my $L = 25;
  26 my $maxHits = 200;
  27 my $phred33 = 0;
  28 my $phred64 = 0;
  29 my $solexa = 0;
  30
  31 my $is_sam = 0;
  32 my $is_bam = 0;
  33 my $fn_list = "";
  34 my $tagName = "XM";
  35
  36 my $probF = 0.5;
  37
  38 my $minL = 1;
  39 my $maxL = 1000;
  40 my $mean = -1;
  41 my $sd = 0;
  42
  43 my $estRSPD = 0;
  44 my $B = 20;
  45
  46 my $nThreads = 1;
  47 my $genBamF = 0;
  48 my $calcCI = 0;
  49 my $quiet = 0;
  50 my $help = 0;
  51
  52 my $paired_end = 0;
  53 my $no_qual = 0;
  54 my $keep_intermediate_files = 0;
  55
  56 my $strand_specific = 0;
  57
  58 my $mTime = 0;
  59 my ($time_start, $time_end, $time_alignment, $time_rsem, $time_ci) = (0, 0, 0, 0, 0);
  60
  61 GetOptions("keep-intermediate-files" => \$keep_intermediate_files,
  62            "no-qualities" => \$no_qual,
  63            "paired-end" => \$paired_end,
  64            "strand-specific" => \$strand_specific,
  65            "sam" => \$is_sam,
  66            "bam" => \$is_bam,
  67            "sam-header-info=s" => \$fn_list,
  68            "tag=s" => \$tagName,
  69            "seed-length=i" => \$L,
  70            "bowtie-path=s" => \$bowtie_path,
  71            "bowtie-n=i" => \$C,
  72            "bowtie-e=i" => \$E,
  73            "bowtie-m=i" => \$maxHits,
  74            "phred33-quals" => \$phred33,
  75            "phred64-quals" => \$phred64, #solexa1.3-quals" => \$phred64,
  76            "solexa-quals" => \$solexa,
  77            "forward-prob=f" => \$probF,
  78            "fragment-length-min=i" => \$minL,
  79            "fragment-length-max=i" => \$maxL,
  80            "fragment-length-mean=f" => \$mean,
  81            "fragment-length-sd=f" => \$sd,
  82            "estimate-rspd" => \$estRSPD,
  83            "num-rspd-bins=i" => \$B,
  84            "p|num-threads=i" => \$nThreads,
  85            "out-bam" => \$genBamF,
  86            "calc-ci" => \$calcCI,
  87            "ci-memory=i" => \$NMB,
  88            "time" => \$mTime,
  89            "q|quiet" => \$quiet,
  90            "h|help" => \$help) or pod2usage(-exitval => 2, -verbose => 2);
  91
  92 pod2usage(-verbose => 2) if ($help == 1);
  93
  94
  95 #check parameters and options
  96
  97 if ($is_sam || $is_bam) {
  98     pod2usage(-msg => "Invalid number of arguments!", -exitval => 2, -verbose => 2) if (scalar(@ARGV) != 3);
  99     pod2usage(-msg => "--sam and --bam cannot be active at the same time!", -exitval => 2, -verbose => 2) if ($is_sam == 1&& $is_bam == 1);
 100     pod2usage(-msg => "--bowtie-path, --bowtie-n, --bowtie-e, --bowtie-m, --phred33-quals, --phred64-quals or --solexa-quals cannot be set if input is SAM/BAM format!", -exitval => 2, -verbose => 2) if ($bowtie_path ne "" || $C != 2 || $E != 99999999 || $maxHits != 200 || $phred33 || $phred64 || $solexa);
 101 }
 102 else {
 103     pod2usage(-msg => "Invalid number of arguments!", -exitval => 2, -verbose => 2) if (!$paired_end && scalar(@ARGV) != 3 || $paired_end && scalar(@ARGV) != 4);
 104     pod2usage(-msg => "Only one of --phred33-quals --phred64-quals/--solexa1.3-quals --solexa-suqls can be active!", -exitval => 2, -verbose => 2) if ($phred33 + $phred64 + $solexa > 1);
 105     podwusage(-msg => "--sam , --bam or --sam-header-info cannot be set if use bowtie aligner to produce alignments!", -exitval => 2, -verbose => 2) if ($is_sam || $is_bam || $fn_list ne "");
 106 }
 107
 108 pod2usage(-msg => "Forward probability should be in [0, 1]!", -exitval => 2, -verbose => 2) if ($probF < 0 || $probF > 1);
 109 pod2usage(-msg => "Min fragment length should be at least 1!", -exitval => 2, -verbose => 2) if ($minL < 1);
 110 pod2usage(-msg => "Min fragment length should be smaller or equal to max fragment length!", -exitval => 2, -verbose => 2) if ($minL > $maxL);
 111 pod2usage(-msg => "The memory allocated for calculating credibility intervals should be at least 1 MB!\n", -exitval => 2, -verbose => 2) if ($NMB < 1);
 112 pod2usage(-msg => "Number of threads should be at least 1!\n", -exitval => 2, -verbose => 2) if ($nThreads < 1);
 113 pod2usage(-msg => "Seed length should be at least 25!\n", -exitval => 2, -verbose => 2) if ($L < 25);
 114
 115 if ($strand_specific) { $probF = 1.0; }
 116
 117 my $mate1_list = "";
 118 my $mate2_list = "";
 119 my $inpF = "";
 120
 121 my ($refName, $sampleName, $sampleToken, $temp_dir, $stat_dir, $imdName) = ();
 122 my $gap = 32;
 123
 124 if ($paired_end) {
 125     if ($no_qual) { $read_type = 2; }
 126     else { $read_type = 3; }
 127 }
 128 else {
 129     if ($no_qual) { $read_type = 0; }
 130     else { $read_type = 1; }
 131 }
 132
 133 if (scalar(@ARGV) == 3) {
 134     if ($is_sam || $is_bam) { $inpF = $ARGV[0]; }
 135     else {$mate1_list = $ARGV[0]; }
 136     $refName = $ARGV[1];
 137     $sampleName = $ARGV[2];
 138 }
 139 else {
 140     $mate1_list = $ARGV[0];
 141     $mate2_list = $ARGV[1];
 142     $refName = $ARGV[2];
 143     $sampleName = $ARGV[3];
 144 }
 145
 146 my $pos = rindex($sampleName, '/');
 147 if ($pos < 0) { $sampleToken = $sampleName; }
 148 else { $sampleToken = substr($sampleName, $pos + 1); }
 149
 150 $temp_dir = "$sampleName.temp";
 151 $stat_dir = "$sampleName.stat";
 152
 153 if (!(-d $temp_dir) && !mkdir($temp_dir)) { print "Fail to create folder $temp_dir.\n"; exit(-1); }
 154 if (!(-d $stat_dir) && !mkdir($stat_dir)) { print "Fail to create folder $stat_dir.\n"; exit(-1); }
 155
 156 $imdName = "$temp_dir/$sampleToken";
 157
 158 if (!$is_sam && !$is_bam && $phred33 + $phred64 + $solexa == 0) { $phred33 = 1; }
 159
 160 my ($mate_minL, $mate_maxL) = (1, $maxL);
 161
 162 if ($bowtie_path ne "") { $bowtie_path .= "/"; }
 163
 164 my ($fn, $dir, $suf) = fileparse($0);
 165 my $command = "";
 166
 167 if (!$is_sam && !$is_bam) {
 168     $command = $bowtie_path."bowtie";
 169     if ($read_type == 0 || $read_type == 2) { $command .= " -f"; }
 170     else { $command .= " -q"; }
 171
 172     if ($phred33) { $command .= " --phred33-quals"; }
 173     elsif ($phred64) { $command .= " --phred64-quals"; }
 174     elsif ($solexa) { $command .= " --solexa-quals"; }
 175     else { print "Oh, no!!!"; exit(2); }
 176
 177     $command .= " -n $C -e $E -l $L";
 178
 179     if ($read_type == 2 || $read_type == 3) { $command .= " -I $minL -X $maxL"; }
 180
 181     if ($strand_specific || $probF == 1.0) { $command .= " --norc"; }
 182     elsif ($probF == 0.0) { $command .= " --nofw"; }
 183
 184     $command .= " -p $nThreads -a -m $maxHits -S";
 185     if ($quiet) { $command .= " --quiet"; }
 186
 187     $command .= " $refName";
 188     if ($read_type == 0 || $read_type == 1) {
 189         $command .= " $mate1_list";
 190     }
 191     else {
 192         $command .= " -1 $mate1_list -2 $mate2_list";
 193     }
 194
 195     $command .= " | gzip > $imdName.sam.gz";
 196     print "$command\n";
 197
 198     if ($mTime) { $time_start = time(); }
 199
 200     $status = system($command);
 201
 202     if ($mTime) { $time_end = time(); $time_alignment = $time_end - $time_start; }
 203
 204     if ($status != 0) {
 205         print "bowtie failed! Please check if you provide correct parameters/options for the pipeline!\n";
 206         exit(-1);
 207     }
 208     print "\n";
 209
 210     $inpF = "$imdName.sam.gz";
 211     $is_sam = 1; # output of bowtie is a sam file
 212 }
 213
 214 if ($mTime) { $time_start = time(); }
 215
 216 $command = $dir."rsem-parse-alignments $refName $sampleName $sampleToken";
 217
 218 my $samInpType;
 219 if ($is_sam) { $samInpType = "s"; }
 220 elsif ($is_bam) { $samInpType = "b"; }
 221
 222 $command .= " $samInpType $inpF -t $read_type";
 223 if ($fn_list ne "") { $command .= " -l $fn_list"; }
 224 if ($tagName ne "") { $command .= " -tag $tagName"; }
 225 if ($quiet) { $command .= " -q"; }
 226
 227 print "$command\n";
 228 $status = system($command);
 229 if ($status != 0) {
 230     print "rsem-parse-alignments failed! Please check if you provide correct parameters/options for the pipeline!\n";
 231     exit(-1);
 232 }
 233 print "\n";
 234
 235 $command = $dir."rsem-build-read-index $gap";
 236 switch($read_type) {
 237     case 0  { $command .= " 0 $quiet $imdName\_alignable.fa"; }
 238     case 1  { $command .= " 1 $quiet $imdName\_alignable.fq"; }
 239     case 2  { $command .= " 0 $quiet $imdName\_alignable_1.fa $imdName\_alignable_2.fa"; }
 240     case 3  { $command .= " 1 $quiet $imdName\_alignable_1.fq $imdName\_alignable_2.fq"; }
 241 }
 242 print "$command\n";
 243 $status = system($command);
 244 if ($status != 0) {
 245     print "rsem-build-read-index failed! Please check if you provide correct parameters/options for the pipeline!\n";
 246     exit(-1);
 247 }
 248 print "\n";
 249
 250 my $doesOpen = open(OUTPUT, ">$imdName.mparams");
 251 if ($doesOpen == 0) { print "Cannot generate $imdName.mparams!\n"; exit(-1); }
 252 print OUTPUT "$minL $maxL\n";
 253 print OUTPUT "$probF\n";
 254 print OUTPUT "$estRSPD\n";
 255 print OUTPUT "$B\n";
 256 print OUTPUT "$mate_minL $mate_maxL\n";
 257 print OUTPUT "$mean $sd\n";
 258 print OUTPUT "$L\n";
 259 close(OUTPUT);
 260
 261 $command = $dir."rsem-run-em $refName $read_type $sampleName $sampleToken -p $nThreads";
 262 if ($genBamF) {
 263     $command .= " -b $samInpType $inpF";
 264     if ($fn_list ne "") { $command .= " 1 $fn_list"; }
 265     else { $command .= " 0"; }
 266 }
 267 if ($calcCI) { $command .= " --gibbs-out"; }
 268 if ($quiet) { $command .= " -q"; }
 269
 270 print "$command\n";
 271 $status = system($command);
 272 if ($status != 0) {
 273     print "rsem-run-em failed! Please check if you provide correct parameters/options for the pipeline!\n";
 274     exit(-1);
 275 }
 276 print "\n";
 277
 278 if ($genBamF) {
 279     $command = $dir."sam/samtools sort $sampleName.bam $sampleName.sorted";
 280     print "$command\n";
 281     $status = system($command);
 282     if ($status != 0) {
 283         print "sam/samtools sort failed! Please check if you provide correct parameters/options for the pipeline!\n";
 284         exit(-1);
 285     }
 286     print "\n";
 287     $command = $dir."sam/samtools index $sampleName.sorted.bam";
 288     print "$command\n";
 289     $status = system($command);
 290     if ($status != 0) {
 291         print "sam/samtools index failed! Please check if you provide correct parameters/options for the pipeline!\n";
 292         exit(-1);
 293     }
 294     print "\n";
 295 }
 296
 297 &collectResults("$imdName.iso_res", "$sampleName.isoforms.results"); # isoform level
 298 &collectResults("$imdName.gene_res", "$sampleName.genes.results"); # gene level
 299
 300 if ($mTime) { $time_end = time(); $time_rsem = $time_end - $time_start; }
 301
 302 if ($mTime) { $time_start = time(); }
 303
 304 if ($calcCI) {
 305     $command = $dir."rsem-run-gibbs $refName $sampleName $sampleToken $BURNIN $CHAINLEN $SAMPLEGAP";
 306 #    $command .= " -p $nThreads";
 307     if ($quiet) { $command .= " -q"; }
 308     print "$command\n";
 309     $status = system($command);
 310     if ($status != 0) {
 311         print "rsem-run-gibbs failed! Please check if you provide correct parameters/options for the pipeline!\n";
 312         exit(-1);
 313     }
 314     print "\n";
 315
 316     system("mv $sampleName.isoforms.results $imdName.isoforms.results.bak1");
 317     system("mv $sampleName.genes.results $imdName.genes.results.bak1");
 318     &collectResults("$imdName.iso_res", "$sampleName.isoforms.results"); # isoform level
 319     &collectResults("$imdName.gene_res", "$sampleName.genes.results"); # gene level
 320
 321     $command = $dir."rsem-calculate-credibility-intervals $refName $sampleName $sampleToken $CONFIDENCE $NSPC $NMB";
 322     if ($quiet) { $command .= " -q"; }
 323     print "$command\n";
 324     $status = system($command);
 325     if ($status != 0) {
 326         print "rsem-calculate-credibility-intervals failed! Please check if you provide correct parameters/options for the pipeline!\n";
 327         exit(-1);
 328     }
 329     print "\n";
 330
 331     system("mv $sampleName.isoforms.results $imdName.isoforms.results.bak2");
 332     system("mv $sampleName.genes.results $imdName.genes.results.bak2");
 333     &collectResults("$imdName.iso_res", "$sampleName.isoforms.results"); # isoform level
 334     &collectResults("$imdName.gene_res", "$sampleName.genes.results"); # gene level
 335 }
 336
 337 if ($mTime) { $time_end = time(); $time_ci = $time_end - $time_start; }
 338
 339 if ($mTime) { $time_start = time(); }
 340
 341 if (!$keep_intermediate_files) {
 342     $status = system("rm -rf $temp_dir");
 343     if ($status != 0) {
 344         print "Fail to delete the temporary folder!\n";
 345         exit(-1);
 346     }
 347 }
 348
 349 if ($mTime) { $time_end = time(); }
 350
 351 if ($mTime) {
 352     open(OUTPUT, ">$sampleName.time");
 353     print OUTPUT "Alignment: $time_alignment s.\n";
 354     print OUTPUT "RSEM: $time_rsem s.\n";
 355     print OUTPUT "CI: $time_ci s.\n";
 356     my $time_del = $time_end - $time_start;
 357     print OUTPUT "Delete: $time_del s.\n";
 358     close(OUTPUT);
 359 }
 360
 361 # inpF, outF
 362 sub collectResults {
 363     my $local_status;
 364     my ($inpF, $outF);
 365     my (@results, @comment) = ();
 366     my $line;
 367     my $cnt;
 368
 369     $inpF = $_[0];
 370     $outF = $_[1];
 371
 372     $local_status = open(INPUT, $inpF);
 373     if ($local_status == 0) { print "Fail to open file $inpF!\n"; exit(-1); }
 374
 375     $cnt = 0;
 376     @results = ();
 377
 378     while ($line = <INPUT>) {
 379         ++$cnt;
 380         chomp($line);
 381         my @local_arr = split(/\t/, $line);
 382         if ($cnt == 4) { @comment = @local_arr; }
 383         else { push(@results, \@local_arr); }
 384     }
 385
 386     push(@results, \@comment);
 387     close(INPUT);
 388
 389     $local_status = open(OUTPUT, ">$outF");
 390     if ($local_status == 0) { print "Fail to create file $outF!\n"; exit(-1); }
 391
 392     my $n = scalar(@results);
 393     my $m = scalar(@{$results[0]});
 394     for (my $i = 0; $i < $m; $i++) {
 395         my @out_arr = ();
 396         for (my $j = 0; $j < $n; $j++) { push(@out_arr, $results[$j][$i]); }
 397         $" = "\t";
 398         print OUTPUT "@out_arr\n";
 399     }
 400     close(OUTPUT);
 401 }
 402
 403
 404 __END__
 405
 406 =head1 NAME
 407
 408 rsem-calculate-expression
 409
 410 =head1 SYNOPSIS
 411
 412 =over
 413
 414  rsem-calculate-expression [options] upstream_read_file(s) reference_name sample_name
 415  rsem-calculate-expression [options] --paired-end upstream_read_file(s) downstream_read_file(s) reference_name sample_name
 416  rsem-calculate-expression [options] --sam/--bam [--paired-end] input reference_name sample_name
 417
 418 =back
 419
 420 =head1 ARGUMENTS
 421
 422 =over
 423
 424 =item B<upstream_read_files(s)>
 425
 426 Comma-separated list of files containing single-end reads or upstream reads for paired-end data.  By default, these files are assumed to be in FASTQ format.  If the --no-qualities option is specified, then FASTA format is expected.
 427
 428 =item B<downstream_read_file(s)>
 429
 430 Comma-separated list of files containing downstream reads which are paired with the upstream reads.  By default, these files are assumed to be in FASTQ format.  If the --no-qualities option is specified, then FASTA format is expected.
 431
 432 =item B<input>
 433
 434 SAM/BAM formatted input file.  If "-" is specified for the filename, SAM/BAM input is instead assumed to come from standard input. RSEM requires all alignments of the same read group together. For paired-end reads, RSEM also requires the two mates of any alignment be adjacent. See Description section for how to make input file obey RSEM's requirements.
 435
 436 =item B<reference_name>
 437
 438 The name of the reference used.  The user must have run 'rsem-prepare-reference' with this reference_name before running this program.
 439
 440 =item B<sample_name>
 441
 442 The name of the sample analyzed. All output files are prefixed by this name (e.g., sample_name.genes.results)
 443
 444 =back
 445
 446 =head1 OPTIONS
 447
 448 =over
 449
 450 =item B<--paired-end>
 451
 452 Input reads are paired-end reads. (Default: off)
 453
 454 =item B<--no-qualities>
 455
 456 Input reads do not contain quality scores. (Default: off)
 457
 458 =item B<--strand-specific>
 459
 460 The RNA-Seq protocol used to generate the reads is strand specific, i.e., all (upstream) reads are derived from the forward strand.  This option is equivalent to --forward-prob=1.0.  With this option set, if RSEM runs the Bowtie aligner, the '--norc' Bowtie option will be used, which disables alignment to the reverse strand of transcripts.  (Default: off)
 461
 462 =item B<--sam>
 463
 464 Input file is in SAM format. (Default: off)
 465
 466 =item B<--bam>
 467
 468 Input file is in BAM format. (Default: off)
 469
 470 =item B<--sam-header-info> <file>
 471
 472 RSEM reads header information from input by default. If this option is on, header information is read from the specified file. For the format of the file, please see SAM official website. (Default: "")
 473
 474 =item B<-p/--num-threads> <int>
 475
 476 Number of threads to use. Both Bowtie and expression estimation will use this many threads. (Default: 1)
 477
 478 =item B<--out-bam>
 479
 480 Generate a BAM file, 'sample_name.bam', with alignments mapped to genomic coordinates and annotated with their posterior probabilities. In addition, RSEM will call samtools (included in RSEM package) to sort and index the bam file. 'sample_name.sorted.bam' and 'sample_name.sorted.bam.bai' will be generated. (Default: off)
 481
 482 =item B<--calc-ci>
 483
 484 Calculate 95% credibility intervals and posterior mean estimates.  (Default: off)
 485
 486 =item B<--seed-length> <int>
 487
 488 Seed length used by the read aligner.  Providing the correct value for this parameter is important for RSEM's accuracy if the data are single-end reads.  If RSEM runs Bowtie, it uses this value for Bowtie's seed length parameter. The minimum value is 25. Any read with its or at least one of its mates' (for paired-end reads) length less than 25 will be ignored. (Default: 25)
 489
 490 =item B<--tag> <string>
 491
 492 The name of the optional field used in the SAM input for identifying a read with too many valid alignments. The field should have the format <tagName>:i:<value>, where a <value> bigger than 0 indicates a read with too many alignments. (Default: "")
 493
 494 =item B<--bowtie-path> <path>
 495
 496 The path to the bowtie executables. (Default: the path to the bowtie executables is assumed to be in the user's PATH environment variable)
 497
 498 =item B<--bowtie-n> <int>
 499
 500 (Bowtie parameter) max # of mismatches in the seed. (Range: 0-3, Default: 2)
 501
 502 =item B<--bowtie-e> <int>
 503
 504 (Bowtie parameter) max sum of mismatch quality scores across the alignment. (Default: 99999999)
 505
 506 =item B<--bowtie-m> <int>
 507
 508 (Bowtie parameter) suppress all alignments for a read if > <int> valid alignments exist. (Default: 200)
 509
 510 =item B<--phred33-quals>
 511
 512 Input quality scores are encoded as Phred+33. (Default: on)
 513
 514 =item B<--phred64-quals>
 515
 516 Input quality scores are encoded as Phred+64 (default for GA Pipeline ver. >= 1.3). (Default: off)
 517
 518 =item B<--solexa-quals>
 519
 520 Input quality scores are solexa encoded (from GA Pipeline ver. < 1.3). (Default: off)
 521
 522 =item B<--forward-prob> <double>
 523
 524 Probability of generating a read from the forward strand of a transcript. Set to 1 for a strand-specific protocol where all (upstream) reads are derived from the forward strand, 0 for a strand-specific protocol where all (upstream) read are derived from the reverse strand, or 0.5 for a non-strand-specific protocol. (Default: 0.5)
 525
 526 =item B<--fragment-length-min> <int>
 527
 528 Minimum read/insert length allowed. This is also the value for the bowtie -I option. (Default: 1)
 529
 530 =item B<--fragment-length-max> <int>
 531
 532 Maximum read/insert length allowed. This is also the value for the bowtie -X option. (Default: 1000)
 533
 534 =item B<--fragment-length-mean> <double>
 535
 536 (single-end data only) The mean of the fragment length distribution, which is assumed to be a Gaussian. (Default: -1, which disables use of the fragment length distribution)
 537
 538 =item B<--fragment-length-sd> <double>
 539
 540 (single-end data only) The standard deviation of the fragment length distribution, which is assumed to be a Gaussian.  (Default: 0, which assumes that all fragments are of the same length, given by the rounded value of B<--fragment-length-mean>)
 541
 542 =item B<--estimate-rspd>
 543
 544 Set this option if you want to estimate the read start position distribution (RSPD) from data. Otherwise, RSEM will use a uniform RSPD. (Default: off)
 545
 546 =item B<--num-rspd-bins> <int>
 547
 548 Number of bins in the RSPD. Only relevant when '--estimate-rspd' is specified.  Use of the default setting is recommended. (Default: 20)
 549
 550 =item B<--ci-memory> <int>
 551
 552 Amount of memory (in MB) RSEM is allowed to use for computing credibility intervals. (Default: 1024)
 553
 554 =item B<--keep-intermediate-files>
 555
 556 Keep temporary files generated by RSEM.  RSEM creates a temporary directory, 'sample_name.temp', into which it puts all intermediate output files. If this directory already exists, RSEM overwrites all files generated by previous RSEM runs inside of it. By default, after RSEM finishes, the temporary directory is deleted.  Set this option to prevent the deletion of this directory and the intermediate files inside of it. (Default: off)
 557
 558 =item B<-q/--quiet>
 559
 560 Suppress the output of logging information. (Default: off)
 561
 562 =item B<-h/--help>
 563
 564 Show help information.
 565
 566 =back
 567
 568 =head1 DESCRIPTION
 569
 570 In its default mode, this program aligns input reads against a reference transcriptome with Bowtie and calculates expression values using the alignments.  RSEM assumes the data are single-end reads with quality scores, unless the '--paired-end' or '--no-qualities' options are specified.  Users may use an alternative aligner by specifying one of the --sam and --bam options, and providing an alignment file in the specified format. However, users should make sure the alignment file satisfies the requirements mentioned in ARGUMENTS section.
 571
 572 One simple way to make the alignment file (e.g. input.sam) satisfying RSEM's requirements (assuming the aligner used put mates in a paired-end read adjacent) is to use the following command:
 573
 574   sort -k 1,1 -s input.sam > input.sorted.sam
 575
 576 The SAM/BAM format RSEM uses is v1.3. However, it is compatible with old SAM/BAM format.
 577
 578 The user must run 'rsem-prepare-reference' with the appropriate reference before using this program.
 579
 580 For single-end data, it is strongly recommended that the user provide the fragment length distribution parameters (--fragment-length-mean and --fragment-length-sd).  For paired-end data, RSEM will automatically learn a fragment length distribution from the data.
 581
 582 Please note that some of the default values for the Bowtie parameters are not the same as those defined for Bowtie itself.
 583
 584 The temporary directory and all intermediate files will be removed when RSEM finishes unless '--keep-intermediate-files' is specified.
 585
 586 With the "--calc-ci" option, 95% credibility intervals and posterior mean estimates will be calculated in addition to maximum likelihood estimates.
 587
 588 =head1 OUTPUT
 589
 590 =over
 591
 592 =item B<sample_name.genes.results>
 593
 594 File containing gene level expression estimates. The format of each
 595 line in this file is:
 596
 597 gene_id expected_counts tau_value [pmc_value tau_pme_value tau_ci_lower_bound tau_ci_upper_bound] transcript_id_list
 598
 599 Fields are separated by the tab character. Fields within "[]" are only
 600 presented if '--calc-ci' is set. pme stands for posterior mean
 601 estimation. pmc stands for posterior mean counts. ci_lower_bound(l)
 602 means the lower bound of the credibility intervals, ci_upper_bound(u)
 603 means the upper bound of the credibility intervals. So the credibility
 604 interval is [l, u]. 'transcript_id_list' is a space-separated list of
 605 transcript_ids belonging to the gene.
 606
 607 =item B<sample_name.isoforms.results>
 608
 609 File containing isoform level expression values. The format of each
 610 line in this file is:
 611
 612 transcript_id expected_counts tau_value [pmc_value tau_pme_value tau_ci_lower_bound tau_ci_upper_bound] other_attributes
 613
 614 Fields are separated by the tab character. 'other_attributes' are all
 615 other attributes after attribute 'transcript_id' field in the GTF
 616 file. If no other attributes are given or no GTF file is provided in
 617 'rsem-prepare-reference', there will be no tab after the
 618 tau_value field.
 619
 620 =item B<sample_name.bam, sample_name.sorted.bam and sample_name.sorted.bam.bai>
 621
 622 Only generated when --out-bam is specified.
 623
 624 'sample_name.bam' is a BAM-formatted file of read alignments in
 625 genomic coordinates. Alignments of reads that have identical genomic
 626 coordinates (i.e., alignments to different isoforms that share the
 627 same genomic region) are collapsed into one alignment.  The MAPQ field
 628 of each alignment is set to min(100, floor(-10 * log10(1.0 - w) +
 629 0.5)), where w is the posterior probability of that alignment being
 630 the true mapping of a read.  In addition, RSEM pads a new tag
 631 ZW:f:value, where value is a single precision floating number
 632 representing the posterior probability.
 633
 634 'sample_name.sorted.bam' and 'sample_name.sorted.bam.bai' are the
 635 sorted BAM file and indices generated by samtools (included in RSEM package).
 636
 637 =item B<sample_name.stat>
 638
 639 This is a folder instead of a file. All model related statistics are stored in this folder. Use 'rsem-plot-model' can generate plots using this folder.
 640
 641 =back
 642
 643 =head1 EXAMPLES
 644
 645 Assume the path to the bowtie executables is in the user's PATH environment variable. Reference files are under '/ref' with name 'mm9'.
 646
 647 1) '/data/mmliver.fq', single-end reads with quality scores. Quality scores are encoded as for 'GA pipeline version >= 1.3'. We want to use 8 threads and generate a BAM file:
 648
 649  rsem-calculate-expression --phred64-quals \
 650                            -p 8 \
 651                            --out-bam \
 652                            /data/mmliver.fq \
 653                            /ref/mm9 \
 654                            mmliver_single_quals
 655
 656 2) '/data/mmliver_1.fq' and '/data/mmliver_2.fq', paired-end reads with quality scores. Quality scores are in SANGER format. We want to use 8 threads and do not generate a BAM file:
 657
 658  rsem-calculate-expression -p 8 \
 659                            --paired-end \
 660                            /data/mmliver_1.fq \
 661                            /data/mmliver_2.fq \
 662                            /ref/mm9 \
 663                            mmliver_paired_end_quals
 664
 665 3) '/data/mmliver.fa', single-end reads without quality scores. We want to use 8 threads and generate a BAM file:
 666
 667  rsem-calculate-expression -p 8 \
 668                            --no-qualities \
 669                            /data/mmliver.fa \
 670                            /ref/mm9 \
 671                            mmliver_single_without_quals
 672
 673 4) Data are the same as 1). We want to take a fragment length distribution into consideration. We set the fragment length mean to 150 and the standard deviation to 35. In addition to a BAM file, we also want to generate credibility intervals.  We allow RSEM to use 1GB of memory for CI calculation.
 674
 675  rsem-calculate-expression --bowtie-path /sw/bowtie \
 676                            --phred64-quals \
 677                            --fragment-length-mean 150.0 \
 678                            --fragment-length-sd 35.0 \
 679                            -p 8 \
 680                            --out-bam \
 681                            --calc-ci \
 682                            --ci-memory 1024 \
 683                            /data/mmliver.fq \
 684                            /ref/mm9 \
 685                            mmliver_single_quals
 686
 687 5) '/data/mmliver_paired_end_quals.bam', paired-end reads with quality scores.  We want to use 8 threads and do not generate a BAM file:
 688
 689  rsem-calculate-expression --paired-end \
 690                            --bam \
 691                            -p 8 \
 692                            /data/mmliver_paired_end_quals.bam \
 693                            /ref/mm9 \
 694                            mmliver_paired_end_quals
 695
 696 =cut
 697