有的GEO数据分组信息找起来比较麻烦
2019-02-18 本文已影响50人
小洁忘了怎么分身
getGEO下载了表达矩阵文件后可获得他的ExpressionSet对象,用pdata函数提取临床信息表格,命名为pd。此处需要耗时端详这个pd表格,找到里面的分组信息列,生成一个grouplist,后续分析会多次用到。
有的直接拿出其中一列即可,有时候则不太好找,这时就考验R语言技巧了。有的pd表格有几十列,先减小一下筛选范围,我写了一个函数,用来汇总表格,不要关心function里面的内容,就是直接用。
dumd <- function(df){
library(dplyr)
colname <- vector("character")
count <- vector("integer")
for(i in 1:ncol(df)){
colname[i] = colnames(df)[[i]]
count[i]=nrow(df[!duplicated(df[,i]),])
}
output <- tibble(colname,count)
print(output)
}
x=dumd(pd)
dumd函数其实就是统计了一下每列有多少非重复值。只有一个的,说明所有行都一样,直接去掉即可。重点就是,所有行对应值都一样的列,是不会有分组信息在里面的,所以直接删掉这部分。如果有>1的个位数的,很有可能就是分组信息咯。
pd=pd[,x$count>1]
还没那么容易,如果没有直接可用的列,grouplist大概可以有这三种生成方式:
#第一类,简单的。
pd$title #tab试一下
nchar(pd$characteristics_ch1[1])
group_list=substring(pd$characteristics_ch1,5,5)
which(colnames(pd)=="characteristics_ch1")
#第二类,复杂的
group_list=ifelse(str_detect(pd$title,"TT2")==TRUE,"TT2","TT3")
group_list=ifelse(str_detect(pd$title,"Control")==TRUE,"contorl","treat")
pd$characteristics_ch1.1
group_list=ifelse(pd$characteristics_ch1.1=="triple-negative status: not TN","notTN","TN")
#第三类,更复杂的,如果分组有三组以上 参考这个代码(来自生信技能树)
g=ifelse(deg$P.Value>0.01,'stable',
ifelse( deg$logFC >logFC_t,'up',
ifelse( deg$logFC < -logFC_t,'down','stable') )
)
ifelse是一个由if循环进化来的函数,原本语法是:
if(boolean_expression) {
// statement(s) will execute if the boolean expression is true.
} else {
// statement(s) will execute if the boolean expression is false.
}
简写成了ifelse(test, yes, no),分情况将数据分组。