# Select top 100 genes ranked by variance in read counts

I have a gene expression data (read counts) table with many samples. Showing some of them below:

Ensembl_ID          A   B   C   D   E   F   G
ENSG00000236601.2   0   0   0   1   0   0   1
ENSG00000237094.12  2   187 442 109 88  144 486
ENSG00000269732.1   0   0   0   0   0   0   0
ENSG00000284733.1   0   0   0   0   0   0   0
ENSG00000233653.3   0   0   0   0   0   0   0
ENSG00000250575.1   1   11  17  5   5   8   17
ENSG00000278757.1   0   0   2   0   1   0   7
ENSG00000230021.9   0   25  30  27  9   8   11
ENSG00000235146.2   0   0   0   0   0   0   0
ENSG00000225972.1   0   7   3016    11  5   14  5
ENSG00000225630.1   1   113 194 76  91  47  94
ENSG00000237973.1   25  1037    9767    1160    886 1321    5220


I want to select top 100 genes with reverse sort of variance. How to do that?

I assume by "reverse sort of variance" you mean "highest variance". Assuming you made a matrix out of that (set the row names to the first column and then remove it) and called it m:

sel = order(apply(m, 1, var), decreasing=TRUE)[1:100]


sel then contains the indices into your matrix (or the original dataframe).

BTW, I hope you're not filtering by variance before performing differential expression. That would be incredibly incorrect. Doing this sort of filtering for performing PCA or something like that can be useful, though you should really be using normalized counts (or a transformation of them).

• Thanks. yes, I'm doing this filtering for PCA not for differential expression. Apr 24 '18 at 8:44
• BTW, before doing this I converted Ensembl_ID to gene names and I have multiple genes now so I did like this U2 <- aggregate(. ~ gene_name, data = U, max) and on that I will be calculating variance. Do you think this is right? Apr 24 '18 at 8:47
• Usually what I do is using the mean of the different transcripts for each gene, but it is right (if you can justify it).
– llrs
Apr 24 '18 at 8:53
• @raju Do not convert to gene names until you are completely finished with the analysis. Apr 24 '18 at 8:56
• Ok. But if there are duplicate Ensembl_ids can I do this U2 <- aggregate(. ~ Ensembl_ID, data = U, max) Apr 24 '18 at 9:03