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;
  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 $CHAINLEN $SAMPLEGAP";
 289     if ($quiet) { $command .= " -q"; }
 290     print "$command\n";
 291     $status = system($command);
 292     if ($status != 0) {
 293         print "rsem-run-gibbs failed! Please check if you provide correct parameters/options for the pipeline!\n";
 294         exit(-1);
 295     }
 296     print "\n";
 297
 298     system("mv $sampleName.isoforms.results $imdName.isoforms.results.bak1");
 299     system("mv $sampleName.genes.results $imdName.genes.results.bak1");
 300     &collectResults("$imdName.iso_res", "$sampleName.isoforms.results"); # isoform level
 301     &collectResults("$imdName.gene_res", "$sampleName.genes.results"); # gene level
 302
 303     $command = $dir."rsem-calculate-credibility-intervals $refName $sampleName $sampleToken $CONFIDENCE $NSPC $NMB";
 304     if ($quiet) { $command .= " -q"; }
 305     print "$command\n";
 306     $status = system($command);
 307     if ($status != 0) {
 308         print "rsem-calculate-credibility-intervals failed! Please check if you provide correct parameters/options for the pipeline!\n";
 309         exit(-1);
 310     }
 311     print "\n";
 312
 313     system("mv $sampleName.isoforms.results $imdName.isoforms.results.bak2");
 314     system("mv $sampleName.genes.results $imdName.genes.results.bak2");
 315     &collectResults("$imdName.iso_res", "$sampleName.isoforms.results"); # isoform level
 316     &collectResults("$imdName.gene_res", "$sampleName.genes.results"); # gene level
 317 }
 318
 319 if (!$keep_intermediate_files) {
 320     $status = system ("rm -rf $temp_dir");
 321     if ($status != 0) {
 322         print "Fail to delete the temporary folder!\n";
 323         exit(-1);
 324     }
 325 }
 326
 327 # inpF, outF
 328 sub collectResults {
 329     my $local_status;
 330     my ($inpF, $outF);
 331     my (@results, @comment) = ();
 332     my $line;
 333     my $cnt;
 334
 335     $inpF = $_[0];
 336     $outF = $_[1];
 337
 338     $local_status = open(INPUT, $inpF);
 339     if ($local_status == 0) { print "Fail to open file $inpF!\n"; exit(-1); }
 340
 341     $cnt = 0;
 342     @results = ();
 343
 344     while ($line = <INPUT>) {
 345         ++$cnt;
 346         chomp($line);
 347         my @local_arr = split(/\t/, $line);
 348         if ($cnt == 4) { @comment = @local_arr; }
 349         else { push(@results, \@local_arr); }
 350     }
 351
 352     push(@results, \@comment);
 353     close(INPUT);
 354
 355     $local_status = open(OUTPUT, ">$outF");
 356     if ($local_status == 0) { print "Fail to create file $outF!\n"; exit(-1); }
 357
 358     my $n = scalar(@results);
 359     my $m = scalar(@{$results[0]});
 360     for (my $i = 0; $i < $m; $i++) {
 361         my @out_arr = ();
 362         for (my $j = 0; $j < $n; $j++) { push(@out_arr, $results[$j][$i]); }
 363         $" = "\t";
 364         print OUTPUT "@out_arr\n";
 365     }
 366     close(OUTPUT);
 367 }
 368
 369
 370 __END__
 371
 372 =head1 NAME
 373
 374 rsem-calculate-expression
 375
 376 =head1 SYNOPSIS
 377
 378 =over
 379
 380  rsem-calculate-expression [options] upstream_read_file(s) reference_name sample_name
 381  rsem-calculate-expression [options] --paired-end upstream_read_file(s) downstream_read_file(s) reference_name sample_name
 382  rsem-calculate-expression [options] --sam/--bam [--paired-end] input reference_name sample_name
 383
 384 =back
 385
 386 =head1 ARGUMENTS
 387
 388 =over
 389
 390 =item B<upstream_read_files(s)>
 391
 392 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.
 393
 394 =item B<downstream_read_file(s)>
 395
 396 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.
 397
 398 =item B<input>
 399
 400 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.
 401
 402 =item B<reference_name>
 403
 404 The name of the reference used.  The user must have run 'rsem-prepare-reference' with this reference_name before running this program.
 405
 406 =item B<sample_name>
 407
 408 The name of the sample analyzed. All output files are prefixed by this name (e.g., sample_name.genes.results)
 409
 410 =back
 411
 412 =head1 OPTIONS
 413
 414 =over
 415
 416 =item B<--paired-end>
 417
 418 Input reads are paired-end reads. (Default: off)
 419
 420 =item B<--no-qualities>
 421
 422 Input reads do not contain quality scores. (Default: off)
 423
 424 =item B<--strand-specific>
 425
 426 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)
 427
 428 =item B<--sam>
 429
 430 Input file is in SAM format. (Default: off)
 431
 432 =item B<--bam>
 433
 434 Input file is in BAM format. (Default: off)
 435
 436 =item B<--sam-header-info> <file>
 437
 438 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: "")
 439
 440 =item B<-p/--num-threads> <int>
 441
 442 Number of threads to use. Both Bowtie and expression estimation will use this many threads. (Default: 1)
 443
 444 =item B<--out-bam>
 445
 446 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)
 447
 448 =item B<--calc-ci>
 449
 450 Calculate 95% credibility intervals and posterior mean estimates.  (Default: off)
 451
 452 =item B<--seed-length> <int>
 453
 454 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)
 455
 456 =item B<--tag> <string>
 457
 458 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: "")
 459
 460 =item B<--bowtie-path> <path>
 461
 462 The path to the bowtie executables. (Default: the path to the bowtie executables is assumed to be in the user's PATH environment variable)
 463
 464 =item B<--bowtie-n> <int>
 465
 466 (Bowtie parameter) max # of mismatches in the seed. (Range: 0-3, Default: 2)
 467
 468 =item B<--bowtie-e> <int>
 469
 470 (Bowtie parameter) max sum of mismatch quality scores across the alignment. (Default: 99999999)
 471
 472 =item B<--bowtie-m> <int>
 473
 474 (Bowtie parameter) suppress all alignments for a read if > <int> valid alignments exist. (Default: 200)
 475
 476 =item B<--phred33-quals>
 477
 478 Input quality scores are encoded as Phred+33. (Default: on)
 479
 480 =item B<--phred64-quals>
 481
 482 Input quality scores are encoded as Phred+64 (default for GA Pipeline ver. >= 1.3). (Default: off)
 483
 484 =item B<--solexa-quals>
 485
 486 Input quality scores are solexa encoded (from GA Pipeline ver. < 1.3). (Default: off)
 487
 488 =item B<--forward-prob> <double>
 489
 490 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)
 491
 492 =item B<--fragment-length-min> <int>
 493
 494 Minimum read/insert length allowed. This is also the value for the bowtie -I option. (Default: 1)
 495
 496 =item B<--fragment-length-max> <int>
 497
 498 Maximum read/insert length allowed. This is also the value for the bowtie -X option. (Default: 1000)
 499
 500 =item B<--fragment-length-mean> <double>
 501
 502 (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)
 503
 504 =item B<--fragment-length-sd> <double>
 505
 506 (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>)
 507
 508 =item B<--estimate-rspd>
 509
 510 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)
 511
 512 =item B<--num-rspd-bins> <int>
 513
 514 Number of bins in the RSPD. Only relevant when '--estimate-rspd' is specified.  Use of the default setting is recommended. (Default: 20)
 515
 516 =item B<--ci-memory> <int>
 517
 518 Amount of memory (in MB) RSEM is allowed to use for computing credibility intervals. (Default: 1024)
 519
 520 =item B<--keep-intermediate-files>
 521
 522 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)
 523
 524 =item B<-q/--quiet>
 525
 526 Suppress the output of logging information. (Default: off)
 527
 528 =item B<-h/--help>
 529
 530 Show help information.
 531
 532 =back
 533
 534 =head1 DESCRIPTION
 535
 536 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.
 537
 538 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:
 539
 540   sort -k 1,1 -s input.sam > input.sorted.sam
 541
 542 The SAM/BAM format RSEM uses is v1.3. However, it is compatible with old SAM/BAM format.
 543
 544 The user must run 'rsem-prepare-reference' with the appropriate reference before using this program.
 545
 546 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.
 547
 548 Please note that some of the default values for the Bowtie parameters are not the same as those defined for Bowtie itself.
 549
 550 The temporary directory and all intermediate files will be removed when RSEM finishes unless '--keep-intermediate-files' is specified.
 551
 552 With the "--calc-ci" option, 95% credibility intervals and posterior mean estimates will be calculated in addition to maximum likelihood estimates.
 553
 554 =head1 OUTPUT
 555
 556 =over
 557
 558 =item B<sample_name.genes.results>
 559
 560 File containing gene level expression estimates. The format of each
 561 line in this file is:
 562
 563 gene_id expected_counts tau_value [pmc_value tau_pme_value tau_ci_lower_bound tau_ci_upper_bound] transcript_id_list
 564
 565 Fields are separated by the tab character. Fields within "[]" are only
 566 presented if '--calc-ci' is set. pme stands for posterior mean
 567 estimation. pmc stands for posterior mean counts. ci_lower_bound(l)
 568 means the lower bound of the credibility intervals, ci_upper_bound(u)
 569 means the upper bound of the credibility intervals. So the credibility
 570 interval is [l, u]. 'transcript_id_list' is a space-separated list of
 571 transcript_ids belonging to the gene.
 572
 573 =item B<sample_name.isoforms.results>
 574
 575 File containing isoform level expression values. The format of each
 576 line in this file is:
 577
 578 transcript_id expected_counts tau_value [pmc_value tau_pme_value tau_ci_lower_bound tau_ci_upper_bound] other_attributes
 579
 580 Fields are separated by the tab character. 'other_attributes' are all
 581 other attributes after attribute 'transcript_id' field in the GTF
 582 file. If no other attributes are given or no GTF file is provided in
 583 'rsem-prepare-reference', there will be no tab after the
 584 tau_value field.
 585
 586 =item B<sample_name.bam, sample_name.sorted.bam and sample_name.sorted.bam.bai>
 587
 588 Only generated when --out-bam is specified.
 589
 590 'sample_name.bam' is a BAM-formatted file of read alignments in
 591 genomic coordinates. Alignments of reads that have identical genomic
 592 coordinates (i.e., alignments to different isoforms that share the
 593 same genomic region) are collapsed into one alignment.  The MAPQ field
 594 of each alignment is set to max(100, floor(-10 * log10(1.0 - w) +
 595 0.5)), where w is the posterior probability of that alignment being
 596 the true mapping of a read.  In addition, RSEM pads a new tag
 597 ZW:f:value, where value is a single precision floating number
 598 representing the posterior probability.
 599
 600 'sample_name.sorted.bam' and 'sample_name.sorted.bam.bai' are the
 601 sorted BAM file and indices generated by samtools (included in RSEM package).
 602
 603 =item B<sample_name.stat>
 604
 605 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.
 606
 607 =back
 608
 609 =head1 EXAMPLES
 610
 611 Assume the path to the bowtie executables is in the user's PATH environment variable. Reference files are under '/ref' with name 'mm9'.
 612
 613 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:
 614
 615  rsem-calculate-expression --phred64-quals \
 616                            -p 8 \
 617                            --out-bam \
 618                            /data/mmliver.fq \
 619                            /ref/mm9 \
 620                            mmliver_single_quals
 621
 622 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:
 623
 624  rsem-calculate-expression -p 8 \
 625                            --paired-end \
 626                            /data/mmliver_1.fq \
 627                            /data/mmliver_2.fq \
 628                            /ref/mm9 \
 629                            mmliver_paired_end_quals
 630
 631 3) '/data/mmliver.fa', single-end reads without quality scores. We want to use 8 threads and generate a BAM file:
 632
 633  rsem-calculate-expression -p 8 \
 634                            --no-qualities \
 635                            /data/mmliver.fa \
 636                            /ref/mm9 \
 637                            mmliver_single_without_quals
 638
 639 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.
 640
 641  rsem-calculate-expression --bowtie-path /sw/bowtie \
 642                            --phred64-quals \
 643                            --fragment-length-mean 150.0 \
 644                            --fragment-length-sd 35.0 \
 645                            -p 8 \
 646                            --out-bam \
 647                            --calc-ci \
 648                            --ci-memory 1024 \
 649                            /data/mmliver.fq \
 650                            /ref/mm9 \
 651                            mmliver_single_quals
 652
 653 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:
 654
 655  rsem-calculate-expression --paired-end \
 656                            --bam \
 657                            -p 8 \
 658                            /data/mmliver_paired_end_quals.bam \
 659                            /ref/mm9 \
 660                            mmliver_paired_end_quals
 661
 662 =cut