-
Notifications
You must be signed in to change notification settings - Fork 7
/
Copy path1.1 download_and_process_GEO_dataset.R
158 lines (121 loc) · 4.37 KB
/
1.1 download_and_process_GEO_dataset.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
# Step0 Before starting your project --------------------------------------
## Remove everything in the working environment, not including loaded libraries.
rm(list = objects( all = TRUE ))
if (!is.null( dev.list() )) dev.off()
clearhistory <- function() {
write( "", file = ".blank" )
loadhistory( ".blank" )
unlink( ".blank" )
}
clearhistory()
## basecal packages
sysPackages <- (.packages())
## data.frame(..., row.names = NULL, check.rows = FALSE,
## check.names = TRUE, fix.empty.names = TRUE,
## stringsAsFactors = default.stringsAsFactors())
options( stringsAsFactors = FALSE )
## Now winInet not supported for use in service, but the default setting of
## download.file.method is "wininet".
## If your system support "libcurl", set the downloaded method to libcurl.
if ( capabilities( "libcurl" ) == T ) {
options( download.file.method = "libcurl" )
}
options()$download.file.method
## Change the library location of the packages
## Even your updated your R, you can still use your packages.
.libPaths( c( "G:/R-packages",
"C:/Program Files/R/R-3.5.2/library") )
.libPaths()
# Step1 download GEO dateset by GEOquery ----------------------------------
## Basic data
GSE_ID <- ''
GSE_file <- paste('./data/', GSE_ID, '.Rdata', sep = "", collapse = NULL)
library( "GEOquery" )
if (!file.exists( GSE_file )) {
## dir.create("./data/")
download_GEO <- function(GSE_ID) {
options( download.file.method.GEOquery = 'libcurl' )
## dir.create("./raw_data/")
gset <- getGEO( GSE_ID, getGPL = T, destdir = "./raw_data/" )
save( gset, file = GSE_file )
}
download_GEO( GSE_ID )
}
# Step2 Data loading ---------------------------------------------------
load( GSE_file )
class( gset )
length( gset )
exprSet <- gset[[1]]
str( exprSet, max.level = 2 )
## assayData <- exprSet@assayData$exprs[1:5,1:6]
assayData <- exprs(exprSet)
dim(assayData)
assayData[,1:4]
## phenoData <- exprSet@phenoData@data
phenoData <- pData(exprSet)
dim(phenoData)
phenoData[1:2,1:6]
# Step3 Grouping by special clinical information --------------------------
pheno_num <- c()
invisible(
lapply(1:ncol(phenoData),
function(col_num){
## Assume that the classification project is between 2 and 4
if (1 < dim(table(phenoData[,col_num])) &
dim(table(phenoData[,col_num])) < 5) {
pheno_num <<- append(pheno_num, col_num, after = length(pheno_num))
}
}
)
)
View(phenoData[, pheno_num])
names(phenoData[, pheno_num])
## choose special pheno
group_list <- phenoData[, "neo-adjuvant chemotherapy:ch1"]
table(group_list)
group_list <- ifelse( group_list == 'resistant',
'resistant','sensitive')
names(group_list) <- rownames(phenoData)
newAssayData <- assayData[, names(group_list)]
dim(newAssayData)
newAssayData[1:5, 1:6]
if (!file.exists( "./data/AssayData_by_group.Rdata" )) {
save( newAssayData, group_list, file = "./data/AssayData_by_group.Rdata" )
}
# Step4 Filt gene ---------------------------------------------------------
load( './data/AssayData_by_group.Rdata' )
load( './data/GSE108565.Rdata' )
exprSet <- gset[[1]]
featureData = exprSet@featureData@data
dim( featureData )
colnames( featureData )
View( featureData )
## b <- apply(featureData, 1,
## function(i){
## featureData[i, 1] == featureData[i, 2]
## })
## table(b)
featureData$max <- apply(newAssayData, 1, max)
featureData <- featureData[order(featureData$ID,
featureData$max,
decreasing = T), ]
dim( featureData )
featureData <- featureData[!duplicated(featureData$ID), ]
dim( featureData )
AssayData <- newAssayData[featureData$ID, ]
dim(AssayData)
AssayData[1:5, 1:6]
# AssayData = log2(AssayData)
# dim(AssayData)
# AssayData[1:5,1:6]
# group_list <- unname(group_list)
save(AssayData, group_list, file = './data/final_AssayData.Rdata')
# Step5 Remove new loaded packages ----------------------------------------
allPackages <- (.packages())
newPackages <- setdiff( allPackages, sysPackages )
lapply( newPackages,
function(package) {
package <- paste('package:', package, sep = "", collapse = NULL)
detach( package, character.only = TRUE )
}
)