TCGA分析-数据下载2

2023-11-01 17:04:12 浏览数 (2)

是TCGA分析-数据整理-2的上一步

https://cloud.tencent.com/developer/article/2353514


title: "xiaohe"

output: html_document

date: "2023-11-01"


R Markdown

代码语言:text复制
### 小何开始运行

#1.数据下载 #从Gene Expression Omnibus (GEO)数据库下载数据
library(GEOquery)
proj = "GSE106899"
eSet <- getGEO(proj, destdir = '.', getGPL = F)
代码语言:txt复制
## Found 1 file(s)
代码语言:txt复制
## GSE106899_series_matrix.txt.gz
代码语言:txt复制
## Using locally cached version: ./GSE106899_series_matrix.txt.gz
代码语言:text复制
class(eSet)
代码语言:txt复制
## [1] "list"
代码语言:text复制
length(eSet)
代码语言:txt复制
## [1] 1
代码语言:text复制
#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)
代码语言:txt复制
## [1] "data.frame"
代码语言:text复制
exp=as.matrix(re2)
代码语言:txt复制
代码语言:text复制
#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)
代码语言:txt复制
## [1] 12277
代码语言:text复制
#### 常用过滤标准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)
代码语言:txt复制
## Found 1 file(s)
代码语言:txt复制
## GSE106899_series_matrix.txt.gz
代码语言:txt复制
## Using locally cached version: ./GSE106899_series_matrix.txt.gz
代码语言:text复制
library(stringr)
colnames(exp)=str_split(fs,"_",simplify=T)[,1]
identical(a$V1,colnames(exp))
代码语言:txt复制
## [1] TRUE
代码语言:text复制
class(a$V2)
代码语言:txt复制
## [1] "character"
代码语言:text复制
Group=ifelse(str_detect(a$V2,"Control"),"control","tumor")
Group=factor(Group,levels=c("control","tumor"))
table(Group)
代码语言:txt复制
## Group
## control   tumor 
##      22      44
代码语言:text复制
group=Group

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

0 人点赞