数据下载

力荐xena

  1. proj = "DHA"
  2. library(GEOquery)
  3. eSet <- getGEO("GSE162550", destdir = '.', getGPL = F)
  4. eSet = eSet[[1]]
  5. pd <- pData(eSet)

检查表达矩阵

  1. dat = read.delim("GSE162550_gene_sample_count_with_symbol.xls.gz")
  2. exp= as.matrix(dat[,4:9])
  3. rownames(exp) = dat$Symbol
  4. exp = exp[rownames(exp)!="---",]
  5. k = !duplicated(rownames(exp));table(k)
  6. #> k
  7. #> FALSE  TRUE 
  8. #>  1626 56142
  9. exp  = exp[k,]

基因过滤

过滤标准不唯一，一般需要过滤那些在很多样本里表达量都为0或很低的基因。

  1. nrow(exp)#过滤之前基因数
  2. #> [1] 56142
  3. #常用过滤标准1：仅去除在所有样本里表达量都为零的基因
  4. exp1 = exp[rowSums(exp)>0,]
  5. nrow(exp1)
  6. #> [1] 30863
  7. #常用过滤标准2(推荐)：仅保留在一半以上样本里表达的基因
  8. exp = exp[apply(exp, 1, function(x) sum(x > 0) >0.5*ncol(exp)), ]
  9. nrow(exp)
  10. #> [1] 22230

分组信息提取

根据样本ID的第14-15位，给样本分组（tumor和normal）

  1. library(stringr)
  2. #> Warning: package 'stringr' was built under R version 4.1.1
  3. Group = rep(c("DMSO","DHA"),each = 3)
  4. Group = factor(Group,levels = c("DMSO","DHA"))
  5. table(Group)
  6. #> Group
  7. #> DMSO  DHA 
  8. #>    3    3

保存数据

TCGA以外的数据没有clinical,surv，从下面代码里去掉。

  1. save(exp,Group,proj,pd,file = paste0(proj,".Rdata"))