step-1

 

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啊？为何画boxplot的时候，这里6个的中值在一条线上数据才是可用的呢？每个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 = ".")
##这里如果提前终止该命令，特别是在网速慢的时候会提前终止，得到的GPL6244.soft文件不完全，但是在windows却也不能删除它，提示被占用。
gpl_Table <- Table(gpl)
##gpl_Table是data.frame，可见Table()就是提取GPL6244.soft这个文件中前述表达矩阵中对应探针的注释部分
##在notepad打开GPL6244.soft这个文件，前面是以!或者#开头的解释行,包括对后续正文的列名的注释，然后是!platform_table_begin，正文，!platform_table_end.
##哪有Table()这个函数的源码呀？
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行数一致，测试一下dat的每个行号是否都在probe2gene$ID中
a = rownames(dat) %in% probe2gene$ID
table(a)
###返回的是 TRUE 33297
##
##有个东西 一直想不起来
##前面有很多行以！或者#开头表示注释，然后只是。。。begin，然后是正文，是一个大的表达矩阵还是啥
##最后有。。。end.......打开一看，不就是这里的GPLxxx.soft文件

##统计一下probe2gene$category这一列的频次，所以每一种都代表的啥意思啊？
table(probe2gene$category)

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

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

##如何知道要下载这个数据集？它跟前面的dat中的行号以及probe2gene中的ID列有啥关系？

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，反正就是排除为空的行！
###下面验证，确实是的
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
###但是下面呢，数据框还可以依据逻辑值访问嘛？
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%一样，返回的是逻辑值，然后用逻辑值访问数据框
##我已开始还以为是挨个判断symbol为不为空，不为空的话就访问这个symbol，然后取出对应的probe2gene添加到新的ids中的行
dim(ids)
ids = ids[ids$probe_id %in% rownames(dat),] ##2步
##只取逻辑值为TRUE的行嘛？前面已经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",],提示下标出界
##我就想访问矩阵有哪些方式：可不可以通过行号访问啊？代码是怎样的？
##然后百度发现可以，就又试了下上面那句，还是不行，那是不是这个行号没了，就head了dat，挑第一二行试了下，可以！
dat["7892501",]
dat["7896759",]
dat[c("7896759","7896761"),]

ids$median = apply(dat,1,median)
##这时候dat和ids的行就一一对应了嘛？直接添加一列median是匹配的吗？
##当然是的，因为现在的dat是由ids$probe_id访问原有的dat的行号得到的,顺序应该和ids$probe_id是一致
##前面对GSMxxx的列话boxplot，这里又对每一行取中值。。。不懂啊
##还有，直接在现在的dat中添加median这一列不行嘛？
dim(ids)
head(ids)
ids = ids[order(ids$symbol,ids$median,decreasing = T),]
##Jimmy说symbol按照median从大到小排列，但是从结果看怎么是median按照symbol倒序排列？？？
##为什么要排序呢？
dim(ids)
head(ids)
##这里order()之后的结果有些迷惑。。。下面见示例
##按示例
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就已经变了，所以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
##到这里可以看出,经过dat=dat[ids$probe_id,]这一步，dat中的行号(也就是ids中的probe_id)
##也跟着ids的顺序更改了，且一一对应,因此才可以有rownames(dat) = ids$symbol这一句。
##不然，应该不能直接把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