10

I see that TCGA RNASeq V2 RSEM data is normalized with upper-quartile normalization.

After doing Quantification with RSEM with the samples I have, I got "genes.results" as output which has gene id, transcript id(s), length, expected count, and FPKM. So, from all the sample output files I got the gene_id and expected_counts [from all samples]. for eg: it looks like following:

Gene_id S1  S2  S3  S4  S5
ENSG00000000003 1355    1121    242 709 1071
ENSG00000000005 5   0   0   0   0
ENSG00000000419 839 1345    48  976 1234
ENSG00000000457 429 1803.08 386 1156    628
ENSG00000000460 404 1444    523 479 1260
ENSG00000000938 294 312 93  189 683
ENSG00000000971 3911    4633.92 264 2863    5194
ENSG00000001036 1546    2276    271 1641    2141
ENSG00000001084 1945    2125    490 980 1533
ENSG00000001167 2054    4069.72 3457    2075    2207
ENSG00000001460 153 339 77  221 273

I want to apply "upper-quartile" normalization on this data. But doesn't know which R package I can use and which function?

Is anyone aware about this?

P.S. This is not for differential analysis

llrs
  • 4,693
  • 1
  • 18
  • 42
stack_learner
  • 1,262
  • 14
  • 26

1 Answers1

9

You can use the quantile function in base R to get the value of a particular quantile (e.g. 0.75 for the upper quartile). This can then be used as a factor for normalisation: divide the observed expression by this number.

## set up simulated data
data.mat <- floor(10^data.frame(sapply(1:5,function(x){rnorm(100)})));
colnames(data.mat) <- paste0("S",1:5);
data.mat[data.mat > 0] <- floor(data.mat[data.mat > 0]);

It's generally a good idea to first remove genes that have no counts in any samples. Rows can be removed using a logical selection based on the total expression in any given row.

## Only retain rows that have a positive total expression per gene
data.mat <- data.mat[rowSums(data.mat) > 0,];

Two methods of upper quartile determination are presented here, one which is based on the total counts, and one based on the expressed counts (i.e. excluding values with no expression). These both use the quantile function of R. The apply function calculates statistics that are based on a single specified dimension (or multiple dimensions) of the data. The first dimension of a 2D matrix in R corresponds to the row, while the second corresponds to the column. By using 2 in the apply function, data are collected along the second dimension (i.e. the column), and statistics generated from that:

data.quantileAll <- apply(data.mat, 2, function(x){quantile(x, 0.75)});
data.quantileExpressed <-
  apply(data.mat, 2, function(x){quantile(x[x>0], 0.75)});

Expressed counts may be preferred when read counts are low, so it's difficult to distinguish between a transcript with very low expression and a transcript with no expression.

Normalisation is applied by dividing the values by the column-based upper quartile statistic. R will repeat values along the last dimension first (which in this case is the column), so the matrix needs to be transposed first (with columns as genes) to get the arithmetic working properly, with a final transposition to return the orientation back to the original direction (with columns as samples).

data.norm <- t(t(data.mat) / data.quantileExpressed);
gringer
  • 14,012
  • 5
  • 23
  • 79