前面的这个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.数据的整理很重要，前期数据整理好了后面的分析就都顺畅了。
#致谢生信技能树和小洁老师！
