前面的這個GEO轉(zhuǎn)錄組數(shù)據(jù)的下載和整理流程是我在生信技能樹學習的筆記，今天打開公眾號發(fā)現(xiàn)有50個瀏覽量，說明有和我一樣需要友好的轉(zhuǎn)錄組數(shù)據(jù)下載和整理示范。前面的那個筆記很粗糙，可能對于一些新手來說還是有解釋不到位的地方，所以今天我把這個筆記重新整理了一下。
1.數(shù)據(jù)下載 具體流程見：https://mp.weixin.qq.com/s/-cC23K4WZIJhJMY5Usc-bw
##1.1 進入GEO官網(wǎng)https://www.ncbi.nlm.nih.gov/geo/輸入要搜索的GSE編號GSE150392，點擊搜索。?
##1.2先看看是芯片數(shù)據(jù)還是轉(zhuǎn)錄組數(shù)據(jù)，轉(zhuǎn)錄組數(shù)據(jù)才能用這個流程。閱讀摘要，了解實驗設計和分組情況。?
##1.3下載原始counts數(shù)據(jù)?
##1.4將下載的數(shù)據(jù)放入到工作目錄下，方便R讀取。?
proj ="GSE150392"#給項目去個名字，方便識別和區(qū)分，不然時間長了不記得是做的什么項目和數(shù)據(jù)了。
2.讀取和整理數(shù)據(jù)
##2.1給項目命名
proj ="GSE150392"#給項目去個名字，方便識別和區(qū)分，不然時間長了不記得是做的什么項目和數(shù)據(jù)了。
2.2 讀取并查看表達矩陣是否正常
#用read.table()讀取后會發(fā)現(xiàn)數(shù)據(jù)框很亂，需要調(diào)整參數(shù)，所以換一個函數(shù)讀取。dat = read.csv("GSE150392_Cov_Mock_Raw_COUNTS.csv")#讀取表達矩陣并檢查數(shù)據(jù)是否正常。讀取后發(fā)現(xiàn)數(shù)據(jù)是3萬多行，是不是數(shù)據(jù)讀取的不對？直接打開csv文件查看是否正確。可見看見有一些ERCC開頭的基因表達量都是0，后面過濾的的步驟可以把這些過濾掉，但是萬一這些基因中有表達量不為0的基因，可能就沒有過濾掉，影響結(jié)果，因為這些基因是外源基因，所以要直接去除這些基因。x列的數(shù)據(jù)是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開頭的基因是外源基因#發(fā)現(xiàn)dat的X列的entrezid和gene symbol以"_"連接在一起，要將它們分開。a = str_split(dat$X,"_",simplify =T)#理論上應該分成2列，但卻分成了4列。#鼠標左鍵單擊新分成的矩陣a的列名右側(cè)的?可按順序查看里面的內(nèi)容。dat$X[24]#發(fā)現(xiàn)X列的數(shù)據(jù)有多個_（基因的名字中也包含_）,所以拆分成4列了。解決辦法是找一個拆分第1個_的函數(shù)就可以了。#> [1] "ENSG00000002586.20_PAR_Y_CD99"b = str_replace(dat$X,"_","http:///")#只替換第一個“_”,并將其替換成基因名稱中不包含的字符就可以，如“/”或“?”。怕基因中包含“/”或“?”，那就多寫幾個“/”或“?”，即“///”或“???”。str_replace_all是替換所有的字符。試了用5個？替代，但是報錯，應該是？被占用了。如果下載的表達矩陣沒有entrezid和gene symbol以"_"連接在一起的情況，就不需要這個分割步驟。注意：GEO的數(shù)據(jù)有些地方?jīng)]有統(tǒng)一的規(guī)范，所以處理數(shù)據(jù)的流程也不是完全一樣的。b = str_replace(dat$X,"_","http:///") %>%? str_split("http:///",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的行名，因為數(shù)據(jù)框和矩陣不允許行名有重復，而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]#檢查表達矩陣的行名和數(shù)據(jù)是否達正常#>? ? ? ? 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)#表達矩陣的類型是數(shù)據(jù)框。#> [1] "data.frame"exp = as.matrix(exp)#將表達矩陣由數(shù)據(jù)框轉(zhuǎn)換為矩陣。class(exp)#檢查表達矩陣的類型已經(jīng)按預期轉(zhuǎn)換成了矩陣的數(shù)據(jù)類型。#> [1] "matrix" "array"
4.基因過濾
需要過濾一下那些在很多樣本里表達量都為0或者表達量很低的基因。過濾標準不唯一。
過濾之前基因數(shù)量：
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")# 提取分組信息，將樣本名稱的數(shù)字去除，\\指數(shù)字的意思。 注意：這里提取分組信息也是因數(shù)據(jù)而異的，不同的數(shù)據(jù)提取分組信息的代碼不同。Group#檢查分組信息的數(shù)字已按預期去除。#> [1] "Cov"? "Cov"? "Cov"? "Mock" "Mock" "Mock"Group = factor(Group,levels = c("Mock","Cov"))#對照組在前，處理組在后，注意別寫反了。table(Group)#再次檢查是否達到預期。#> Group#> Mock? Cov #>? ? 3? ? 3
6.保存數(shù)據(jù)
save(exp,Group,proj,file = paste0(proj,".Rdata"))
6.導出數(shù)據(jù)
write.csv(exp,file = paste0(proj,".csv"))
其他數(shù)據(jù)下載方法
gdc-client
是從官方網(wǎng)站下載，代碼在gdc文件夾，難度較大，結(jié)果與xena整理的一致，作為補充學習瀏覽一下即可。
GDCRNAtools
http://bioconductor.org/packages/devel/bioc/vignettes/GDCRNATools/inst/doc/GDCRNATools.html
注意事項：
1.不同的GEO數(shù)據(jù)其數(shù)據(jù)的整理代碼不同，要根據(jù)數(shù)據(jù)的實際情況將代碼進行調(diào)整，沒有通用的代碼； 2.要有這個代碼里的一步一檢查的思想，有時候代碼不報錯，但運行的結(jié)果如果沒有達到預期，后續(xù)的分析會錯的離譜； 3.數(shù)據(jù)的整理很重要，前期數(shù)據(jù)整理好了后面的分析就都順暢了。
#致謝生信技能樹和小潔老師！
