GEO数据挖掘用dat[ids$probe_id,] 取子集都是

GSE44077

数据集是GSE44077，pd信息经过取冗余处理后得到了226行，11列的数据。

 pd <-  pd[,apply(pd, 2, function(x){
    length(unique(x))>1
  })]

第一个疑问，如何从pd$characteristics_ch1这一列中提取NSCLC和tumor都含有的子集呢

> table(grepl(c('NSCLC','tumor'),pd$characteristics_ch1))

FALSE  TRUE 
  171    55 
Warning message:
In grepl(c("NSCLC", "tumor"), pd$characteristics_ch1) :
  变数'pattern'的长度大于一，因此只能用其第一元素

得到的pd中的tissue这一列如下图

<img src="https://tva1.sinaimg.cn/large/006y8mN6gy1g95pvms85cj305s0p0jsu.jpg" alt="image-20191119205419711" style="zoom:50%;" />

上图中的应该是NSCLC和tumor是一组的，我需要选择这二者所在的行，一起提取出来，用了上面的的table(grepl(c('NSCLC','tumor'),pd$characteristics_ch1))的代码后，出现了报错，然后我又尝试如下，

pd_tmp <- pd
x <-  c(grepl('NSCLC',pd$`tissue:ch1`) + grepl('tumor',pd$`tissue:ch1`))
pd_tmp <- pd[x,]
# table(x)

不过我觉得一定是可以优雅地提出子集来，但是我不会！不过不能纠结，就放一小下，下次学了那个函数，再做一下！

第二个疑问，用名字提取自己dat_tmp[ids$probe_id,]时，出现NA

image-20191119200341934

解决如下，发现了一个很好的坑，哈哈哈，我遇到过一次，当dat的行名也就是探针名（ID/probe_id）都是数字的时候，这个时候用dat_tmp=dat_tmp[ids$probe_id,]就会全部变成NA，而用match函数解决了

if(T){
  ids <- probe2gene
  dat <- as.data.frame(dat)
  ids$median=apply(dat,1,median) #ids新建median这一列，列名为median，同时对dat这个矩阵按行操作，取每一行的中位数，将结果给到median这一列的每一行
  ids=ids[order(ids$symbol,ids$median,decreasing = T),]#对ids$symbol按照ids$median中位数从大到小排列的顺序排序，将对应的行赋值为一个新的ids
  ids=ids[!duplicated(ids$symbol),]#将symbol这一列取取出重复项，'!'为否，即取出不重复的项，去除重复的gene ，保留每个基因最大表达量结果s
  head(ids)
  dat[1:4,1:4]
  dat_tmp <- dat
  dat_tmp[1:4,1:4]
  dat_tmp=dat_tmp[ids$probe_id,] #新的ids取出probe_id这一列，将dat按照取出的这一列中的每一行组成一个新的dat
  dat_tmp[1:4,1:4]#踩坑啊
  #重新复制dat_tmp
  dat_tmp <- dat_tmp[match(ids$probe_id,rownames(dat_tmp)),]
  dat_tmp[1:4,1:4]
  identical(ids$probe_id,rownames(dat_tmp))
  class(ids$probe_id)
  class(rownames(dat))
  identical(as.character(ids$probe_id),rownames(dat_tmp))
  dat <- dat_tmp
  rownames(dat)=ids$symbol#把ids的symbol这一列中的每一行给dat作为dat的行名
  dat[1:4,1:4]  #保留每个基因ID第一次出现的信息
  dim(dat)
  boxplot(dat)
  dat=log(dat+1)
  dat[1:4,1:4] 
  
  save(dat,ids,group_list,file = 'step1-getdata.Rdata')
  load(file = 'step1-getdata.Rdata')
  
}

最后友情宣传生信技能树

• 生物信息学“义诊”

• 生物信息学"拍卖会"

• 全国巡讲:R基础，Linux基础和RNA-seq实战演练 : 预告：12月28-30长沙站

• 广州珠江新城GEO数据挖掘滚动开班

• DNA及RNA甲基化数据分析与课题设计