|
今天這個(gè)帖子，對(duì)我來說很有意義，我在最后要介紹一位心中的男神。 我每天努力學(xué)習(xí)就是為了消除心中的疑惑。 幾年前，當(dāng)我跟著別人的代碼跑了個(gè)流程后，突破不了的是兩個(gè)事情 那個(gè)時(shí)候，R語言基礎(chǔ)很差，處理不了數(shù)據(jù)，很有挫敗感，所以就停止了R語言的學(xué)習(xí)。 直到我碰到了那個(gè)男人。 現(xiàn)在我解決了這個(gè)問題，分享給大家。 第一種，我們可以直接用平臺(tái)的數(shù)據(jù)進(jìn)入官網(wǎng) https://www.ncbi.nlm./geo/  
知道平臺(tái)是GPL6244  這時(shí)候我們進(jìn)入R語言，用GEOquery中的getGEO可以獲得探針和基因名的信息 網(wǎng)絡(luò)不好的可以直接略過 library(GEOquery)
GPL6244 <>'GPL6244',destdir ='.')
轉(zhuǎn)換成數(shù)據(jù)框形式,有3萬行，12列 GPL6244_anno <->Table(GPL6244)
查看內(nèi)容，我們發(fā)現(xiàn)基因名稱藏在了gene_assignment這一列的中間  所以我們要把他和第一列id提取出來 library(dplyr)
library(tidyr)
probe2symbol_df <- gpl6244_anno="" %="">%
select(ID,gene_assignment) %>%
filter(gene_assignment != '---') %>%
separate(gene_assignment,c('drop','symbol'),sep='//') %>%
select(-drop)
看一下，數(shù)據(jù)已經(jīng)被提取出來了。  假如getGEO這一步網(wǎng)絡(luò)不好呢 library(GEOquery)
GPL6244 <>'GPL6244',destdir ='.')
我們?cè)谶@個(gè)一開始的這個(gè)頁面下載平臺(tái)的soft文件 點(diǎn)擊soft文件  下載解壓 然后用data.table這個(gè)包中的fread即可閱讀進(jìn)來，注意，skip這個(gè)參數(shù)十分重要?。?/span> GPL6244_anno <>'./GSE42872_family.soft/GSE42872_family.soft',skip ='ID')
得到了GPL6244_anno,我們又可以運(yùn)行下面的代碼提出探針和基因名稱對(duì)應(yīng)的關(guān)系了。 library(dplyr)
library(tidyr)
probe2symbol_df <- gpl6244_anno="" %="">%
select(ID,gene_assignment) %>%
filter(gene_assignment != '---') %>%
separate(gene_assignment,c('drop','symbol'),sep='//') %>%
select(-drop)
第二種，使用R包R有很多注釋包，首先我們要知道什么平臺(tái)用什么R包。 我這里整理了一個(gè)對(duì)應(yīng)關(guān)系  讀入R語言 platformMap <->'platformMap.txt')
為了方便查詢，我寫了以下的代碼，每次只要修改平臺(tái)的名稱即可 index <->'GPL6244'
paste0(platformMap$bioc_package[grep(index,platformMap$gpl)],'.db')
輸出的是， 'hugene10sttranscriptcluster.db'
安裝并加載這個(gè)包 BiocInstaller::biocLite('hugene10sttranscriptcluster.db')
library(hugene10sttranscriptcluster.db)
使用toTable函數(shù)找到對(duì)應(yīng)關(guān)系 probe2symbol_df <->get('hugene10sttranscriptclusterSYMBOL'))
那探針I(yè)D轉(zhuǎn)換的兩種方法就講完了 正式進(jìn)行探針I(yè)D之間的轉(zhuǎn)換我們現(xiàn)在用的是矩陣文件exprSet,  我們要達(dá)到的效果是這個(gè)樣的  其實(shí)很簡(jiǎn)單，就是把兩個(gè)數(shù)據(jù)框按照探針名稱合并就可以了， 同時(shí)在這里，我們還需要把重復(fù)的探針去掉，方法就是計(jì)算每一個(gè)探針的平均值，留下最大的那個(gè)。 先調(diào)整兩個(gè)數(shù)據(jù)框的探針名字和類型，方便合并 names(exprSet)[1] <->1]
exprSet$probe_id <->as.character(exprSet$probe_id)
然后， 重點(diǎn)就來了，下面的代碼恢弘大氣，一次性完成這兩個(gè)事情。library(dplyr)
exprSet <- exprset="" %="">%
inner_join(probe2symbol_df,by='probe_id') %>% #合并探針的信息
select(-probe_id) %>% #去掉多余信息
select(symbol, everything()) %>% #重新排列，
mutate(rowMean =rowMeans(.[grep('GSM', names(.))])) %>% #求出平均數(shù)(這邊的.真的是畫龍點(diǎn)睛)
arrange(desc(rowMean)) %>% #把表達(dá)量的平均值按從大到小排序
distinct(symbol,.keep_all = T) %>% # symbol留下第一個(gè)
select(-rowMean) %>% #反向選擇去除rowMean這一列
tibble::column_to_rownames(colnames(.)[1]) # 把第一列變成行名并刪除
真的很skr！這個(gè)代碼的產(chǎn)生過程，以后再說。 這里的%>% 相當(dāng)于Linux中的xargs， Hadley Wickham真的是天才！ 如果不是他的這個(gè)神包tidyverse以及健明給的鼓勵(lì) 我可能到現(xiàn)在還入不了R的門！
|
