GEO转录组数据下载GSE150392

前面的这个GEO转录组数据的下载和整理流程是我在生信技能树学习的笔记,今天打开公众号发现有50个浏览量,说明有和我一样需要友好的转录组数据下载和整理示范。前面的那个笔记很粗糙,可能对于一些新手来说还是有解释不到位的地方,所以今天我把这个笔记重新整理了一下。

1.数据下载 具体流程见:https://mp.weixin.qq.com/s/-cC23K4WZIJhJMY5Usc-bw

##1.1 进入GEO官网https://www.ncbi.nlm.nih.gov/geo/输入要搜索的GSE编号GSE150392,点击搜索。 

##1.2先看看是芯片数据还是转录组数据,转录组数据才能用这个流程。阅读摘要,了解实验设计和分组情况。 

##1.3下载原始counts数据 

##1.4将下载的数据放入到工作目录下,方便R读取。 

proj ="GSE150392"#给项目去个名字,方便识别和区分,不然时间长了不记得是做的什么项目和数据了。

2.读取和整理数据

##2.1给项目命名

proj ="GSE150392"#给项目去个名字,方便识别和区分,不然时间长了不记得是做的什么项目和数据了。

2.2 读取并查看表达矩阵是否正常

#用read.table()读取后会发现数据框很乱,需要调整参数,所以换一个函数读取。dat = read.csv("GSE150392_Cov_Mock_Raw_COUNTS.csv")#读取表达矩阵并检查数据是否正常。读取后发现数据是3万多行,是不是数据读取的不对?直接打开csv文件查看是否正确。可见看见有一些ERCC开头的基因表达量都是0,后面过滤的的步骤可以把这些过滤掉,但是万一这些基因中有表达量不为0的基因,可能就没有过滤掉,影响结果,因为这些基因是外源基因,所以要直接去除这些基因。x列的数据是ENSG00000000003.15_TSPAN6,既包括ensmbol和symbol,需要提取symbol.

library(stringr)#> Warning: package 'stringr' was built under R version 4.1.2k1 = !str_detect(dat$X,"ERCC");table(k1)#ERCC开头的基因是外源基因,看看有多少个外源基因。#> k1#> FALSE  TRUE #>  104 36837dat = dat[k1,]#去除ERCC开头的基因是外源基因#发现dat的X列的entrezid和gene symbol以"_"连接在一起,要将它们分开。a = str_split(dat$X,"_",simplify =T)#理论上应该分成2列,但却分成了4列。#鼠标左键单击新分成的矩阵a的列名右侧的▶可按顺序查看里面的内容。dat$X[24]#发现X列的数据有多个_(基因的名字中也包含_),所以拆分成4列了。解决办法是找一个拆分第1个_的函数就可以了。#> [1] "ENSG00000002586.20_PAR_Y_CD99"b = str_replace(dat$X,"_","///")#只替换第一个“_”,并将其替换成基因名称中不包含的字符就可以,如“/”或“?”。怕基因中包含“/”或“?”,那就多写几个“/”或“?”,即“///”或“???”。str_replace_all是替换所有的字符。试了用5个?替代,但是报错,应该是?被占用了。如果下载的表达矩阵没有entrezid和gene symbol以"_"连接在一起的情况,就不需要这个分割步骤。注意:GEO的数据有些地方没有统一的规范,所以处理数据的流程也不是完全一样的。b = str_replace(dat$X,"_","///") %>%  str_split("///",simplify =T)dat$X[24]#检查并确认基因里的“_”并没有被///替代,别将基因里的“_”都替代了。#> [1] "ENSG00000002586.20_PAR_Y_CD99"exp = dat[,2:7]#提取dat里的第2-7列。exp[1:4,1:4]#>  Cov1 Cov2 Cov3 Mock1#> 1 1270 1013  848  4571#> 2  121  125  49  3220#> 3 1023 1045  674  2404#> 4 1049  713  692  854#rownames(dat) = b[,2]#将b的基因名称作为dat的行名,因为数据框和矩阵不允许行名有重复,而b的genen symbol是有重复的,所以运行这个代码会报错,运行这个代码之前要先去重。k2 = !duplicated(b[,2]);table(k2)#查看b第2列的基因名称有多少不重复和重复的。#> k2#> FALSE  TRUE #>    23 36814exp = dat[,2:7]nrow(b)#都是36849行。#> [1] 36837nrow(exp)#都是36849行。#> [1] 36837b = b[k2,]exp = exp[k2,]rownames(exp) = b[,2]exp[1:4,1:4]#检查表达矩阵的行名和数据是否达正常#>        Cov1 Cov2 Cov3 Mock1#> TSPAN6 1270 1013  848  4571#> TNMD    121  125  49  3220#> DPM1  1023 1045  674  2404#> SCYL3  1049  713  692  854class(exp)#表达矩阵的类型是数据框。#> [1] "data.frame"exp = as.matrix(exp)#将表达矩阵由数据框转换为矩阵。class(exp)#检查表达矩阵的类型已经按预期转换成了矩阵的数据类型。#> [1] "matrix" "array"

4.基因过滤

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

过滤之前基因数量:

nrow(exp)#> [1] 36814

常用过滤标准1:

仅去除在所有样本里表达量都为零的基因

exp1 = exp[rowSums(exp)>0,]nrow(exp1)#> [1] 35210

常用过滤标准2(推荐):

仅保留在一半以上样本里表达的基因

exp = exp[apply(exp,1,function(x) sum(x >0) >0.5*ncol(exp)), ]nrow(exp)#> [1] 28197

5.分组信息获取

colnames(exp)#> [1] "Cov1"  "Cov2"  "Cov3"  "Mock1" "Mock2" "Mock3"Group = str_remove(colnames(exp),"\\d")# 提取分组信息,将样本名称的数字去除,\\指数字的意思。 注意:这里提取分组信息也是因数据而异的,不同的数据提取分组信息的代码不同。Group#检查分组信息的数字已按预期去除。#> [1] "Cov"  "Cov"  "Cov"  "Mock" "Mock" "Mock"Group = factor(Group,levels = c("Mock","Cov"))#对照组在前,处理组在后,注意别写反了。table(Group)#再次检查是否达到预期。#> Group#> Mock  Cov #>    3    3

6.保存数据

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

6.导出数据

write.csv(exp,file = paste0(proj,".csv"))

其他数据下载方法

gdc-client

是从官方网站下载,代码在gdc文件夹,难度较大,结果与xena整理的一致,作为补充学习浏览一下即可。

GDCRNAtools

http://bioconductor.org/packages/devel/bioc/vignettes/GDCRNATools/inst/doc/GDCRNATools.html

注意事项:

1.不同的GEO数据其数据的整理代码不同,要根据数据的实际情况将代码进行调整,没有通用的代码; 2.要有这个代码里的一步一检查的思想,有时候代码不报错,但运行的结果如果没有达到预期,后续的分析会错的离谱; 3.数据的整理很重要,前期数据整理好了后面的分析就都顺畅了。

#致谢生信技能树和小洁老师!

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
平台声明:文章内容(如有图片或视频亦包括在内)由作者上传并发布,文章内容仅代表作者本人观点,简书系信息发布平台,仅提供信息存储服务。

推荐阅读更多精彩内容