TCGA分析-数据下载2

### 小何开始运行

#1.数据下载 #从Gene Expression Omnibus (GEO)数据库下载数据
library(GEOquery)
proj = "GSE106899"
eSet <- getGEO(proj, destdir = '.', getGPL = F)

## Found 1 file(s)

## GSE106899_series_matrix.txt.gz

## Using locally cached version: ./GSE106899_series_matrix.txt.gz

class(eSet)

## [1] "list"

length(eSet)

## [1] 1

#eSet = eSet[[1]]
#在R语言中，eSet = eSet[[1]] 这句代码是用来提取 eSet 数据框中的第一列数据。
#eSet 通常是一个包含多个数据集的对象，这些数据集可能来自一个生物实验。在这些数据集中，第一列数据可能是样本的标识符、组别、条件、处理方式等表型数据。通过将 eSet 数据框中的第一列赋值给新的变量 eSet，可以方便地对这些数据进行后续的分析和处理。
#上述代码提取表达矩阵，但是提取出来是0行，不存在。
#2.提取表达矩阵
#clinical<- pData(eSet)
#具体来说，pData()函数是从eSet中提取“数据”部分,即提取临床信息。这个函数通常与setNames()函数一起使用，后者为数据框的列设置名称。
#phenoData的全称是表型数据。在生物信息学中，它通常指的是描述样本信息的临床数据，如年龄、性别、治疗手段等。
#.提取表达矩阵 read.delim函数用于读取以制表符为分隔符的文本文件，并将其解析为数据框（data frame）对象。它通常用于读取以  .txt  或  .tsv  格式保存的数据文件。 row.names  参数设置为 1，您可以指定数据框中的第一列作为行名。
#隐式循环+自定义函数
#第一种方法
fs=dir("GSE106899_RAW/")
re = list() 
for (i in 1:length(fs)){
  re[[i]]=read.delim(paste0("GSE106899_RAW/",fs[i]),row.names = 1)
}
re2=do.call(cbind,re)
class(re2)

## [1] "data.frame"

exp=as.matrix(re2)

#strsplit(fs, "_", simplify=T) 是将字符串 fs 按下划线 _ 分割，并将结果转化为简单数据结构（如向量）。simplify=T 参数是为了将结果转化为数据结构，而不是列表。然后取分割后向量的第一列

#第二种方法
#re3=lapply(fs,function(f){
#  read.delim(paste0("GSE106899_RAW/",f),row.names = 1)
#})
#re4=do.call(cbind,re3)
#以上是将列表中的元素合并成一个数据框
#re=list()
# 3.基因过滤
##需要过滤一下那些在很多样本里表达量都为0或者表达量很低的基因。过滤标准不唯一。
#过滤之前基因数量：
#### 常用过滤标准1：
#仅去除在所有样本里表达量都为零的基因
exp4 = exp[rowSums(exp)>0,]
nrow(exp4)

## [1] 12277

#### 常用过滤标准2(推荐)：
#仅保留在一半以上样本里表达的基因
exp5 = exp4[apply(exp4, 1, function(x) sum(x > 0) > 0.5*ncol(exp4)), ]
exp6 = exp5
#在R语言中，若要把fun应用到x的每一列，margin参数应该设置为1。 #1，函数会应用于矩阵的每一列（即，横向）。 #2，函数会应用于矩阵的每一行（即，纵向）。
#常用的过滤基因的标准

### 4.分组信息获取 一般使control在前 treat在后 要变成因子型 才具有顺序
#header=F参数表示该文件的第一行不是列名，即该文件没有标题行。
a=read.delim("Group.txt",header=F)
library(GEOquery)
a1=getGEO("GSE106899",destdir =".",getGPL=F)

## Found 1 file(s)

## GSE106899_series_matrix.txt.gz

## Using locally cached version: ./GSE106899_series_matrix.txt.gz

library(stringr)
colnames(exp)=str_split(fs,"_",simplify=T)[,1]
identical(a$V1,colnames(exp))

## [1] TRUE

class(a$V2)

## [1] "character"

Group=ifelse(str_detect(a$V2,"Control"),"control","tumor")
Group=factor(Group,levels=c("control","tumor"))
table(Group)

## Group
## control   tumor 
##      22      44

group=Group

### 5.保存数据
save(exp6,group,proj,file = paste0(proj,".Rdata"))