箱线图走一个
2019-02-11 本文已影响5人
小洁忘了怎么分身
# 从内置数据集的表达矩阵中找TP53基因的表达量
rm(list=ls())
options(stringsAsFactors = F)
if(!suppressMessages(require("hgu95av2.db")))BiocManager::install("hgu95av2.db")
if(!suppressMessages(require("CLL")))BiocManager::install("CLL")
if(!suppressMessages(require(tidyverse)))install.packages("tidyverse")
if(!suppressMessages(require(ggpubr)))install.packages("ggpubr")
suppressMessages(library(hgu95av2.db))
suppressMessages(library(CLL))
suppressMessages(library(ggpubr))
data(sCLLex)
# sCLLex
exprSet <- exprs(sCLLex) #探针的表达量
exprSet[1:4,1:4]
#> CLL11.CEL CLL12.CEL CLL13.CEL CLL14.CEL
#> 1000_at 5.743132 6.219412 5.523328 5.340477
#> 1001_at 2.285143 2.291229 2.287986 2.295313
#> 1002_f_at 3.309294 3.318466 3.354423 3.327130
#> 1003_s_at 1.085264 1.117288 1.084010 1.103217
pd <- pData(sCLLex) #sampleID与disease的对应关系
head(pd)
#> SampleID Disease
#> CLL11.CEL CLL11 progres.
#> CLL12.CEL CLL12 stable
#> CLL13.CEL CLL13 progres.
#> CLL14.CEL CLL14 progres.
#> CLL15.CEL CLL15 progres.
#> CLL16.CEL CLL16 progres.
p2s <- toTable(hgu95av2SYMBOL) #探针与symbol的对应关系
head(p2s)
#> probe_id symbol
#> 1 1000_at MAPK3
#> 2 1001_at TIE1
#> 3 1002_f_at CYP2C19
#> 4 1003_s_at CXCR5
#> 5 1004_at CXCR5
#> 6 1005_at DUSP1
p3 <- filter(p2s,symbol=='TP53')
# boxplot [find TP53 has 3 probe IDs]
probe_tp53 <- p3$probe_id
for(i in 1:3){
boxplot(exprSet[probe_tp53[i],] ~ pd$Disease)
}



#用ggpubr作图
#http://www.sthda.com/english/articles/24-ggpubr-publication-ready-p#lots/
expd <- rownames_to_column(as.data.frame(exprSet))
expd[1:4,1:4]
#> rowname CLL11.CEL CLL12.CEL CLL13.CEL
#> 1 1000_at 5.743132 6.219412 5.523328
#> 2 1001_at 2.285143 2.291229 2.287986
#> 3 1002_f_at 3.309294 3.318466 3.354423
#> 4 1003_s_at 1.085264 1.117288 1.084010
expd2 <- gather(expd,
key = 'sample',
value = 'exp',-1)
pd <- rownames_to_column(pd)
head(pd)
#> rowname SampleID Disease
#> 1 CLL11.CEL CLL11 progres.
#> 2 CLL12.CEL CLL12 stable
#> 3 CLL13.CEL CLL13 progres.
#> 4 CLL14.CEL CLL14 progres.
#> 5 CLL15.CEL CLL15 progres.
#> 6 CLL16.CEL CLL16 progres.
expd3 <- inner_join(expd2,pd,by=c('sample'='rowname'))
i=1 ###可换1,2
for(i in 1:3){
p <- ggboxplot(filter(expd3,rowname==probe_tp53[i]),
x = 'Disease',
y = 'exp',
color = "Disease", palette =c("#00AFBB", "#E7B800", "#FC4E07"),
add = "jitter", shape = "Disease")
print(p)
}


