查看: 530|回复: 10

[软件使用] 如何对老数据进行“翻新”?

[复制链接]

迅猛龙

Rank: 8Rank: 8

主题
257
注册时间
2020.6.16
在线时间
165 小时

发表于 7 天前 | 显示全部楼层 |阅读模式
本帖最后由 基迪奥-Jt桃 于 2021.6.8 09:10 编辑

有时候,我们从数据库下载的转录组表达量数据比较老,比如接下来用于差异分析的一个案例数据是发表于2010年的,10多年过去了,有些基因的ID 、注释信息等已经与当前数据库不匹配。

如何对老数据进行“翻新”呢?例如,本案例是人的转录组数据,我们可以利用“org.Hs.eg.db”包提供的NCBI annotation信息,只保留当前仍在使用的id,同时删除重复的基因条目。为了方便大家练习使用,本文用到的范例数据已上传到OmicShare,可直接下载使用。

范例数据下载:
https://www.omicshare.com/forum/thread-6985-1-1.html

●1.读入范例数据●

  1. rawdata <- read.delim("case01data.txt", check.names=FALSE, stringsAsFactors=FALSE)
  2. #查看前6行;
  3. head(rawdata)
复制代码


  1. #载入edgeR包;
  2. library(edgeR)
  3. #创建edgeR储存数据的DGElist对象;
  4. y <- DGEList(counts=rawdata[,4:9], genes=rawdata[,1:3])
  5. #创建的list对象结构比较简单,包括counts、samples、genes这3个表格;
  6. y
复制代码


●2.基因ID更新●

  1. #安装人的基因信息注释包“org.Hs.eg.db”;
  2. #BiocManager::install("org.Hs.eg.db")
  3. #载入人的基因信息注释包;
  4. library(org.Hs.eg.db)
  5. #提去org.Hs.eg.db包的RefSeq ID;
  6. refID <- mappedRkeys(org.Hs.egREFSEQ)
  7. head(refID)

  8. #判断数据的id是否在refID中;
  9. idfound <- y$genes$RefSeqID %in% refID
  10. #保留匹配的基因;
  11. y <- y[idfound,]
  12. dim(y)
  13. #提取注释信息表;
  14. egREFSEQ <- toTable(org.Hs.egREFSEQ)
  15. head(egREFSEQ)
复制代码


  1. #获取匹配索引,将RefSeq ID转成NCBI gene id;
  2. m <- match(y$genes$RefSeqID, egREFSEQ$accession)
  3. #将相对应的gene id 添加到DGElist对象中;
  4. y$genes$EntrezGene <- egREFSEQ$gene_id[m]

  5. #同样的方法更新基因的symbol;
  6. egSYMBOL <- toTable(org.Hs.egSYMBOL)
  7. head(egSYMBOL)
  8. m <- match(y$genes$EntrezGene, egSYMBOL$gene_id)
  9. y$genes$Symbol <- egSYMBOL$symbol[m]
  10. #查看更新后的基因信息;
  11. head(y$genes)
复制代码


Tips
  1. #match()返回第1个向量在第2个向量中的位置(索引)信息;
  2. a <- c(1,2,3,7,8)
  3. b <- c(1,2,3,4,5)
  4. c <- c(1,2,3,4,5,6,7,8)
  5. match(a,b)
  6. #[1] 1 2 3 NA NA
  7. match(a,c)
  8. #[1] 1 2 3 7 8
  9. #判断向量或数据框的元素(行)是否与前面的重复;
  10. #结果中第1个重复值为FALSE,之后为TRUE;
  11. d <- c(1,2,3,3,4,5,5,3)
  12. duplicated(d)
  13. #[1] FALSE FALSE FALSE TRUE FALSE FALSE TRUE FALSE
  14. del <- duplicated(d)
  15. #提取重复制;
  16. d[del]
  17. #去掉后面的重复值,保留第1个;
  18. d[!del]
复制代码

●3.重复数据过滤●

通过观察原始数据,我们发现一个symbol号有时会对应多个RefSeqID,比如第5行和第6行的RefSeq ID不同而Symbol号却相同。

  1. #去掉重复的RefSeqID行,只保留count值最大的行;
  2. #按照基因count之和,对数据进行降序排列;
  3. o <- order(rowSums(y$counts), decreasing=TRUE)
  4. y <- y[o,]
  5. #查看排序后的结果;
  6. head(y$counts)
  7. tail(y$counts)
复制代码


一般情况下需要过滤掉低表达的基因,但这里两个分组的总counts数都大于50,不做过滤。

  1. #提取重复值信息;
  2. d <- duplicated(y$genes$Symbol)
  3. #去掉后面的重复值,保留第1个(即count值最大的);
  4. y <- y[!d,]
  5. #查看去重后的基因数量;
  6. dim(y)
  7. #重新计算lib.size;
  8. y$samples$lib.size <- colSums(y$counts)
  9. #统一将gene id设为数据的行名;
  10. rownames(y$counts) <- rownames(y$genes) <- y$genes$EntrezGene
  11. #去掉EntrezGene列;
  12. y$genes$EntrezGene <- NULL

  13. a <- y$genes
  14. b <- y$counts
  15. cc <- data.frame(a,b)
  16. head(cc)
复制代码


  1. #导出处理后的数据;
  2. write.csv(cc,file = "case01data_filtered.csv")
复制代码

今天的内容就先到这里啦,更多实操案例欢迎查看往期文章。

基于edgeR的差异分析案例实操
如何进行差异表达分析?


本文作者:基迪奥-莫北

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
新的一天加油!
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
0
注册时间
2020.11.21
在线时间
42 小时

发表于 7 天前 | 显示全部楼层
坚持就是胜利!
回复

使用道具 举报

中华鲟

Rank: 5Rank: 5

主题
0
注册时间
2020.2.12
在线时间
51 小时

发表于 7 天前 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

中华鲟

Rank: 5Rank: 5

主题
0
注册时间
2020.2.12
在线时间
51 小时

发表于 6 天前 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
3
注册时间
2017.9.8
在线时间
40 小时

发表于 6 天前 | 显示全部楼层
加油,加油!
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
3
注册时间
2017.9.8
在线时间
40 小时

发表于 5 天前 | 显示全部楼层
加油,加油!
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
3
注册时间
2017.9.8
在线时间
40 小时

发表于 4 天前 | 显示全部楼层
加油,加油!
回复

使用道具 举报

帝王蝶

Rank: 4

主题
0
注册时间
2018.11.22
在线时间
13 小时

发表于 4 天前 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

中华鲟

Rank: 5Rank: 5

主题
0
注册时间
2020.2.12
在线时间
51 小时

发表于 3 天前 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
3
注册时间
2017.9.8
在线时间
40 小时

发表于 前天 09:12 | 显示全部楼层
加油,加油!
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
0
注册时间
2020.11.21
在线时间
42 小时

发表于 昨天 10:22 | 显示全部楼层
坚持就是胜利!
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

快速回复 返回顶部 返回列表