gset <- getGEO("GSE42872",destdir = ".",AnnotGPL = F,getGPL = F)
if(F){
  gset <- getGEO("GSE42872",destdir = ".",AnnotGPL = F,getGPL = F)
  save(gset,file  = "GSE42872_eSet.Rdata")
}

if(T){
  gset <- getGEO("GSE42872",destdir = ".",AnnotGPL = F,getGPL = F)
  save(gset,file  ="GSE42872_eSet.Rdata")
}
if(F){
  gset <- getGEO("GSE42872",destdir = ".",AnnotGPL = F,getGPL = F)
  save(gset,file  ="GSE42872_eSet.Rdata")
}

gset <- getGEO("GSE42872",destdir = ".",AnnotGPL = F,getGPL = F)
class(gset)
length(gset)
class(gset[[1]])
gset

a  =  getGEO(file = "GSE42872_series_matrix.txt.gz",AnnotGPL = F,getGPL = F)
class(a)
length(a)
a


if(T){
  gpl <- getGEO("GPL6244",destdir = ".")
  colnames(Table(gpl))
  head(Table(gpl)[,c(1,12)])
  probe2gene = Table(gpl)[,c(1,12)]
  head(probe2gene)
  library(stringr)
  save(probe2gene,file = "probe2gene.Rdata")
}

a  = gset[[1]]
dat = exprs(a)
?exprs
class(a)
class(dat)
dim(dat)
dat[1:4,1:4]
boxplot(dat,las = 2)
##dat中的GSMxxx是不是探針的id?。繛楹萎媌oxplot的時候，這里6個的中值在一條線上數(shù)據(jù)才是可用的呢？每個GSMxxx的列中所有行的值的分布
boxplot(dat,las = 1)
boxplot(dat,las = 2)
boxplot(dat[,1])
boxplot(dat[,1],las = 1)
boxplot(dat[,1:2],las = 1)
##las = 1表示平行，2表示垂直與x軸
pd = pData(a)
class(pd)
dim(pd)
library(stringr)
pd[1:4,1:4]
str_split(pd$title, " ")
class(str_split(pd$title, " "))
class(str_split(pd$title, " ",simplify = T))
str_split(pd$title, " ",simplify = T)
?str_split
group_list = str_split(pd$title," ",simplify  =T)[,4]
class(str_split(pd$title," ",simplify  =T)) #得到的是矩陣
group_list
table(group_list)


gpl <- getGEO("GPL6244",destdir = ".")
##這里如果提前終止該命令，特別是在網(wǎng)速慢的時候會提前終止，得到的GPL6244.soft文件不完全，但是在windows卻也不能刪除它，提示被占用。
gpl_Table <- Table(gpl)
##gpl_Table是data.frame，可見Table()就是提取GPL6244.soft這個文件中前述表達(dá)矩陣中對應(yīng)探針的注釋部分
##在notepad打開GPL6244.soft這個文件，前面是以!或者#開頭的解釋行,包括對后續(xù)正文的列名的注釋，然后是!platform_table_begin，正文，!platform_table_end.
##哪有Table()這個函數(shù)的源碼呀？
dim(gpl_Table)##33297行 12列
dim(dat) ##都是33297行
head(gpl_Table)
head(dat)
head(gpl_Table)[,c(1,12)]
tail(gpl_Table)[,c(1,12)]
probe2gene = gpl_Table[,c(1,12)]
head(probe2gene)
save(probe2gene,file = "probe2gene.Rdata")

dim(probe2gene)
dim(dat)
###既然probe2gene與dat行數(shù)一致，測試一下dat的每個行號是否都在probe2gene$ID中
a = rownames(dat) %in% probe2gene$ID
table(a)
###返回的是 TRUE 33297
##
##有個東西 一直想不起來
##前面有很多行以！或者#開頭表示注釋，然后只是。。。begin，然后是正文，是一個大的表達(dá)矩陣還是啥
##最后有。。。end.......打開一看，不就是這里的GPLxxx.soft文件

##統(tǒng)計(jì)一下probe2gene$category這一列的頻次，所以每一種都代表的啥意思??？
table(probe2gene$category)

library(BiocManager)
BiocManager::install("hugene10sttranscriptcluster.db")

suppressMessages(library(hugene10sttranscriptcluster.db))
ids = toTable(hugene10sttranscriptclusterSYMBOL)

##如何知道要下載這個數(shù)據(jù)集？它跟前面的dat中的行號以及probe2gene中的ID列有啥關(guān)系？

head(ids)
dim(ids)
head(probe2gene)
dim(probe2gene)
##這里得到的probe2gene后面似乎沒有用到啊，那為什么要取它？

##rownames(dat)有13483個是ids$probe2gene中沒有的：FALSE 13483 TRUE 19814
b = rownames(dat) %in% ids$probe_id
table(b)

##反過來,ids$probe_id都在dat的行號或者probe2gene$ID中：TRUE 19814
c = ids$probe_id %in% rownames(dat)
table(c)

colnames(ids) = c("probe_id","symbol") ##不明白這里，ids的colnames就是probe_idh和symbol,為啥這里又再賦值一次？
ids = ids[ids$symbol !="",] ##難道說，symbo這一列還有為空的行，也就是有probe_id，沒有symbol，反正就是排除為空的行！
###下面驗(yàn)證，確實(shí)是的
x1 = data.frame(x = c(1,4,"",23,46),y = c(1,1,1,1,1))
x1
y1 = a[a$x != "",]
y1
z1 = a[a$y != "",]
z1
###但是下面呢，數(shù)據(jù)框還可以依據(jù)邏輯值訪問嘛？
ids_no_symbol = ids$symbol != ""
head(ids_no_symbol)
table(ids_no_symbol)
ids[TRUE,]
ids[FALSE,]
ids[,FALSE]

ids = ids[ids$symbol != "",] ##1步
##跟下面%in%一樣，返回的是邏輯值，然后用邏輯值訪問數(shù)據(jù)框
##我已開始還以為是挨個判斷symbol為不為空，不為空的話就訪問這個symbol，然后取出對應(yīng)的probe2gene添加到新的ids中的行
dim(ids)
ids = ids[ids$probe_id %in% rownames(dat),] ##2步
##只取邏輯值為TRUE的行嘛？前面已經(jīng)table了，所有的都為TURE
dim(ids)

dat[1:4,1:4]
head(ids)
dat = dat[ids$probe_id,] ##3步
##1先去除symbol為空的行，然后2取ids中probe_id和dat中rownames相同的行;然后再3將dat中非ids$probe_id的行去除
dim(dat) ##dat只有19814行了

##通過行名訪問矩陣，加上引號
##一開始用的是dat["7892501",],提示下標(biāo)出界
##我就想訪問矩陣有哪些方式：可不可以通過行號訪問??？代碼是怎樣的？
##然后百度發(fā)現(xiàn)可以，就又試了下上面那句，還是不行，那是不是這個行號沒了，就head了dat，挑第一二行試了下，可以！
dat["7892501",]
dat["7896759",]
dat[c("7896759","7896761"),]

ids$median = apply(dat,1,median)
##這時候dat和ids的行就一一對應(yīng)了嘛？直接添加一列median是匹配的嗎？
##當(dāng)然是的，因?yàn)楝F(xiàn)在的dat是由ids$probe_id訪問原有的dat的行號得到的,順序應(yīng)該和ids$probe_id是一致
##前面對GSMxxx的列話boxplot，這里又對每一行取中值。。。不懂啊
##還有，直接在現(xiàn)在的dat中添加median這一列不行嘛？
dim(ids)
head(ids)
ids = ids[order(ids$symbol,ids$median,decreasing = T),]
##Jimmy說symbol按照median從大到小排列，但是從結(jié)果看怎么是median按照symbol倒序排列？？？
##為什么要排序呢？
dim(ids)
head(ids)
##這里order()之后的結(jié)果有些迷惑。。。下面見示例
##按示例
t1 = data.frame(x  =c("a","b","c","d","e"),y = c(8,2,4,10,6))
t1
t2  = t1[order(t1$x,t1$y,decreasing = T),]
t2


##跑了這一步，ids就已經(jīng)變了，所以ids中duplicated的子集是什么樣的？不加！試試
dim(ids)
dim(dat)
dat = dat[ids$probe_id,]
dim(dat)
dat[1:4,1:4]
ids[1:4,1:3] 
rownames(dat) = ids$symbol
##到這里可以看出,經(jīng)過dat=dat[ids$probe_id,]這一步，dat中的行號(也就是ids中的probe_id)
##也跟著ids的順序更改了，且一一對應(yīng),因此才可以有rownames(dat) = ids$symbol這一句。
##不然，應(yīng)該不能直接把ids中symbol直接賦值dat的行名
dat[1:4,1:4]
ls()
save(dat,ids,group_list,gset,probe2gene,pd,file = "step1-output-version2.Rdata")

參考：https://github.com/jmzeng1314/GEO/tree/master/GSE42872_main

總結(jié)

以上是生活随笔為你收集整理的step-1的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。