Load library
library(GEOquery)
## Loading required package: Biobase
## Loading required package: BiocGenerics
##
## Attaching package: 'BiocGenerics'
## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs
## The following objects are masked from 'package:base':
##
## anyDuplicated, aperm, append, as.data.frame, basename, cbind,
## colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
## get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
## match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
## Position, rank, rbind, Reduce, rownames, sapply, setdiff, sort,
## table, tapply, union, unique, unsplit, which.max, which.min
## Welcome to Bioconductor
##
## Vignettes contain introductory material; view with
## 'browseVignettes()'. To cite Bioconductor, see
## 'citation("Biobase")', and for packages 'citation("pkgname")'.
## Setting options('download.file.method.GEOquery'='auto')
## Setting options('GEOquery.inmemory.gpl'=FALSE)
library(pheatmap)
library(limma)
##
## Attaching package: 'limma'
## The following object is masked from 'package:BiocGenerics':
##
## plotMA
library(readr)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following object is masked from 'package:Biobase':
##
## combine
## The following objects are masked from 'package:BiocGenerics':
##
## combine, intersect, setdiff, union
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(ggplot2)
library(ggrepel)
setwd("D:\\R.data\\CTCL")
get GEO data
gse<-getGEO("GSE59307")
## Found 1 file(s)
## GSE59307_series_matrix.txt.gz
check how many platform and get unite data for each platform
length(gse)
## [1] 1
gse<-gse[[1]]
sampleInfo<-pData(gse)
features<-fData(gse)
sampleInfo<-select(sampleInfo,"condition:ch1","Stage:ch1")
colnames(sampleInfo)<-c("Group","Stage")
features<-select(features,"Gene Symbol","ENTREZ_GENE_ID")
colnames(features)<-c("Symbol","Entrez_Gene_ID")
features<-tidyr::separate(features,'Symbol',into=c("Symbol_1","Symbol_2"),sep=" /// ")
## Warning: Expected 2 pieces. Additional pieces discarded in 682 rows [258, 259, 335, 404,
## 432, 768, 779, 849, 881, 905, 990, 1130, 1131, 1178, 1191, 1428, 1432, 1794,
## 1922, 1991, ...].
## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 51879 rows [2, 3, 4, 5,
## 7, 8, 9, 10, 11, 12, 13, 14, 15, 17, 18, 19, 20, 21, 22, 23, ...].
Scale (0 - 16) and Normalize
summary(exprs(gse))
## GSM1434743 GSM1434744 GSM1434745 GSM1434746
## Min. : 2.285 Min. : 2.285 Min. : 2.285 Min. : 2.285
## 1st Qu.: 2.287 1st Qu.: 2.285 1st Qu.: 2.285 1st Qu.: 2.285
## Median : 3.142 Median : 3.112 Median : 3.176 Median : 3.150
## Mean : 4.251 Mean : 4.318 Mean : 4.331 Mean : 4.323
## 3rd Qu.: 5.721 3rd Qu.: 5.969 3rd Qu.: 5.974 3rd Qu.: 5.966
## Max. :15.951 Max. :15.911 Max. :15.994 Max. :15.663
## GSM1434747 GSM1434748 GSM1434749 GSM1434750
## Min. : 2.285 Min. : 2.285 Min. : 2.285 Min. : 2.285
## 1st Qu.: 2.285 1st Qu.: 2.285 1st Qu.: 2.285 1st Qu.: 2.285
## Median : 3.110 Median : 3.177 Median : 3.107 Median : 3.196
## Mean : 4.332 Mean : 4.340 Mean : 4.350 Mean : 4.322
## 3rd Qu.: 6.016 3rd Qu.: 6.020 3rd Qu.: 6.037 3rd Qu.: 5.925
## Max. :15.894 Max. :15.912 Max. :15.943 Max. :15.876
## GSM1434751 GSM1434752 GSM1434753 GSM1434754
## Min. : 2.285 Min. : 2.285 Min. : 2.285 Min. : 2.285
## 1st Qu.: 2.300 1st Qu.: 2.285 1st Qu.: 2.285 1st Qu.: 2.285
## Median : 3.284 Median : 3.162 Median : 3.136 Median : 3.139
## Mean : 4.181 Mean : 4.325 Mean : 4.331 Mean : 4.335
## 3rd Qu.: 5.363 3rd Qu.: 5.962 3rd Qu.: 5.992 3rd Qu.: 6.025
## Max. :15.911 Max. :15.950 Max. :15.928 Max. :15.845
## GSM1434755 GSM1434756 GSM1434757 GSM1434758
## Min. : 2.285 Min. : 2.285 Min. : 2.285 Min. : 2.285
## 1st Qu.: 2.285 1st Qu.: 2.286 1st Qu.: 2.286 1st Qu.: 2.286
## Median : 3.165 Median : 3.225 Median : 3.253 Median : 3.240
## Mean : 4.338 Mean : 4.332 Mean : 4.342 Mean : 4.340
## 3rd Qu.: 6.006 3rd Qu.: 5.960 3rd Qu.: 5.995 3rd Qu.: 6.007
## Max. :15.897 Max. :15.950 Max. :15.967 Max. :15.971
## GSM1434759 GSM1434760 GSM1434761 GSM1434762
## Min. : 2.285 Min. : 2.285 Min. : 2.285 Min. : 2.285
## 1st Qu.: 2.285 1st Qu.: 2.285 1st Qu.: 2.286 1st Qu.: 2.285
## Median : 3.213 Median : 3.197 Median : 3.267 Median : 3.226
## Mean : 4.335 Mean : 4.346 Mean : 4.327 Mean : 4.345
## 3rd Qu.: 5.996 3rd Qu.: 6.029 3rd Qu.: 5.953 3rd Qu.: 6.038
## Max. :15.920 Max. :15.965 Max. :15.965 Max. :15.913
## GSM1434763 GSM1434764
## Min. : 2.285 Min. : 2.285
## 1st Qu.: 2.285 1st Qu.: 2.285
## Median : 3.225 Median : 3.224
## Mean : 4.335 Mean : 4.336
## 3rd Qu.: 6.012 3rd Qu.: 5.998
## Max. :15.973 Max. :15.963
boxplot(exprs(gse),outline=FALSE)
exprs(gse)<-normalizeQuantiles(exprs(gse))
boxplot(exprs(gse),outline=FALSE)
pheatmap for corelation matrix
corMatrix<-cor(exprs(gse),use="c")
pheatmap(corMatrix,annotation_col=sampleInfo,show_colnames = FALSE, show_rownames = FALSE)
Pricipal Component Analysis
pca<-prcomp(t(exprs(gse)))
cbind(sampleInfo, pca$x) %>% ggplot(aes(x = PC1, y = PC2, col=Group)) +geom_point()
Differential Expression Analasis
design<-model.matrix(~0+sampleInfo$Group)
design
## sampleInfo$GroupCutaneus Tcell Lymphoma sampleInfo$GroupHealthy
## 1 1 0
## 2 1 0
## 3 1 0
## 4 1 0
## 5 1 0
## 6 1 0
## 7 1 0
## 8 1 0
## 9 1 0
## 10 1 0
## 11 1 0
## 12 1 0
## 13 1 0
## 14 1 0
## 15 0 1
## 16 0 1
## 17 0 1
## 18 0 1
## 19 0 1
## 20 0 1
## 21 0 1
## 22 0 1
## attr(,"assign")
## [1] 1 1
## attr(,"contrasts")
## attr(,"contrasts")$`sampleInfo$Group`
## [1] "contr.treatment"
colnames(design)<-c("CTCL","Normal")
cutoff<-median(exprs(gse))
is_expressed<-exprs(gse)>cutoff
keep<-rowSums(is_expressed)>2
gse<-gse[keep,]
aw<-arrayWeights(exprs(gse),design)
fit<-lmFit(exprs(gse),design,weights=aw)
contrasts<-makeContrasts(CTCL-Normal,levels=design)
fit2<-contrasts.fit(fit,contrasts)
fit2<-eBayes(fit2)
fit2$genes<-features
full_results<-topTable(fit2,number = Inf)
full_results<-tibble::rownames_to_column(full_results,"ID")
write_csv(full_results,path="GSE59307_full_results.CSV")
## Warning: The `path` argument of `write_csv()` is deprecated as of readr 1.4.0.
## ℹ Please use the `file` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
full_results %>% mutate(Significant = adj.P.Val<0.05, abs(logFC)>1) %>% mutate(Rank=1:n(), Label = ifelse(Rank<20,Symbol_1,"")) %>% ggplot(aes(x=logFC, y = B, col=Significant, label = Label)) + geom_point() + geom_text_repel(col="black")
## Warning: ggrepel: 8 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
filter(full_results, adj.P.Val < 0.05, abs(logFC)>1) %>% write_csv(path="GSE59307_filtered_results.csv")