1
$\begingroup$

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.

$\endgroup$
2
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.