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