美文网首页
2022-05-09

2022-05-09

作者: 玄武Ken | 来源:发表于2022-05-09 00:11 被阅读0次

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_X_affy <- ifelse( annotLookuphgnc_symbol == "",
annotLookupaffy_hg_u95av2, 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

相关文章

  • 橙子的ScalersTalk第六轮新概念朗读持续力训练Day 2

    练习材料:[Day 2783 2022-05-09] L51-3:Predicting the future Ba...

  • 养家夜读重点

    尾盘降温 原创 不易 养家不易 2022-05-09 20:49 情绪有周期 总是周而复始的循环 低潮-升温-发酵...

  • 礼多人不怪-学习85种感谢的说法

    2022-05-09礼多人不怪。感谢的话多说没错。还得说好。https://www.berries.com/blo...

  • 2022-05-09

    Cell Stem Cell | 培养心肌细胞获突破,可用于替代衰老心肌 图灵基因图灵基因2022-05-09 1...

  • 0189教练|感受你的感受

    2022-05-09 北京 晴天☀️ 今天决定要重新回归跑步,今天跑鞋也到了。我来继续挑战徒步108km,每天早上...

  • 审稿怎能断章取义?

    2022-05-09 多云 星期一 文字现在都那么功利了吗? 如果写文章是直奔市场而去的,这样的文章已经不是...

  • 2022-05-09

    2022-05-09坚持分享第1317天 读《知行合一2》P127-135,感悟:阳明心学命运——分裂。左派心学洪...

  • 和宝贝聊作文

    2022-05-09 阴 周一 宝贝的期中考试上周就考完了。 我本不想询问她考后的感觉,但是上...

  • 督导收获

    中原焦点团队 坚持分享1395天 2022-05-09 今晚,是刘老师对我们高七和高六的联合督导。通过刘老师的...

  • 饺面店的第八十二天:申请外卖

    2022-05-09 一天又将结束,拖着疲惫的身体,突然感觉脑袋里空空的:这一天,我做了什么?稀里糊涂、忙忙碌...

网友评论

      本文标题:2022-05-09

      本文链接:https://www.haomeiwen.com/subject/pgvcurtx.html