library(readxl)
library(ggplot2)
library(plotly)
##
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
library(lattice)
library(colorspace)
library(readr)
library(SAM)
## Loading required package: splines
library(splines)
library(splines2)
library(siggenes)
## 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, 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")'.
## Loading required package: multtest
library(Biobase)
library(BiocGenerics)
library(parallel)
library(multtest)
#library(MESS)
#library(cp4p)
library(qvalue)
library(EnhancedVolcano)
## Loading required package: ggrepel
#library(wordcloud)
#library(wordcloud2)
library(devtools)
## Loading required package: usethis
library(rlang)
##
## Attaching package: 'rlang'
## The following object is masked from 'package:Biobase':
##
## exprs
library(basicPlotteR)
library(ggrepel)
library(qvalue)
#library(fdrci)
library(readr)
library(readxl)
library(tidyverse)
## ── Attaching packages
## ───────────────────────────────────────
## tidyverse 1.3.2 ──
## ✔ tibble 3.1.8 ✔ dplyr 1.0.10
## ✔ tidyr 1.2.1 ✔ stringr 1.5.0
## ✔ purrr 1.0.0 ✔ forcats 0.5.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ purrr::%@%() masks rlang::%@%()
## ✖ dplyr::combine() masks Biobase::combine(), BiocGenerics::combine()
## ✖ rlang::exprs() masks Biobase::exprs()
## ✖ dplyr::filter() masks plotly::filter(), stats::filter()
## ✖ purrr::flatten() masks rlang::flatten()
## ✖ purrr::flatten_chr() masks rlang::flatten_chr()
## ✖ purrr::flatten_dbl() masks rlang::flatten_dbl()
## ✖ purrr::flatten_int() masks rlang::flatten_int()
## ✖ purrr::flatten_lgl() masks rlang::flatten_lgl()
## ✖ purrr::flatten_raw() masks rlang::flatten_raw()
## ✖ purrr::invoke() masks rlang::invoke()
## ✖ dplyr::lag() masks stats::lag()
## ✖ BiocGenerics::Position() masks ggplot2::Position(), base::Position()
## ✖ purrr::splice() masks rlang::splice()
install.packages("devtools")
## Warning: package 'devtools' is in use and will not be installed
devtools::install_github("JosephCrispell/basicPlotteR")
## WARNING: Rtools is required to build R packages, but is not currently installed.
##
## Please download and install Rtools 4.2 from https://cran.r-project.org/bin/windows/Rtools/ or https://www.r-project.org/nosvn/winutf8/ucrt3/.
## Skipping install of 'basicPlotteR' from a github remote, the SHA1 (4628140c) has not changed since last install.
## Use `force = TRUE` to force installation
if (!requireNamespace('BiocManager', quietly = TRUE))
install.packages('BiocManager')
BiocManager::install('EnhancedVolcano')
## Bioconductor version 3.15 (BiocManager 1.30.19), R 4.2.2 (2022-10-31 ucrt)
## Warning: package(s) not installed when version(s) same as or greater than current; use
## `force = TRUE` to re-install: 'EnhancedVolcano'
## Installation paths not writeable, unable to update packages
## path: C:/Program Files/R/R-4.2.2/library
## packages:
## boot, foreign, survival
## Old packages: 'BiocParallel', 'Biostrings', 'cli', 'data.table', 'IRanges',
## 'Matrix', 'nlme', 'plyr', 'png', 'progressr', 'purrr', 'rbibutils', 'Rcpp',
## 'RCurl'
BiocManager::install("siggenes")
## Bioconductor version 3.15 (BiocManager 1.30.19), R 4.2.2 (2022-10-31 ucrt)
## Warning: package(s) not installed when version(s) same as or greater than current; use
## `force = TRUE` to re-install: 'siggenes'
## Installation paths not writeable, unable to update packages
## path: C:/Program Files/R/R-4.2.2/library
## packages:
## boot, foreign, survival
## Old packages: 'BiocParallel', 'Biostrings', 'cli', 'data.table', 'IRanges',
## 'Matrix', 'nlme', 'plyr', 'png', 'progressr', 'purrr', 'rbibutils', 'Rcpp',
## 'RCurl'
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("qvalue")
## Bioconductor version 3.15 (BiocManager 1.30.19), R 4.2.2 (2022-10-31 ucrt)
## Warning: package(s) not installed when version(s) same as or greater than current; use
## `force = TRUE` to re-install: 'qvalue'
## Installation paths not writeable, unable to update packages
## path: C:/Program Files/R/R-4.2.2/library
## packages:
## boot, foreign, survival
## Old packages: 'BiocParallel', 'Biostrings', 'cli', 'data.table', 'IRanges',
## 'Matrix', 'nlme', 'plyr', 'png', 'progressr', 'purrr', 'rbibutils', 'Rcpp',
## 'RCurl'
data_table <- read_delim("Z:/TSever/Orbitrap/20200824_lgmn_kidney_WT_44_45_46_KO_58_60_68_v2p32/MaxQuant/combined/txt/additional tables/data table.txt",
"\t", escape_double = FALSE, trim_ws = TRUE)
## New names:
## Rows: 1129 Columns: 9
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "\t" chr
## (1): ...1 dbl (8): WT44, WT45, WT46, KO58, KO60, KO68, difference, p-value
## ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
## Specify the column types or set `show_col_types = FALSE` to quiet this message.
## • `` -> `...1`
#View(data_table)
dt3 <- read_excel("Z:/TSever/Orbitrap/20200824_lgmn_kidney_WT_44_45_46_KO_58_60_68_v2p32/MaxQuant/combined/txt/additional tables/data table3.xlsx",
sheet = "data_table1 (3)")
## New names:
## • `` -> `...18`
## • `` -> `...20`
## • `` -> `...21`
## • `` -> `...22`
## • `` -> `...23`
## • `` -> `...24`
## • `` -> `...25`
## • `` -> `...26`
names(data_table)[names(data_table) == "p-value"] <- 'p'
labels<- read_excel("Z:/TSever/Orbitrap/20200824_lgmn_kidney_WT_44_45_46_KO_58_60_68_v2p32/MaxQuant/combined/txt/additional tables/data table3.xlsx",
sheet = "data1")
## New names:
## • `` -> `...18`
labels<- labels[ !(labels$Name %in% c('TAU','VINC','ACACB')), ]
q0<- read_excel("Z:/TSever/Orbitrap/20200824_lgmn_kidney_WT_44_45_46_KO_58_60_68_v2p32/MaxQuant/combined/txt/additional tables/data table3.xlsx",
sheet = "data2")
ggplot(data_table, aes(x=difference , y=p))+
geom_point(size=0.7, shape=19)
r<- dt3$`Test statistic`
s<-dt3$rank
fudge2(r, s, alpha = seq(0, 1, 0.05), include.zero = T)
## $alpha.hat
## [1] 1
##
## $s.zero
## 100%
## 1129
##
## $vec.cv
## [1] 3.5339473 3.2810390 1.1598744 0.9172485 0.7961094 0.7201776 0.6673186
## [8] 0.6282224 0.5981610 0.5744053 0.5552447 0.5395404 0.5265040 0.5155669
## [15] 0.5063088 0.4984116 0.4916301 0.4857725 0.4806867 0.4762504 0.4723646
## [22] 0.4689482
##
## $msg
## [1] "s0 = 1129 (The 100 % quantile of the s values.) \n \n"
# Function to compute the smooth threshold curve
smooth.threshold=function(x,ta,s0,df){
xp=x[x>(ta*s0)]; xn=x[x<(-ta*s0)]; dp=xp/ta-s0; dn=xn/(-ta)-s0;
dp=s0/dp; dp=ta*(1+dp); dn=s0/dn; dn=ta*(1+dn);
fp=pt(dp,df=df); fn=pt(dn,df=df); yp=-log10(2*(1-fp)); yn=-log10(2*(1-fn));
return(cbind(c(xn,xp),c(yn,yp)));
}
ff <- fudge2(dt3$difference, dt3$sd, include.zero = F)
ff
## $alpha.hat
## [1] 0
##
## $s.zero
## 0%
## 0.08827611
##
## $vec.cv
## [1] 0.3879591 0.4309007 0.4469155 0.4595445 0.4673658 0.4760025 0.4842464
## [8] 0.4929749 0.5008420 0.5103027 0.5174509 0.5266599 0.5362250 0.5445787
## [15] 0.5561395 0.5656553 0.5743912 0.5878559 0.6014645 0.6183764 0.6687815
##
## $msg
## [1] "s0 = 0.0883 (The 0 % quantile of the s values.) \n \n"
df<- (3+3-2)
confidence_level<- 0.05
smoothcurve =smooth.threshold(seq(3,3,by=0.0001),ta=qt(confidence_level,df=df),0.0883,df=df)
ggplot(dt3, aes(x=difference , y=p))+
geom_jitter()
#geom_point(size=0.7, shape=19)
#geom_point(aes(difference[(Protein)],p[(Protein)]))+
#smoothcurve =smooth.threshold(seq(3,3,by=0.0001),ta=qt(confidence_level,df=df),0.0883,df=df)+
#geom_line(aes(smoothcurve, col='green',lwd=2))
abs.inf=-120
abs.sup=120
c<-(-1:1)
c1<-as.vector(labels$difference)
c2<-(labels$p)
x=dt3$difference
y=dt3$p
plot(x=dt3$difference, y=dt3$p, xlab="Difference wt-ko", ylab="-log p", col=ifelse(x< -1 & y>2.086 | x>1 & y>2.086, "red", "black"), pch=ifelse(x< -1 & y>2.086 | x>1 & y>2.086, 19,20), xlim=c(-5,5))
#xlim=c(abs.inf, abs.sup)
#pnt <- identify(x, y, plot = F)
#points(x[pnt], y[pnt], col = "red")
confidence_level=0.990; ff=0.0883;
smoothcurve=smooth.threshold(seq(abs.inf,abs.sup,by=0.0001),ta=qt(confidence_level,df=df),s0=ff,df=df);
lines(smoothcurve, col='#b21111',lwd=2)
#abline(v=c(-1,1), h=-log10(0.01))
#text(labels$difference,labels$p, labels=labels$NameAAposition, cex=0.7, font=2, pos=4, offset = 0.5, col='#000000', family='sans')
addTextLabels(labels$difference, labels$p, labels$NameAAposition, cex.label=0.7, col.label='black', col.background=NULL, avoidPoints=TRUE)
points(q0$difference, q0$p, col='green', type='p')
#par(new=TRUE)
#plot(q0$difference, q0$p, type="p", col="green" )
plot(x=dt3$difference, y=dt3$p, xlab="Difference wt-ko", ylab="-log p", col=ifelse(x< -1 & y>2.0 | x>1 & y>2.0, "red", "black"), pch=ifelse(x< -1 & y>2.0 | x>1 & y>2.0, 19,20), xlim=c(-5,5), xaxt='n')
#xlim=c(abs.inf, abs.sup)
#pnt <- identify(x, y, plot = F)
#points(x[pnt], y[pnt], col = "red")
#confidence_level=0.990; ff=0.0883;
#smoothcurve=smooth.threshold(seq(abs.inf,abs.sup,by=0.0001),ta=qt(confidence_level,df=df),s0=ff,df=df);
#lines(smoothcurve, col='#b21111',lwd=2)
abline(v=c(-1,1), h=-log10(0.01))
#text(labels$difference,labels$p, labels=labels$NameAAposition, cex=0.7, font=2, pos=4, offset = 0.5, col='#000000', family='sans')
addTextLabels(labels$difference, labels$p, labels$NameAAposition, cex.label=0.7, col.label='black', col.background=NULL, avoidPoints=FALSE)
points(q0$difference, q0$p, col='green', type='p', pch=8)
#par(new=TRUE)
#plot(q0$difference, q0$p, type="p", col="green" )
axis(1, at = seq(-5, 5, by = 1), las=1)
plot(x=dt3$difference, y=dt3$p, xlab="Difference wt-ko", ylab="-log p", col=ifelse(x< -1 & y>2.0 | x>1 & y>2.0, "red", "black"), pch=ifelse(x< -1 & y>2.0 | x>1 & y>2.0, 19,20), xlim=c(-5,5), xaxt='n')
#xlim=c(abs.inf, abs.sup)
#pnt <- identify(x, y, plot = F)
#points(x[pnt], y[pnt], col = "red")
#confidence_level=0.990; ff=0.0883;
#smoothcurve=smooth.threshold(seq(abs.inf,abs.sup,by=0.0001),ta=qt(confidence_level,df=df),s0=ff,df=df);
#lines(smoothcurve, col='#b21111',lwd=2)
abline(v=c(-1,1), h=-log10(0.01))
#text(labels$difference,labels$p, labels=labels$NameAAposition, cex=0.7, font=2, pos=4, offset = 0.5, col='#000000', family='sans')
addTextLabels(labels$difference, labels$p, labels$NameAAposition, cex.label=0.99, lwd=1,col.label='black', col.background=NULL, avoidPoints=T)
# points(q0$difference, q0$p, col='green', type='p', pch=8)
#par(new=TRUE)
#plot(q0$difference, q0$p, type="p", col="green" )
axis(1, at = seq(-5, 5, by = 1), las=1)
ggplot(dt3, aes(x=difference, y=p))+
geom_point()+
geom_label(
data=dt3 %>% filter(p>2 & (difference>1 | difference<1)),
aes(label=NameAAposition),
nudge_x = 0.25, nudge_y = 0.25,
label.padding = unit(0.001, "lines"), # Rectangle size around label
label.size = 0.001,
label.r = unit(0.1,'lines'),
color = "black",
fill="#bada55",
check_overlap=TRUE
)
## Warning in geom_label(data = dt3 %>% filter(p > 2 & (difference > 1 | difference
## < : Ignoring unknown parameters: `check_overlap`
ggplot(dt3, aes(x=difference, y=p))+
geom_point()+
geom_text(
data=dt3 %>% filter(p>2 & (difference>1 | difference<1)),
aes(label=NameAAposition),
nudge_x = 0.25, nudge_y = 0.25,
label.padding = unit(0.001, "lines"), # Rectangle size around label
label.size = 0.001,
label.r = unit(0.1,'lines'),
color = "black",
fill="#bada55",
check_overlap=FALSE
)+
theme_classic()
## Warning in geom_text(data = dt3 %>% filter(p > 2 & (difference > 1 | difference < : Ignoring unknown parameters: `label.padding`, `label.size`, `label.r`, and
## `fill`
abs.inf=-6
abs.sup=6
plot(x=dt3$difference, y=dt3$p, xlab="difference wt-ko", ylab="-log p")
xlim=c(abs.inf, abs.sup)
confidence_level=0.99; ff=0.1;
smoothcurve=smooth.threshold(seq(abs.inf,abs.sup,by=0.0001),ta=qt(confidence_level,df=df),s0=ff,df=df);
lines(smoothcurve, col='pink',lwd=2)
abline(v=c(-1,1), h=-log10(0.01))
text(dt3$difference,dt3$p, labels=dt3$Name, cex=0.7, font=2, pos=4, col='#bada55')
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("qvalue")
## Bioconductor version 3.15 (BiocManager 1.30.19), R 4.2.2 (2022-10-31 ucrt)
## Warning: package(s) not installed when version(s) same as or greater than current; use
## `force = TRUE` to re-install: 'qvalue'
## Installation paths not writeable, unable to update packages
## path: C:/Program Files/R/R-4.2.2/library
## packages:
## boot, foreign, survival
## Old packages: 'BiocParallel', 'Biostrings', 'cli', 'data.table', 'IRanges',
## 'Matrix', 'nlme', 'plyr', 'png', 'progressr', 'purrr', 'rbibutils', 'Rcpp',
## 'RCurl'
options(repr.plot.width=2.5, repr.plot.height=2.5)
hist(dt3$pp, nclass=20)
qobj<-qvalue(p=dt3$pp, fdr.level=0.01)
summary(qobj)
##
## Call:
## qvalue(p = dt3$pp, fdr.level = 0.01)
##
## pi0: 0.6591035
##
## Cumulative number of significant calls:
##
## <1e-04 <0.001 <0.01 <0.025 <0.05 <0.1 <1
## p-value 1 9 38 73 116 215 1129
## q-value 0 0 0 0 2 11 1129
## local FDR 0 0 0 0 2 8 1129
options(repr.plot.width=6, repr.plot.height=5)
plot(qobj)
if (!requireNamespace('BiocManager', quietly = TRUE))
install.packages('BiocManager')
BiocManager::install('EnhancedVolcano')
## Bioconductor version 3.15 (BiocManager 1.30.19), R 4.2.2 (2022-10-31 ucrt)
## Warning: package(s) not installed when version(s) same as or greater than current; use
## `force = TRUE` to re-install: 'EnhancedVolcano'
## Installation paths not writeable, unable to update packages
## path: C:/Program Files/R/R-4.2.2/library
## packages:
## boot, foreign, survival
## Old packages: 'BiocParallel', 'Biostrings', 'cli', 'data.table', 'IRanges',
## 'Matrix', 'nlme', 'plyr', 'png', 'progressr', 'purrr', 'rbibutils', 'Rcpp',
## 'RCurl'
library(EnhancedVolcano)
v1<-EnhancedVolcano(dt3, lab = dt3$NameAAposition, x= 'difference', y= 'pp', pCutoff = 10e-3, FCcutoff = 1,
ylim = c(0, 5),
colAlpha = 0.9,
title = 'Volcano plot of Legumain ko/wt murine kidney phosphoproteomic analysis',
legendPosition = 'right',
subtitle = 'Differential expression',
colGradient = c('red3', 'royalblue'),
drawConnectors = TRUE,
hline = c(10e-8),
widthConnectors = 0.5,
gridlines.major = FALSE,
gridlines.minor = FALSE)+
ggplot2::coord_cartesian(xlim=c(-5.5, 5.5)) +
ggplot2::scale_x_continuous(
breaks=seq(-5,5, 1))
## Coordinate system already present. Adding new coordinate system, which will
## replace the existing one.
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
v1
## Warning: Removed 1 rows containing missing values (`geom_hline()`).
EnhancedVolcano(dt3, lab = dt3$Name, x= 'difference', y= 'pp', pCutoff = 1*10e-3, FCcutoff = 1,
ylim = c(0, 5), xlab = 'difference (WT-KO)',
colAlpha = 0.9,
drawConnectors = FALSE,
gridlines.major = FALSE,
gridlines.minor = FALSE)
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("qvalue")
## Bioconductor version 3.15 (BiocManager 1.30.19), R 4.2.2 (2022-10-31 ucrt)
## Warning: package(s) not installed when version(s) same as or greater than current; use
## `force = TRUE` to re-install: 'qvalue'
## Installation paths not writeable, unable to update packages
## path: C:/Program Files/R/R-4.2.2/library
## packages:
## boot, foreign, survival
## Old packages: 'BiocParallel', 'Biostrings', 'cli', 'data.table', 'IRanges',
## 'Matrix', 'nlme', 'plyr', 'png', 'progressr', 'purrr', 'rbibutils', 'Rcpp',
## 'RCurl'
EnhancedVolcano(dt3, lab = dt3$Name, x= 'difference', y= 'pp', pCutoff = 1*10e-3, FCcutoff = 1,
ylim = c(0, 5), xlab = 'difference (WT-KO)',
colAlpha = 0.9,
drawConnectors = FALSE,
gridlines.major = FALSE,
gridlines.minor = FALSE)
