text: user3541: I am running the STAR-RSEM pipeline, I have a simple shell script(*), which I use to loop over fasta files corresponding to different samples. From the output, I am interested in the *.genes.results files, for which the columns are:
  > names(my_df1)
  [1] "gene_id" "transcript_id(s)" "length" "effective_length" "expected_count" "TPM" "FPKM"            
  
  To my surprise, whereas the gene_id and transcript_id(s) columns are the same, I see different length and effective_length values for the same genes in different sample results (omitting the transcript_id(s) column for formatting purposes):
  > my_df1 <- fread(paste0("alignment/", my_files[1]))
  > head(my_df1[1:5,c(1,3:6)])
                gene_id  length effective_length expected_count   TPM
  1: ENSG00000000003.14 2292.23          2242.23            978 26.69
  2:  ENSG00000000005.5 1339.00          1289.00             11  0.52
  3: ENSG00000000419.12 1024.45           974.45            469 29.45
  4: ENSG00000000457.13 4204.00          4154.00            241  3.55
  5: ENSG00000000460.16 2953.40          2903.40             62  1.31
  > 
  > my_df2 <- fread(paste0("alignment/", my_files[2]))
  > head(my_df2[1:5,c(1,3:6)])
                gene_id  length effective_length expected_count    TPM
  1: ENSG00000000003.14 2260.38          2210.38       10052.00 117.78
  2:  ENSG00000000005.5 1269.22          1219.22         132.00   2.80
  3: ENSG00000000419.12 1058.13          1008.13        4386.00 112.68
  4: ENSG00000000457.13 4175.05          4125.05         873.00   5.48
  5: ENSG00000000460.16 2783.92          2733.92         780.01   7.39
  > 
  > my_df3 <- fread(paste0("alignment/", my_files[3]))
  > head(my_df3[1:5,c(1,3:6)])
                gene_id  length effective_length expected_count  TPM
  1: ENSG00000000003.14 2662.59          2612.59             49 2.90
  2:  ENSG00000000005.5  940.50           890.50              0 0.00
  3: ENSG00000000419.12  956.87           906.87             58 9.89
  4: ENSG00000000457.13 6364.00          6314.00             66 1.62
  5: ENSG00000000460.16 2661.00          2611.00             15 0.89
  
  The length and effective_length are defined per gene (and its transcripts) and should not differ across samples. Below is from the manual:
  
  A gene's 'length' and 'effective_length' are defined as the weighted
    average of its transcripts' lengths and effective lengths (weighted by
    'IsoPct')
  
  I would appreciate any hints pointing out what I am missing.
  * The script I am using:
  for file in ...; \
  do rsem-calculate-expression \
  -p 8 \
  --star \
  # --output-genome-bam \
  --star-path path-to-star-executable \
  --star-gzipped-read-file \
  --forward-prob=0 \
  $file \
  path-to-genome \ # should end with the string used for the "reference_name" parameter of "rsem-prepare-reference"
  ${file:11:7}; \ # sample prefix
  done
  
  EDIT: For those who will be interested, this entry from the RSEM Users Google Group demonstrates how to calculate gene length from weighted transcript length.
  
  user77: If there are differences in the expression of individual transcripts between samples then both the length and effective length will differ between samples. It is commonly the case that the ratio of longer and short transcript expression varies a bit between samples, which gives rise to this.
  
original_text: 3541: <p>I am running the <code>STAR-RSEM</code> pipeline, I have a simple shell script(*), which I use to loop over fasta files corresponding to different samples. From the output, I am interested in the <code>*.genes.results</code> files, for which the columns are:</p>
  
  <pre><code>&gt; names(my_df1)
  [1] "gene_id" "transcript_id(s)" "length" "effective_length" "expected_count" "TPM" "FPKM"            
  </code></pre>
  
  <p>To my surprise, whereas the <code>gene_id</code> and <code>transcript_id(s)</code> columns are the same, I see different <code>length</code> and <code>effective_length</code> values for the same genes in different sample results (omitting the <code>transcript_id(s)</code> column for formatting purposes):</p>
  
  <pre><code>&gt; my_df1 &lt;- fread(paste0("alignment/", my_files[1]))
  &gt; head(my_df1[1:5,c(1,3:6)])
                gene_id  length effective_length expected_count   TPM
  1: ENSG00000000003.14 2292.23          2242.23            978 26.69
  2:  ENSG00000000005.5 1339.00          1289.00             11  0.52
  3: ENSG00000000419.12 1024.45           974.45            469 29.45
  4: ENSG00000000457.13 4204.00          4154.00            241  3.55
  5: ENSG00000000460.16 2953.40          2903.40             62  1.31
  &gt; 
  &gt; my_df2 &lt;- fread(paste0("alignment/", my_files[2]))
  &gt; head(my_df2[1:5,c(1,3:6)])
                gene_id  length effective_length expected_count    TPM
  1: ENSG00000000003.14 2260.38          2210.38       10052.00 117.78
  2:  ENSG00000000005.5 1269.22          1219.22         132.00   2.80
  3: ENSG00000000419.12 1058.13          1008.13        4386.00 112.68
  4: ENSG00000000457.13 4175.05          4125.05         873.00   5.48
  5: ENSG00000000460.16 2783.92          2733.92         780.01   7.39
  &gt; 
  &gt; my_df3 &lt;- fread(paste0("alignment/", my_files[3]))
  &gt; head(my_df3[1:5,c(1,3:6)])
                gene_id  length effective_length expected_count  TPM
  1: ENSG00000000003.14 2662.59          2612.59             49 2.90
  2:  ENSG00000000005.5  940.50           890.50              0 0.00
  3: ENSG00000000419.12  956.87           906.87             58 9.89
  4: ENSG00000000457.13 6364.00          6314.00             66 1.62
  5: ENSG00000000460.16 2661.00          2611.00             15 0.89
  </code></pre>
  
  <p>The <code>length</code> and <code>effective_length</code> are defined per gene (and its transcripts) and should not differ across samples. Below is from the <a href="http://deweylab.biostat.wisc.edu/rsem/rsem-calculate-expression.html" rel="nofollow noreferrer">manual</a>:</p>
  
  <blockquote>
    <p>A gene's 'length' and 'effective_length' are defined as the weighted
    average of its transcripts' lengths and effective lengths (weighted by
    'IsoPct')</p>
  </blockquote>
  
  <p>I would appreciate any hints pointing out what I am missing.</p>
  
  <p>* The script I am using:</p>
  
  <pre><code>for file in ...; \
  do rsem-calculate-expression \
  -p 8 \
  --star \
  # --output-genome-bam \
  --star-path path-to-star-executable \
  --star-gzipped-read-file \
  --forward-prob=0 \
  <span class="math-container">$file \
  path-to-genome \ # should end with the string used for the "reference_name" parameter of "rsem-prepare-reference"
  $</span>{file:11:7}; \ # sample prefix
  done
  </code></pre>
  
  <p>EDIT: For those who will be interested, <a href="https://groups.google.com/d/msg/rsem-users/gDh0OBLK6_4/z-3ZoETxn5sJ" rel="nofollow noreferrer">this entry</a> from the RSEM Users Google Group demonstrates how to calculate gene length from weighted transcript length.</p>
   77: <p>If there are differences in the expression of individual transcripts between samples then both the length and effective length will differ between samples. It is commonly the case that the ratio of longer and short transcript expression varies a bit between samples, which gives rise to this.</p>
  