高级题链接：http://www.bio-info-trainee.com/3409.html
第一题
#1
#切换镜像
rm(list = ls())
options()$repos
options()$BioC_mirror
options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")
options("repos" = c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
#对于CRAN包，install.package
#对于bioconductor包，BiocManager::install("package_name",ask = F,update = F)
第二题
#2
suppressPackageStartupMessages(library(CLL))
data(sCLLex) #sCLLex是依赖于CLL这个package的一个对象
sCLLex
dat=exprs(sCLLex)
dim(dat)
dat[1:5,1:5] #提取表达矩阵并查看
samples=sampleNames(sCLLex) #提取样本信息
pdata=pData(sCLLex)
group_list=as.character(pdata[,2])
group_list #提取分组信息并查看
第三题
#3
str(dat)
head(dat)
第四题
#4
#BiocManager::install('hgu95av2.db')
suppressPackageStartupMessages(library(hgu95av2.db))
ls("package:hgu95av2.db")
第五题
#5
ids = toTable(hgu95av2SYMBOL) #用toTable函数提取探针和基因的对应关系
head(ids)
ids[ids$symbol=='TP53',] #找到TP53对应的探针
第六题
#6
table(ids$symbol) #table很关键
sort(table(ids$symbol))
sort(table(ids$symbol),decreasing = T)[1] #最多对应8个探针
plot(table(sort(table(ids$symbol)))) #直接画图可以看到整体情况
probe_num.png
第七题
#7
table(rownames(dat) %in% ids$probe_id) #可以看到有1165个探针没有对应的基因名
第八题
#8
dim(dat)
dat[1:4,1:4]
dat = dat[rownames(dat) %in% ids$probe_id,] #过滤表达矩阵
dat = dat[ids$probe_id,] #一样的效果
dim(dat)
dat[1:4,1:4] #处理前后检查一下表达矩阵，事半功倍
第九题
#9
#以下是我摸索match函数用法的过程
a = c('a','b','c')
b = c('b','c','a')
match(a,b)
match(b,a)
rm(a,b)
ids$mean = apply(dat,1,mean)
ids = ids[order(ids$symbol,ids$mean,decreasing = T),]
ids = ids[!duplicated(ids$symbol),]
dat = dat[ids$probe_id,] #常规的过滤步骤
第十题
#10
rownames(dat) = ids$symbol
第十一题
#11
dat[,1]
boxplot(dat[,1])
hist(dat[,1])
plot(density(dat[,1]))
library(reshape2)
dat_L=melt(dat) #宽矩阵变长矩阵
colnames(dat_L)=c('symbol','sample','value')
dat_L$group=rep(group_list,each=nrow(dat)) #特别注意加上each参数
head(dat_L)
library(ggplot2)
p=ggplot(dat_L,aes(x=sample,y=value,fill=group))+geom_boxplot() #ggplot的“映射”概念，按照group来涂色
print(p)
p=ggplot(dat_L,aes(x=sample,y=value,fill=group))+geom_violin()
print(p)
#这两张图横轴上字太多了，略丑，后面会调整
p=ggplot(dat_L,aes(value,fill=group))+geom_histogram(bins = 200)+facet_wrap(~sample, nrow = 4) #facet_wrap一张画布画多图
print(p)
p=ggplot(dat_L,aes(value,col=group))+geom_density()+facet_wrap(~sample, nrow = 4)
print(p)
p=ggplot(dat_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(p)
ggsave(filename = "test.png")
#高级箱线图，添加了表示均值的红点，改变了背景，调整了横轴标签
多图预警！
boxplot.png
histogram.png
density.png
ggplot-boxplot.png
ggplot-violin.png
上面这两张图横轴上字太多了，略丑，后面会调整
ggplot-hist.png
ggplot-density.png
ggplot-boxplot-advanced.png
第十二题
#12
# apply(dat,1,mean)
# apply(dat,1,max)
# apply(dat,1,min)
# apply(dat,1,median)
# apply(dat,1,sd)
# apply(dat,1,var)
# apply(dat,1,mad) #这些本质都是一样的
top = sort(apply(dat,1,mad),decreasing = T)[1:50]
#tail(sort(apply(exprSet,1,mad)),50) #一样的效果
str(top)
#题目中说这题不合规，我猜测可能是因为1.拿已筛选的表达矩阵来做这步处理以偏概全；2.得到的不是列表
第十三题
#13
dat2=dat[rownames(dat) %in% names(top),]
pheatmap::pheatmap(t(scale(t(dat2))))
pheatmap::pheatmap(dat2,scale = 'row') #两种画法本质一样
#暂时没有采用另外几个热图画法
heatmap.png
第十四题
#14
g_mean <- tail(sort(apply(dat,1,mean)),50)
g_median <- tail(sort(apply(dat,1,median)),50)
g_max <- tail(sort(apply(dat,1,max)),50)
g_min <- tail(sort(apply(dat,1,min)),50)
g_sd <- tail(sort(apply(dat,1,sd)),50)
g_var <- tail(sort(apply(dat,1,var)),50)
g_mad <- tail(sort(apply(dat,1,mad)),50)
#首先取这一堆统计量的top50基因
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)))
df=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(df,nsets = 7) #以后可以套用的常规流程
overlap.png
第十五题
#15
#之前的pdata就是表型数据，还从中取了一个group_list
第十六题
#16
#清理一下工作环境
save(dat,pdata,group_list,file = 'input.Rdata')
rm(list = ls())
options(stringsAsFactors = F)
load(file = 'input.Rdata')
colnames(dat) = paste(group_list,1:22,sep='') #参考答案的，想不出这步说明我还没有彻底适应R语言的“向量”概念
#从这步开始就是聚类分析常规步骤
nodePar <- list(lab.cex = 0.6, pch = c(NA, 19),
cex = 0.8, col = "red") #定义聚类图里的点长啥样
hc=hclust(dist(t(dat))) #用dist计算“距离”
par(mar=c(3,2,2,3)) #调整画板参数（页边距）
plot(as.dendrogram(hc), nodePar = nodePar, horiz = TRUE)
hclust.png
第十七题
#17
#清理一下工作环境
save(dat,pdata,group_list,file = 'input.Rdata')
rm(list = ls())
options(stringsAsFactors = F)
load(file = 'input.Rdata')
#PCA的基本步骤
#BiocManager::install('ggfortify')
library(ggfortify)
df=as.data.frame(t(dat))
df$group=group_list
autoplot(prcomp(df[,1:(ncol(df)-1)]), data = df, colour = 'group') #prcomp是进行主成分分析的函数
library("FactoMineR")
library("factoextra") #画主成分分析图需要加载这两个包
df=as.data.frame(t(dat))
dat.pca <- PCA(df, graph = FALSE) #之前prcomp的东西是带有分组信息的，现在重新赋值一个dat.pca，不带分组信息。
fviz_pca_ind(dat.pca,
geom.ind = "point", # show points only (nbut not "text")
col.ind = group_list, # color by groups
# palette = c("#00AFBB", "#E7B800"),
addEllipses = TRUE, # Concentration ellipses
legend.title = "Groups"
)
PCA1.png
PCA2.png
第十八题
#18
#清理一下工作环境
save(dat,pdata,group_list,file = 'input.Rdata')
rm(list = ls())
options(stringsAsFactors = F)
load(file = 'input.Rdata')
#t检验
pvals = apply(dat, 1, function(x){
t.test(as.numeric(x)~group_list)$p.value
}) #要尽量避免循环，多用apply
p.adj = p.adjust(pvals, method = "BH")
#以下都是参考答案的，这才知道虽然题目只说了“批量t检验”，但既然是差异分析就要有fold change
group1 = which(group_list == "progres.")
group2 = which(group_list == "stable") #获取分组索引
dat1 = dat[, group1]
dat2 = dat[, group2] #分离表达矩阵
# 答案中是化为因子后再操作的，其实可以直接做
# 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)
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]),] #按照原始p值排列
DEG_t.test=as.data.frame(DEG_t.test)
head(DEG_t.test)
第十九题
#19
#DEG by limma
#流程：design-contrast.matrix-step1/2/3
library(limma)
design <- model.matrix(~0 + factor(group_list))
colnames(design)=levels(factor(group_list))
rownames(design)=colnames(dat)
design
contrast.matrix <- makeContrasts(paste0(unique(group_list), collapse = "-"), levels = design)
contrast.matrix
#step1
fit <- lmFit(dat,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)
DEG = na.omit(tempOutput)
#write.csv(nrDEG2,"limma_notrend.results.csv",quote = F)
head(DEG)
#以下是画火山图的代码
logFC_cutoff <- with(DEG,mean(abs(logFC)) + 2*sd(abs(logFC))) #先找一个分界线
DEG$change = as.factor(ifelse(DEG$P.Value < 0.05 & abs(DEG$logFC) > logFC_cutoff,
ifelse(DEG$logFC > logFC_cutoff ,'UP','DOWN'),'NOT')
) #定义上调、下调和无变化
this_tile <- paste0('Cutoff for logFC is ',round(logFC_cutoff,3),
'\nThe number of up gene is ',nrow(DEG[DEG$change =='UP',]) ,
'\nThe number of down gene is ',nrow(DEG[DEG$change =='DOWN',])
) #编辑火山图的标题
this_tile
head(DEG)
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')) #颜色是根据DEG$change来涂的
print(g)
volcano.png
第二十题
#20
head(DEG)
head(DEG_t.test)
DEG_t.test=DEG_t.test[rownames(DEG),] #首先使两个数据框的行名顺序一致
plot(DEG_t.test[,3],DEG[,1]) #可以看到两种方法得到的fold change负相关，这是因为两次的谁和谁比不一样
plot(DEG_t.test[,4],DEG[,4]) #可以看到使用limma包和t.test的p值差异不大
plot(-log10(DEG_t.test[,4]),-log10(DEG[,4]))
logFC.png
p-value.png
-lg p.png
