Differential Expression Analysis
Van-Tuy Nguyen
2024-04-25
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':
##
## plotMAlibrary(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, unionlibrary(ggplot2)
library(ggrepel)
setwd("D:\\R.data\\CTCL")get GEO data
gse<-getGEO("GSE9479")## Found 1 file(s)## GSE9479_series_matrix.txt.gzcheck how many platform and get unite data for each platform
length(gse)## [1] 1gse<-gse[[1]]
sampleInfo<-pData(gse)
features<-fData(gse)view sampleInfo and features
View(sampleInfo)
View(features)modify sampleInfo and features
sampleInfo<-select(sampleInfo,"title","characteristics_ch1")
colnames(sampleInfo)<-c("Title","Stage")
features<-select(features,"Gene Symbol","Entrez Gene")
colnames(features)<-c("Symbol","Entrez_Gene_ID")check for Scale (0 - 16) and Normalize
summary(exprs(gse))## GSM240567 GSM240568 GSM240569 GSM240570
## Min. : 16.66 Min. : 18.22 Min. : 17.43 Min. : 16.43
## 1st Qu.: 35.90 1st Qu.: 36.19 1st Qu.: 42.74 1st Qu.: 38.43
## Median : 119.03 Median : 118.46 Median : 97.28 Median : 116.12
## Mean : 851.98 Mean : 852.64 Mean : 747.26 Mean : 820.07
## 3rd Qu.: 602.20 3rd Qu.: 601.56 3rd Qu.: 469.22 3rd Qu.: 565.10
## Max. :36682.15 Max. :37164.06 Max. :32760.59 Max. :37793.04
## GSM240571 GSM240572 GSM240573 GSM240574
## Min. : 15.75 Min. : 16.80 Min. : 17.66 Min. : 16.78
## 1st Qu.: 36.97 1st Qu.: 37.94 1st Qu.: 37.59 1st Qu.: 38.48
## Median : 120.18 Median : 118.26 Median : 121.95 Median : 116.65
## Mean : 837.71 Mean : 835.44 Mean : 820.14 Mean : 831.75
## 3rd Qu.: 587.93 3rd Qu.: 582.42 3rd Qu.: 598.40 3rd Qu.: 583.07
## Max. :36466.94 Max. :38936.93 Max. :38227.76 Max. :38949.55
## GSM240575 GSM240576 GSM240577 GSM240578
## Min. : 18.31 Min. : 18.16 Min. : 18.30 Min. : 18.64
## 1st Qu.: 38.07 1st Qu.: 37.53 1st Qu.: 36.76 1st Qu.: 37.34
## Median : 119.33 Median : 118.57 Median : 120.81 Median : 120.45
## Mean : 829.69 Mean : 830.12 Mean : 848.02 Mean : 829.79
## 3rd Qu.: 584.05 3rd Qu.: 591.91 3rd Qu.: 584.48 3rd Qu.: 593.42
## Max. :37529.80 Max. :36589.35 Max. :37222.34 Max. :37810.38
## GSM240579 GSM240580 GSM240581 GSM240582
## Min. : 17.92 Min. : 19.06 Min. : 17.38 Min. : 16.23
## 1st Qu.: 37.05 1st Qu.: 36.82 1st Qu.: 36.58 1st Qu.: 36.64
## Median : 120.63 Median : 121.18 Median : 123.15 Median : 121.60
## Mean : 848.24 Mean : 843.98 Mean : 839.95 Mean : 840.39
## 3rd Qu.: 588.75 3rd Qu.: 600.54 3rd Qu.: 586.18 3rd Qu.: 602.24
## Max. :39606.90 Max. :37072.75 Max. :36305.89 Max. :38640.12
## GSM240583 GSM240584 GSM240585 GSM240586
## Min. : 18.08 Min. : 17.85 Min. : 17.72 Min. : 17.17
## 1st Qu.: 36.17 1st Qu.: 36.25 1st Qu.: 36.23 1st Qu.: 36.41
## Median : 123.66 Median : 118.50 Median : 117.73 Median : 108.18
## Mean : 852.33 Mean : 838.03 Mean : 848.34 Mean : 856.97
## 3rd Qu.: 595.91 3rd Qu.: 571.09 3rd Qu.: 581.53 3rd Qu.: 609.28
## Max. :38246.14 Max. :36758.39 Max. :38120.46 Max. :38524.72
## GSM240587 GSM240588 GSM240589 GSM240590
## Min. : 17.73 Min. : 19.18 Min. : 17.28 Min. : 16.87
## 1st Qu.: 36.29 1st Qu.: 38.44 1st Qu.: 37.32 1st Qu.: 37.11
## Median : 124.00 Median : 118.62 Median : 119.81 Median : 123.00
## Mean : 831.62 Mean : 806.17 Mean : 839.82 Mean : 825.45
## 3rd Qu.: 586.27 3rd Qu.: 551.07 3rd Qu.: 585.28 3rd Qu.: 585.07
## Max. :38050.30 Max. :38421.16 Max. :40839.20 Max. :37785.62
## GSM240591 GSM240592 GSM240593 GSM240594
## Min. : 17.39 Min. : 18.86 Min. : 19.12 Min. : 19.17
## 1st Qu.: 37.45 1st Qu.: 39.21 1st Qu.: 39.66 1st Qu.: 38.25
## Median : 121.58 Median : 119.33 Median : 115.39 Median : 120.40
## Mean : 821.32 Mean : 802.20 Mean : 791.86 Mean : 836.64
## 3rd Qu.: 577.58 3rd Qu.: 544.94 3rd Qu.: 543.15 3rd Qu.: 585.85
## Max. :37834.10 Max. :38498.11 Max. :37673.82 Max. :38817.71
## GSM240595 GSM240596 GSM240597 GSM240598
## Min. : 19.71 Min. : 17.73 Min. : 17.07 Min. : 18.31
## 1st Qu.: 37.52 1st Qu.: 37.43 1st Qu.: 38.41 1st Qu.: 36.30
## Median : 123.35 Median : 120.01 Median : 117.91 Median : 115.70
## Mean : 838.06 Mean : 824.96 Mean : 789.90 Mean : 845.46
## 3rd Qu.: 584.50 3rd Qu.: 578.03 3rd Qu.: 553.75 3rd Qu.: 599.34
## Max. :38199.49 Max. :34616.28 Max. :39156.83 Max. :37683.90
## GSM240599 GSM240600 GSM240601 GSM240602
## Min. : 19.56 Min. : 18.42 Min. : 18.92 Min. : 18.34
## 1st Qu.: 38.79 1st Qu.: 36.41 1st Qu.: 36.41 1st Qu.: 38.82
## Median : 109.65 Median : 118.26 Median : 122.67 Median : 106.12
## Mean : 793.54 Mean : 853.37 Mean : 847.02 Mean : 818.31
## 3rd Qu.: 530.05 3rd Qu.: 592.71 3rd Qu.: 590.16 3rd Qu.: 556.04
## Max. :34438.39 Max. :37926.33 Max. :40003.97 Max. :34198.19
## GSM240603 GSM240604 GSM240605 GSM240606
## Min. : 18.20 Min. : 16.46 Min. : 17.94 Min. : 17.69
## 1st Qu.: 36.76 1st Qu.: 37.05 1st Qu.: 36.85 1st Qu.: 37.33
## Median : 120.70 Median : 117.68 Median : 124.13 Median : 118.31
## Mean : 839.83 Mean : 832.92 Mean : 844.01 Mean : 833.55
## 3rd Qu.: 578.81 3rd Qu.: 577.05 3rd Qu.: 595.23 3rd Qu.: 587.45
## Max. :37834.10 Max. :36625.31 Max. :39471.64 Max. :37127.23
## GSM240607 GSM240608 GSM240609 GSM240610
## Min. : 17.83 Min. : 17.53 Min. : 18.51 Min. : 19.23
## 1st Qu.: 41.92 1st Qu.: 37.69 1st Qu.: 40.16 1st Qu.: 38.92
## Median : 102.56 Median : 118.40 Median : 112.49 Median : 115.89
## Mean : 764.16 Mean : 828.33 Mean : 795.31 Mean : 809.00
## 3rd Qu.: 482.66 3rd Qu.: 579.26 3rd Qu.: 527.25 3rd Qu.: 557.19
## Max. :32725.14 Max. :39678.33 Max. :37950.77 Max. :38845.15
## GSM240611 GSM240612 GSM240613 GSM240614
## Min. : 19.55 Min. : 16.51 Min. : 17.31 Min. : 17.05
## 1st Qu.: 41.92 1st Qu.: 37.38 1st Qu.: 35.77 1st Qu.: 35.85
## Median : 109.56 Median : 120.42 Median : 123.77 Median : 122.27
## Mean : 767.50 Mean : 834.86 Mean : 831.01 Mean : 847.48
## 3rd Qu.: 500.00 3rd Qu.: 584.04 3rd Qu.: 595.33 3rd Qu.: 590.39
## Max. :38331.67 Max. :39395.87 Max. :36800.89 Max. :39185.63
## GSM240615 GSM240616 GSM240617 GSM240618
## Min. : 19.05 Min. : 19.08 Min. : 18.39 Min. : 18.39
## 1st Qu.: 36.74 1st Qu.: 38.29 1st Qu.: 36.23 1st Qu.: 37.12
## Median : 116.77 Median : 114.53 Median : 120.01 Median : 114.88
## Mean : 852.05 Mean : 832.91 Mean : 860.43 Mean : 813.83
## 3rd Qu.: 591.70 3rd Qu.: 557.02 3rd Qu.: 595.01 3rd Qu.: 546.93
## Max. :38127.55 Max. :38852.22 Max. :38696.40 Max. :38072.58
## GSM240619 GSM240620 GSM240621 GSM240622
## Min. : 17.91 Min. : 16.31 Min. : 17.99 Min. : 17.46
## 1st Qu.: 37.29 1st Qu.: 36.73 1st Qu.: 36.51 1st Qu.: 36.27
## Median : 118.64 Median : 120.60 Median : 121.63 Median : 120.50
## Mean : 846.11 Mean : 852.04 Mean : 859.75 Mean : 858.06
## 3rd Qu.: 573.04 3rd Qu.: 585.35 3rd Qu.: 595.63 3rd Qu.: 600.90
## Max. :38327.45 Max. :38191.13 Max. :38500.64 Max. :37954.34
## GSM240623 GSM240624 GSM240625 GSM240626
## Min. : 19.14 Min. : 19.76 Min. : 18.23 Min. : 17.68
## 1st Qu.: 37.06 1st Qu.: 37.31 1st Qu.: 36.51 1st Qu.: 36.38
## Median : 117.51 Median : 117.70 Median : 124.53 Median : 120.28
## Mean : 842.21 Mean : 838.12 Mean : 848.70 Mean : 868.37
## 3rd Qu.: 583.56 3rd Qu.: 580.43 3rd Qu.: 606.83 3rd Qu.: 583.93
## Max. :39471.64 Max. :39395.66 Max. :36587.34 Max. :37534.33
## GSM240627 GSM240628 GSM240629
## Min. : 18.12 Min. : 16.88 Min. : 18.08
## 1st Qu.: 36.13 1st Qu.: 36.86 1st Qu.: 37.30
## Median : 125.42 Median : 117.23 Median : 103.54
## Mean : 854.13 Mean : 857.79 Mean : 866.89
## 3rd Qu.: 597.14 3rd Qu.: 583.94 3rd Qu.: 600.84
## Max. :38501.54 Max. :37813.12 Max. :39960.03boxplot(exprs(gse),outline=FALSE)
Normalize
exprs(gse)<-normalizeQuantiles(exprs(gse))
boxplot(exprs(gse),outline=FALSE)
## choose phenotype for pheatmap
stage<-select(sampleInfo,Stage)pheatmap for corelation matrix
corMatrix<-cor(exprs(gse),use="c")
pheatmap(corMatrix,annotation_col=stage,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=Stage)) +geom_point()
## Create classification for comparation
sampleInfo<-sampleInfo %>% mutate(group=recode(Stage,"Stage IA"="Localized", "Stage IB"="Localized", "Stage IIB"="Advanced", "Stage III"="Advanced"))Name group for design
design<-model.matrix(~0+sampleInfo$group)
design## sampleInfo$groupAdvanced sampleInfo$groupLocalized
## 1 0 1
## 2 1 0
## 3 0 1
## 4 0 1
## 5 1 0
## 6 0 1
## 7 0 1
## 8 0 1
## 9 0 1
## 10 1 0
## 11 0 1
## 12 0 1
## 13 0 1
## 14 1 0
## 15 1 0
## 16 0 1
## 17 0 1
## 18 1 0
## 19 0 1
## 20 1 0
## 21 0 1
## 22 1 0
## 23 0 1
## 24 0 1
## 25 0 1
## 26 0 1
## 27 0 1
## 28 0 1
## 29 0 1
## 30 0 1
## 31 1 0
## 32 1 0
## 33 0 1
## 34 0 1
## 35 1 0
## 36 0 1
## 37 0 1
## 38 0 1
## 39 0 1
## 40 0 1
## 41 0 1
## 42 1 0
## 43 1 0
## 44 0 1
## 45 0 1
## 46 1 0
## 47 0 1
## 48 1 0
## 49 1 0
## 50 0 1
## 51 0 1
## 52 1 0
## 53 0 1
## 54 1 0
## 55 0 1
## 56 0 1
## 57 0 1
## 58 0 1
## 59 0 1
## 60 1 0
## 61 0 1
## 62 0 1
## 63 1 0
## attr(,"assign")
## [1] 1 1
## attr(,"contrasts")
## attr(,"contrasts")$`sampleInfo$group`
## [1] "contr.treatment"Differential Expression Analasis
colnames(design)<-c("Advanced","Localized")
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(Advanced-Localized,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="GSE9479_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,"")) %>% ggplot(aes(x=logFC, y = B, col=Significant, label = Label)) + geom_point() + geom_text_repel(col="black")## Warning: ggrepel: 12 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
filter(full_results, adj.P.Val < 0.05, abs(logFC)>1) %>% write_csv(path="GSE9479_filtered_results.csv")