研究需要,不得不断的学习新的软件包,R的突出特点即在于开源,和众多的软件包,其学习的难点也在于需要包太多,每用新的包都需要学习,了解软件的使用参数。该文记录genefu包的学使用,我自己构建的模型已经有了一定的结果,最后希望与现有的多基因芯片模型比较,其实我也不知道自己基于测序的模型能否与现有的芯片模型比较,我也完全可以到现在就整理数据,告一段落,然后继续学习其它的套路。只是因为我想深入一些,真正的做一些有意义的工作,也许可以在时间的长河里留下痕迹,提醒自己时间的力量。
####genefu包的使用
##修改镜像
options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")
if(FALSE){BiocManager::install("genefu")
}
##
library(genefu)
library(xtable)
library(rmeta)
library(Biobase)
library(caret)
##
if(FALSE){BiocManager::install("breastCancerMAINZ")
BiocManager::install("breastCancerTRANSBIG")
BiocManager::install("breastCancerUPP")
BiocManager::install("breastCancerUNT")
BiocManager::install("breastCancerNKI")
}
##加载各数据集
library(breastCancerMAINZ)
library(breastCancerTRANSBIG)
library(breastCancerUPP)
library(breastCancerUNT)
library(breastCancerNKI)
##case1比较风险预测模型
##5个芯片数据集已整合为R包
##
data(breastCancerData)
##查看数据结构
str(nki7g)
dim(nki7g)
cinfo <- colnames(pData(mainz7g))
##整合数据到向量
data.all <- c("transbig7g"=transbig7g, "unt7g"=unt7g, "upp7g"=upp7g,
"mainz7g"=mainz7g, "nki7g"=nki7g)
##mainz与NKI数据集无重复
idtoremove.all <- NULL
duplres <- NULL
##关注transbig7g数据集, unt7g upp7g
dim(pData(transbig7g))
dim(pData(unt7g))
dim(pData(upp7g))
##可观察到都为21列的临床信息,sample数不同
##rbind整合三个临床信息,
demo.all <- rbind(pData(transbig7g), pData(unt7g), pData(upp7g))
dn2 <- c("TRANSBIG", "UNT", "UPP")
## 检索VDXKIU, KIU, UPPU 三个series
ds2 <- c("VDXKIU", "KIU", "UPPU")
##is.element(A,B) A中元素属于B返回TRUE or FALSE
##有完整series信息并且留下想要的
demot <- demo.all[complete.cases(demo.all[ , c("series")]) & is.element(demo.all[ , "series"], ds2), ]
dim(demot)#可见筛选留下了290个sample
## 找到 series中相同的病人
##因demot中包含21列信息
colnames(demot)
# 找到 series id中重复的id
duplid <- sort(unique(demot[duplicated(demot[ , "id"]), "id"]))
duplrest <- NULL
##for循环中提取出游完整id dataset,且满足要求
for(i in 1:length(duplid)) {
tt <- NULL
for(k in 1:length(dn2)) {
myx <- sort(row.names(demot)[complete.cases(demot[ , c("id", "dataset")]) &
demot[ , "id"] == duplid[i] & demot[ , "dataset"] == dn2[k]])
if(length(myx) > 0) { tt <- c(tt, myx) }
}
duplrest <- c(duplrest, list(tt))
}
##只要属于我们关注的3个数据集
names(duplrest) <- duplid
duplres <- c(duplres, duplrest)
## 检索 VVDXOXFU, OXFU series
ds2 <- c("VDXOXFU", "OXFU")
demot <- demo.all[complete.cases(demo.all[ , c("series")]) & is.element(demo.all[ , "series"], ds2), ]
# Find the duplicated patients in that series
duplid <- sort(unique(demot[duplicated(demot[ , "id"]), "id"]))
duplrest <- NULL
for(i in 1:length(duplid)) {
tt <- NULL
for(k in 1:length(dn2)) {
myx <- sort(row.names(demot)[complete.cases(demot[ , c("id", "dataset")]) &
demot[ , "id"] == duplid[i] & demot[ , "dataset"] == dn2[k]])
if(length(myx) > 0) { tt <- c(tt, myx) }
}
duplrest <- c(duplrest, list(tt))
}
names(duplrest) <- duplid
duplres <- c(duplres, duplrest)
#### 得到完整的重复patients
##代码较复杂拆解查看下
a<-lapply(duplres, function(x) { return(x[-1]) })##对列表每个元素操作
unlist(a)##解开列表
duPL <- sort(unlist(lapply(duplres, function(x) { return(x[-1]) } )))
str(duPL)#duPL属于character