1
Looks okay to me. You should always check the output of each command so that you are sure that you know what it is doing. Here is a fully reproducible example:
[1] Download study data from GEO that used AFFY HG U95Av2
library(Biobase)
library(GEOquery)
human <- getGEO("GSE713", GSEMatrix =TRUE, getGPL=FALSE)
if (length(human) > 1) idx <- grep("GPL8300", attr(human, "names")) else idx <- 1
human <- human[[idx]]
human is an ExpressionSet object with Affymetrix probe IDs as rownames
human
ExpressionSet (storageMode: lockedEnvironment)
assayData: 12623 features, 20 samples
element names: exprs
protocolData: none
phenoData
sampleNames: 18.3 24.2 ... 23.2 (20 total)
varLabels: title geo_accession ... data_row_count (30 total)
varMetadata: labelDescription
featureData: none
experimentData: use 'experimentData(object)'
pubMedIds: 12907719
Annotation: GPL8300
rownames(human)[1:8]
[1] "1002_f_at" "1003_s_at" "1004_at" "1005_at" "1006_at" "1007_s_at"
[7] "1008_f_at" "1009_at"
[2] Look up the probe IDs with biomaRt
library("biomaRt")
mart <- useMart("ENSEMBL_MART_ENSEMBL")
mart <- useDataset("hsapiens_gene_ensembl", mart)
annotLookup <- getBM(
mart=mart,
attributes=c("affy_hg_u95av2","hgnc_symbol"),
filter = "affy_hg_u95av2",
values = rownames(human),
uniqueRows=TRUE)
biomaRt never returns data in the same order as it's input:
head(annotLookup)
affy_hg_u95av2 hgnc_symbol
1 1787_at CDKN1C
2 157_at PRAME
3 1227_g_at ADAM17
4 1226_at ADAM17
5 181_g_at MBOAT7
6 180_at MBOAT7
So, we have to align annotLookup to our probe ID rownames.
[3] Align annotLookup to our probe ID rownames
indicesLookup <- match(rownames(human), annotLookup$affy_hg_u95av2)
Does order match?
head(data.frame(rownames(human), annotLookup[indicesLookup,]))
rownames.human. affy_hg_u95av2 hgnc_symbol
178 1002_f_at 1002_f_at
247 1003_s_at 1003_s_at CXCR5
248 1004_at 1004_at CXCR5
184 1005_at 1005_at DUSP1
659 1006_at 1006_at MMP10
table(rownames(human) == annotLookup[indicesLookup,'affy_hg_u95av2'])
TRUE
11611
Great. However, indicesLookup will contain a NA value in situations where a probe ID was not even found via biomaRt.
You'll also notice that some probe IDs have no HGNC symbols, as indicated by empty values in the hgnc_symbol column. For a data matrix object like that used within ExpressionSets, empty rownames / non-unique rownames are tolerated, but they are not tolerated in data-frames.
For these empty values, we can use the original probe ID, if we want, and create a new column in annotLookup:
annotLookuphgnc_symbol == "",
annotLookuphgnc_symbol)
head(annotLookup)
affy_hg_u95av2 hgnc_symbol hgnc_X_affy
1 1787_at CDKN1C CDKN1C
2 157_at PRAME PRAME
3 1227_g_at ADAM17 ADAM17
4 1226_at ADAM17 ADAM17
5 181_g_at MBOAT7 MBOAT7
6 180_at MBOAT7 MBOAT7
tail(annotLookup)
affy_hg_u95av2 hgnc_symbol hgnc_X_affy
13530 736_f_at 736_f_at
13531 896_at 896_at
13532 733_at 733_at
13533 854_at 854_at
13534 41782_g_at 41782_g_at
13535 944_s_at 944_s_at
Finally, to ensure complete compliance with data-frames, we can ensure that each gene has a unique identifier (optional):
[4] Ensure that each gene has a unique identifier to avoid duplicate issues
newnames <- paste(
annotLookup[indicesLookup, "hgnc_X_affy"],
c(1:length(indicesLookup)),
sep="_")
rownames(human) <- newnames
exprs(human)[1:8,1:5]
18.3 24.2 36.9 24.7 47.5
1002_f_at_1 6.8 4.1 2.8 5.5 5.8
CXCR5_2 34.3 27.8 25.0 11.0 31.8
CXCR5_3 6.3 7.4 37.9 11.7 36.0
DUSP1_4 636.4 565.0 418.1 454.4 180.1
MMP10_5 8.5 5.0 2.1 15.6 1.9
DDR1_6 142.5 181.4 231.6 222.1 185.7
EIF2AK2_7 809.6 892.2 1243.7 774.7 1041.8
HINT1_8 807.9 674.6 990.0 670.3 822.0
tail(exprs(human))[,1:5]
18.3 24.2 36.9 24.7 47.5
NA_12618 4.3 2.7 3.0 14.3 23.3
NA_12619 1.6 2.2 1.8 1.7 5.7
NA_12620 4.9 1.3 17.8 1.3 15.5
NA_12621 1.3 1.4 11.1 1.9 3.8
