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