When I run RSEM and then join together the different .isoforms.results files with per sample expression counts into a raw counts matrix (using Trinity script abundance_estimates_to_matrix.pl), this matrix is unsorted and does not carry over the same sorting as the isoforms.results files.
For instance, a raw counts matrix I have starts with the following:
control1 control2 control3 control4 exp1 exp2 exp3 exp4
kf_t00010653-RA 28.00 31.00 34.00 33.00 92.00 100.00 75.00 97.00
kf_t00014892-RA 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
kf_t00015024-RA 14.45 6.00 6.00 0.00 4.00 2.00 0.00 1.00
The RSEM *isoforms.results file has this sorting:
transcript_id gene_id length effective_length expected_count TPM FPKM IsoPct
kf_t00000001-RA kf_g00000001 5182 4933.94 191.91 5.86 5.34 100.00
kf_t00000002-RA kf_g00000002 2390 2141.94 125.88 8.85 8.07 100.00
kf_t00000003-RA kf_g00000003 654 405.97 0.00 0.00 0.00 0.00
This means that if you casually import the gene lengths into the DGEList object in edgeR like this (length is a table generated from the above Trinity script containing lengths in the column length):
y <- DGEList(counts=data.set, group=treatment)
y$genes = data.frame(length$length)
Does this mean that the lengths for kf_t00000001-RA will be used for kf_t00010653-RA, leading to nonsense results later on when rpkm() from EdgeR is used with transcript lengths?
After all, y$genes is just a list of lengths with no labels attached to them with the first length in the vector being "5182" but the first transcript in the y$counts is "kf_t00010653-RA".
I have tried sorting the raw counts matrix before putting it into a DGEList object like this:
data.set <- data.set[order(row.names(data.set)), ]
so that the y$counts matrix has the same order as the file used to extract gene lists. Then the lengths should match the correct transcript.
Does this look like a viable solution?
