这是生信技能树关于生信人的20个R语言习题的答案：
1 安装R包
数据包： ALL, CLL, pasilla, airway
软件包：limma，DESeq2，clusterProfiler
工具包：reshape2
绘图包：ggplot2
#################################
1.安装R包
#################################
source(“https://bioconductor.org/biocLite.R”)
options(BioC_mirror=“http://mirrors.ustc.edu.cn/bioc/”)
biocLite(c(“ALL”,“CLL”, “pasilla”, “airway”)) #数据包
biocLite(c(“limma”,“DESeq2”, “clusterProfiler”)) #软件包
install.packages(“reshape2”) #工具包
install.packages(“ggplot2”) #绘图包
#另外可以检测某个包是否存在，只有不存在时才会安装
if (! require (‘CLL’)){
options(BioC_mirror=“http://mirrors.ustc.edu.cn/bioc/”)
BiocInstaller::biocLite(‘CLL’,ask = F, suppressUpdates = T)
}
#意思就是设定安装镜像，然后隐藏安装信息，中间过程不询问
biocLite(“hgu95av2.db”)
library(hgu95av2.db)
ls(“package:hgu95av2.db”)#产生36个映射数据（探针id转为36种主流id）
capture.output(hgu95av2()) #将hgu95av2包含的具体信息输出为字符串【只用hgu95av2()可能输出，但不易整理】
#结果可以看到，每一个子集都是有keys-value构成的，就相当于list的结构。因此对于hgu95av2的操作都要使用as.list()
as.list(hgu95av2SYMBOL[1]) #变成列表进行查看第一组元素
#结果得到的就是
#$1000_at
#[1] “MAPK3”
#不懂hgu95av2SYMBOL是什么意思，就查找
?hgu95av2SYMBOL
#结果得到：hgu95av2SYMBOL is an R object that provides mappings between manufacturer identifiers and gene abbreviations【就是基因id与探针号的关系，这里显示的就是MAPK3这个基因对应的探针id是1000_at】
ids=toTable(hgu95av2SYMBOL)
length(unique(ids
s
y
m
b
o
l
)
)
t
a
i
l
(
s
o
r
t
(
t
a
b
l
e
(
i
d
s
symbol)) tail(sort(table(ids
symbol))tail(sort(table(idssymbol)))
table(sort(table(ids
s
y
m
b
o
l
)
)
)
p
l
o
t
(
t
a
b
l
e
(
s
o
r
t
(
t
a
b
l
e
(
i
d
s
symbol))) plot(table(sort(table(ids
symbol)))plot(table(sort(table(idssymbol))))
table(rownames(exprSet) %in% ids
p
r
o
b
e
i
d
)
e
x
p
r
S
e
t
=
e
x
p
r
S
e
t
[
r
o
w
n
a
m
e
s
(
e
x
p
r
S
e
t
)
probe_id) exprSet=exprSet[rownames(exprSet) %in% ids
probeid)exprSet=exprSet[rownames(exprSet)probe_id,]
dim(exprSet)
ids=ids[match(rownames(exprSet),ids
p
r
o
b
e
i
d
)
,
]
h
e
a
d
(
i
d
s
)
e
x
p
r
S
e
t
[
1
:
5
,
1
:
5
]
t
m
p
=
b
y
(
e
x
p
r
S
e
t
,
i
d
s
probe_id),] head(ids) exprSet[1:5,1:5] tmp = by(exprSet,ids
probeid),]head(ids)exprSet[1:5,1:5]tmp=by(exprSet,idssymbol,function(x) rownames(x)[which.max(rowMeans(x))] )
probes = as.character(tmp)
exprSet=exprSet[rownames(exprSet) %in% probes ,]
dim(exprSet)
rownames(exprSet)=ids[match(rownames(exprSet),ids
p
r
o
b
e
i
d
)
,
2
]
e
x
p
r
S
e
t
[
1
:
5
,
1
:
5
]
l
i
b
r
a
r
y
(
r
e
s
h
a
p
e
2
)
e
x
p
r
S
e
t
L
=
m
e
l
t
(
e
x
p
r
S
e
t
)
c
o
l
n
a
m
e
s
(
e
x
p
r
S
e
t
L
)
=
c
(
′
p
r
o
b
e
′
,
′
s
a
m
p
l
e
′
,
′
v
a
l
u
e
′
)
e
x
p
r
S
e
t
L
probe_id),2] exprSet[1:5,1:5] library(reshape2) exprSet_L=melt(exprSet) colnames(exprSet_L)=c('probe','sample','value') exprSet_L
probeid),2]exprSet[1:5,1:5]library(reshape2)exprSetL=melt(exprSet)colnames(exprSetL)=c(′probe′,′sample′,′value′)exprSetLgroup=rep(group_list,each=nrow(exprSet))
head(exprSet_L)
ggplot2
library(ggplot2)
p=ggplot(exprSet_L,aes(x=sample,y=value,fill=group))+geom_boxplot()
print§
p=ggplot(exprSet_L,aes(x=sample,y=value,fill=group))+geom_violin()
print§
p=ggplot(exprSet_L,aes(value,fill=group))+geom_histogram(bins = 200)+facet_wrap(~sample, nrow = 4)
print§
p=ggplot(exprSet_L,aes(value,col=group))+geom_density()+facet_wrap(~sample, nrow = 4)
print§
p=ggplot(exprSet_L,aes(value,col=group))+geom_density()
print§
p=ggplot(exprSet_L,aes(x=sample,y=value,fill=group))+geom_boxplot()
p=p+stat_summary(fun.y=“mean”,geom=“point”,shape=23,size=3,fill=“red”)
p=p+theme_set(theme_set(theme_bw(base_size=20)))
p=p+theme(text=element_text(face=‘bold’),axis.text.x=element_text(angle=30,hjust=1),axis.title=element_blank())
print§
mean,median,max,min,sd,var,mad
g_mean <- tail(sort(apply(exprSet,1,mean)),50)
g_median <- tail(sort(apply(exprSet,1,median)),50)
g_max <- tail(sort(apply(exprSet,1,max)),50)
g_min <- tail(sort(apply(exprSet,1,min)),50)
g_sd <- tail(sort(apply(exprSet,1,sd)),50)
g_var <- tail(sort(apply(exprSet,1,var)),50)
g_mad <- tail(sort(apply(exprSet,1,mad)),50)
g_mad
names(g_mad)
heatmap
library(pheatmap)
choose_gene=names(tail(sort(apply(exprSet,1,mad)),50))
choose_matrix=exprSet[choose_gene,]
choose_matrix=t(scale(t(choose_matrix)))
pheatmap(choose_matrix)
UpSetR
https://cran.r-project.org/web/packages/UpSetR/README.html
library(UpSetR)
g_all <- unique(c(names(g_mean),names(g_median),names(g_max),names(g_min),
names(g_sd),names(g_var),names(g_mad) ))
dat=data.frame(g_all=g_all,
g_mean=ifelse(g_all %in% names(g_mean) ,1,0),
g_median=ifelse(g_all %in% names(g_median) ,1,0),
g_max=ifelse(g_all %in% names(g_max) ,1,0),
g_min=ifelse(g_all %in% names(g_min) ,1,0),
g_sd=ifelse(g_all %in% names(g_sd) ,1,0),
g_var=ifelse(g_all %in% names(g_var) ,1,0),
g_mad=ifelse(g_all %in% names(g_mad) ,1,0)
)
upset(dat,nsets = 7)
pdata=pData(sCLLex)
group_list=as.character(pdata[,2])
group_list
dim(exprSet)
exprSet[1:5,1:5]
hclust
colnames(exprSet)=paste(group_list,1:22,sep=’’)
Define nodePar
nodePar <- list(lab.cex = 0.6, pch = c(NA, 19),
cex = 0.7, col = “blue”)
hc=hclust(dist(t(exprSet)))
par(mar=c(5,5,5,10))
plot(as.dendrogram(hc), nodePar = nodePar, horiz = TRUE)
PCA
library(ggfortify)
df=as.data.frame(t(exprSet))
df$group=group_list
autoplot(prcomp( df[,1:(ncol(df)-1)] ), data=df,colour = ‘group’)
t.test
dat = exprSet
group_list=as.factor(group_list)
group1 = which(group_list == levels(group_list)[1])
group2 = which(group_list == levels(group_list)[2])
dat1 = dat[, group1]
dat2 = dat[, group2]
dat = cbind(dat1, dat2)
pvals = apply(exprSet, 1, function(x){
t.test(as.numeric(x)~group_list)$p.value
})
p.adj = p.adjust(pvals, method = “BH”)
avg_1 = rowMeans(dat1)
avg_2 = rowMeans(dat2)
log2FC = avg_2-avg_1
DEG_t.test = cbind(avg_1, avg_2, log2FC, pvals, p.adj)
DEG_t.test=DEG_t.test[order(DEG_t.test[,4]),]
DEG_t.test=as.data.frame(DEG_t.test)
head(DEG_t.test)
DEG by limma
suppressMessages(library(limma))
design <- model.matrix(~0+factor(group_list))
colnames(design)=levels(factor(group_list))
rownames(design)=colnames(exprSet)
design
contrast.matrix<-makeContrasts(paste0(unique(group_list),collapse = “-”),levels = design)
contrast.matrix
##这个矩阵声明，我们要把progres.组跟stable进行差异分析比较
##step1
fit <- lmFit(exprSet,design)
##step2
fit2 <- contrasts.fit(fit, contrast.matrix) ##这一步很重要，大家可以自行看看效果
fit2 <- eBayes(fit2) ## default no trend !!!
##eBayes() with trend=TRUE
##step3
tempOutput = topTable(fit2, coef=1, n=Inf)
nrDEG = na.omit(tempOutput)
#write.csv(nrDEG2,“limma_notrend.results.csv”,quote = F)
head(nrDEG)
volcano plot
DEG=nrDEG
logFC_cutoff <- with(DEG,mean(abs( logFC)) + 2*sd(abs( logFC)) )
DEG
c
h
a
n
g
e
=
a
s
.
f
a
c
t
o
r
(
i
f
e
l
s
e
(
D
E
G
change = as.factor(ifelse(DEG
change=as.factor(ifelse(DEGP.Value < 0.05 & abs(DEG
l
o
g
F
C
)
>
l
o
g
F
C
c
u
t
o
f
f
,
i
f
e
l
s
e
(
D
E
G
logFC) > logFC_cutoff, ifelse(DEG
logFC)>logFCcutoff,ifelse(DEGlogFC > logFC_cutoff ,‘UP’,‘DOWN’),‘NOT’)
)
this_tile <- paste0('Cutoff for logFC is ',round(logFC_cutoff,3),
'\nThe number of up gene is ',nrow(DEG[DEGKaTeX parse error: Expected 'EOF', got '\nThe' at position 21: …e =='UP',]) , '\̲n̲T̲h̲e̲ ̲number of down …change ==‘DOWN’,])
)
g = ggplot(data=DEG, aes(x=logFC, y=-log10(P.Value), color=change)) +
geom_point(alpha=0.4, size=1.75) +
theme_set(theme_set(theme_bw(base_size=20)))+
xlab(“log2 fold change”) + ylab("-log10 p-value") +
ggtitle( this_tile ) + theme(plot.title = element_text(size=15,hjust = 0.5))+
scale_colour_manual(values = c(‘blue’,‘black’,‘red’)) ## corresponding to the levels(res$change)
print(g)
different P values
head(nrDEG)
head(DEG_t.test)
DEG_t.test=DEG_t.test[rownames(nrDEG),]
plot(DEG_t.test[,3],nrDEG[,1])
plot(DEG_t.test[,4],nrDEG[,4])
plot(-log10(DEG_t.test[,4]),-log10(nrDEG[,4]))
exprSet[‘GAPDH’,]
exprSet[‘ACTB’,]
exprSet[‘DLEU1’,]
library(ggplot2)
library(ggpubr)
my_comparisons <- list(
c(“stable”, “progres.”)
)
dat=data.frame(group=group_list,
sampleID= names(exprSet[‘DLEU1’,]),
values= as.numeric(exprSet[‘DLEU1’,]))
ggboxplot(
dat, x = “group”, y = “values”,
color = “group”,
add = “jitter”
)+
stat_compare_means(comparisons = my_comparisons, method = “t.test”)
heatmap
library(pheatmap)
choose_gene=head(rownames(nrDEG),25)
choose_matrix=exprSet[choose_gene,]
choose_matrix=t(scale(t(choose_matrix)))
pheatmap(choose_matrix)
生信技能树
扫一扫
2526
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
