http://www.bio-info-trainee.com/3409.html

1、安装一些R包

数据包： ALL, CLL, pasilla, airway
软件包：limma，DESeq2，clusterProfiler
工具包：reshape2
绘图包：ggplot2

options("repos"="https://mirrors.ustc.edu.cn/CRAN/")
if(!require("BiocManager")) install.packages("BiocManager",update = F,ask = F)
options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")
#包的名字
cran_packages <- c("reshape2",'ggplot2')
Biocductor_packages <- c("limma","DESeq2","clusterProfiler","ALL","CLL","pasilla","airway")

for (pkg in Biocductor_packages){
  if (! require(pkg,character.only=T) ) {
    BiocManager::install(pkg,ask = F,update = F)
    require(pkg,character.only=T) 
  }
}

for (pkg in cran_packages) {
  if(!require(pkg,character.only = T)){
  install.packages(pkg,ask = F,update = F)
  require(pkg,character.only=T)   
    }
  }
 for (pkg in c(cran_packages,Biocductor_packages)) {require(pkg,character.only=T) }

只要最后一条代码运行没错就完全🆗，不用管warning message.

2、了解ExpressionSet对象，比如CLL包里面就有data(sCLLex) ，找到它包含的元素，提取其表达矩阵(使用exprs函数)，查看其大小

http://www.bio-info-trainee.com/1510.html
这个很好的解释了什么是ExpressionSet对象

简单来说就是一个包含了表达矩阵及样本分组的一个封装，就是一大堆信息，像一个大箱子，里面装了很多你分析数据需要用的东西。

library(CLL)
print(data(package = "CLL"))
?sCLLex
data <- data.frame(exprs(sCLLex))
dim(data)

3.了解 str,head,help函数，作用于 第二步提取到的表达矩阵

head(data)#查看前6行
help(sCLLex)#帮助文档
library(stringr)
str_length(rownames(data)[3])#str函数有很多

4.安装并了解 hgu95av2.db 包,看看 ls("package:hgu95av2.db") 后 显示的那些变量

#BiocManager::install("hgu95av2.db")
library(hgu95av2.db)
ls("package:hgu95av2.db")#里面是注释包的一些内容，这个包就是一个盒子，里面很多包裹

5.理解 head(toTable(hgu95av2SYMBOL)) 的用法，找到 TP53 基因对应的探针ID

SYMBOL <- toTable(hgu95av2SYMBOL)
?toTable#toTable() is an S4 generic function provided as an alternative to as.data.frame().
SYMBOL2 <- as.data.frame(hgu95av2SYMBOL)
head(SYMBOL)
SYMBOL[SYMBOL$symbol == "TP53",]

6.理解探针与基因的对应关系，总共多少个基因，基因最多对应多少个探针，是哪些基因，是不是因为这些基因很长，所以在其上面设计多个探针呢？

length(unique(SYMBOL$symbol))#8582个基因
length(unique(SYMBOL$probe_id))#11457个探针
#为什么多个探针对应一个基因，我得理解是因为一个基因有很多编码区片段，自然就会有很多探针，所以取舍可以平均值、最大值、中位值之类的；
#还有多个基因对应一个探针，那就是序列重叠，这几个基因都可以用，或者直接取一个。

7.第二步提取到的表达矩阵是12625个探针在22个样本的表达量矩阵，找到那些不在 hgu95av2.db 包收录的对应着SYMBOL的探针。

#提示：有1165个探针是没有对应基因名字的
data2 <- data[!(rownames(data) %in% SYMBOL$probe_id),]#记住！的用法就是取相反的，然后%in%巨无敌好用

8.过滤表达矩阵，删除那1165个没有对应基因名字的探针。

dat <- data[(rownames(data) %in% SYMBOL$probe_id),]

9. 整合表达矩阵，多个探针对应一个基因的情况下，只保留在所有样本里面平均表达量最大的那个探针。

# 提示，理解 tapply,by,aggregate,split 函数 , 首先对每个基因找到最大表达量的探针。
# 然后根据得到探针去过滤原始表达矩阵
dat <- data[(rownames(data) %in% SYMBOL$probe_id),]
SYMBOL <- SYMBOL[match(rownames(dat),SYMBOL$probe_id),]#调整顺序，让SYMBOL的和dat的一样
SYMBOL$median <- apply(dat, 1, median)

SYMBOL2 <- SYMBOL[order(SYMBOL$probe_id,SYMBOL$median,decreasing = T),]#先按照基因字母顺序排序，同样的基因就在一起，同样的基因一起之后按照中位数排序最大的在最上面
SYMBOL2 <- SYMBOL2[!duplicated(SYMBOL2$symbol),]
dat <- dat[(rownames(dat) %in% SYMBOL2$probe_id),]

10.把过滤后的表达矩阵更改行名为基因的symbol，因为这个时候探针和基因是一对一关系了。

rownames(dat) <- SYMBOL2$symbol

11.对第10步得到的表达矩阵进行探索，先画第一个样本的所有基因的表达量的boxplot,hist,density ， 然后画所有样本的这些图

col = rainbow(15)
boxplot(dat$CLL11.CEL,col = "yellow")#箱线图
hist(dat$CLL11.CEL,col = col,main = paste("Histogram of" , "CLL11.CEL"),xlab = "CLL11.CEL")#彩虹的柱状图嘻嘻
#密度图、柱状图
hist(dat$CLL11.CEL,col = col,main = paste("Histogram of" , "CLL11.CEL"),xlab = "CLL11.CEL",freq = F)
lines(density(dat$CLL11.CEL),col = "red")
#ggplot2
library(ggplot2)
library(reshape2)
##整理数据
pd <- pData(sCLLex)
group_list <- factor(pd$Disease,levels = c("stable","progres."))
dat$prob_id <- rownames(dat) 
dat_l =melt(dat)
colnames(dat_l)=c('probe','sample','value')
dat_l$group=rep(group_list,each=nrow(dat))
head(dat_l)
#画图
p <- ggplot(dat_l,aes(x = sample,y = value,fill = group))
p + geom_boxplot()#箱线图
p <- ggplot(dat_l,aes(value,fill = group))
p + geom_histogram(bins = 200)+facet_wrap(~sample, nrow = 4)
p <- ggplot(dat_l,aes(value,col = group))+ geom_density()+facet_wrap(~sample, nrow = 4)

image.png

image.png

image.png

12.理解统计学指标mean,median,max,min,sd,var,mad并计算出每个基因在所有样本的这些统计学指标,最后按照mad值排序，取top 50 mad值的基因，得到列表。

dat <- dat[,-ncol(dat)]#刚刚加了一列，现在去掉，恢复原文件
me <- apply(dat, 1, mean)
med <- apply(dat, 1, median)
max <- apply(dat, 1, max)
min <- apply(dat,1,min)
sd <- apply(dat, 1, sd)
var <- apply(dat, 1, var)
mad <- apply(dat, 1, mad)#就是先求出给定数据的中位数(注意并非均值)然后原数列的每个值与这个中位数求出绝对差,然后新数列的中位值就是MAD
all <- data.frame(me,med,max,min,sd,var,mad,row.names = rownames(dat))
all_ord <- all[order(all$mad,decreasing = T),]
all_ord <- all_ord[1:50,]

13.根据第12步骤得到top 50 mad值的基因列表来取表达矩阵的子集，并且热图可视化子表达矩阵。试试看其它5种热图包的不同效果。

library(pheatmap)
library(ggplot2)
heatmap50 <-all_ord 
class(heatmap50)
heatmap50 <- as.matrix(heatmap50)#数据框直接做提示Error in heatmap(all_ord) : 'x'必需为数值矩阵
heatmap50 <- scale(heatmap50)#R语言中scale函数，可以对数据进行处理，标准化(归一化)在一定的范围，比较适合大范围变化数据归一化处理从而观察数据变化趋势
pheatmap(all_ord)#1
heatmap(all_ord)#2
library(gplots)
heatmap.2(heatmap50,scale = "none", col=bluered(100), 
          trace = "none", density.info = "none")#3
#整理数据
all_gg <- all_ord
all_gg$ID <- rownames(all_gg)
all_gg <- melt(all_gg)
all_gg$value <- scale(all_gg$value)
p <- ggplot(all_gg,aes(x=variable, y = ID,fill = value)) + geom_tile() + scale_fill_gradient(low = "white", high = "red")
p#4

image.png

image.png

image.png

image.png

14.取不同统计学指标mean,median,max,mean,sd,var,mad的各top50基因列表，使用UpSetR包来看他们之间的overlap情况。

if (!require("UpSetR")) {install.packages(UpSetR)
  
}
library(UpSetR)#集合可视化的神包，像韦恩图，但是可以提供更多的集合信息
mean <- rownames(data.frame(sort(me,decreasing = T)[1:50]))
median <- rownames(data.frame(sort(med,decreasing = T)[1:50]))
max<- rownames(data.frame(sort(max,decreasing = T)[1:50]))
sd <- rownames(data.frame(sort(sd,decreasing = T)[1:50]))
var <-rownames(data.frame(sort(var,decreasing = T)[1:50]))
mad <- rownames(data.frame(sort(mad,decreasing = T)[1:50]))
n_all <- unique(c(mean,median,max,sd,var,mad))
data_upsetR <- data.frame(mean = ifelse(n_all %in% mean ,1,0),
                          median = ifelse(n_all %in% median ,1,0),
                          max = ifelse(n_all %in% max ,1,0),
                          sd = ifelse(n_all %in% sd ,1,0),
                          var = ifelse(n_all %in% var ,1,0),
                          mad = ifelse(n_all %in% mad ,1,0))
rownames(data_upsetR) <- n_all
upset(data_upsetR,nsets = 6)

image.png

15.在第二步的基础上面提取CLL包里面的data(sCLLex) 数据对象的样本的表型数据。

pd <- pData()

16.对所有样本的表达矩阵进行聚类并且绘图，然后添加样本的临床表型数据信息(更改样本名)

#准备数据
head(data)
pd <- pData(sCLLex)
colnames(data) <- paste0(pd$Disease,str_sub(pd$SampleID,4))
mydata <- na.omit(data) 
mydata <- t(scale(mydata))#转换数据变成行是样本，列是基因
d <- dist(mydata)#计算距离
fit2 <- hclust(d,method = "ward.D")#聚类的主要函数
plot(as.dendrogram(fit2))#画图

image.png

17.对所有样本的表达矩阵进行PCA分析并且绘图，同样要添加表型信息。

https://www.sohu.com/a/308627551_120055884

（推荐一个写的可以的）

#数据整理
data[1:3,1:3]
pd <- pData(sCLLex)
colnames(data) <- paste0(pd$Disease,str_sub(pd$SampleID,4))
mydata <- na.omit(data) 
mydata <- data.frame(t(mydata))#转换数据变成行是样本，列是基因
mydata[1:3,1:3]
mydata$group <- group_list
data.pca <- prcomp(mydata[,1:(ncol(mydata)-1)])
summary(data.pca)
head(data.pca$x)
#可视化
##
library(devtools)
library(ggplot2)
install_github('sinhrks/ggfortify')
library(ggfortify)
autoplot(data.pca,data = mydata,colour = 'group',label = TRUE)

image.png

18.根据表达矩阵及样本分组信息进行批量T检验，得到检验结果表格

#数据整理
pd <- pData(sCLLex)
group_list <- factor(pd$Disease,levels = c("stable","progres."))
data <- data.frame(exprs(sCLLex))
data_stable <- data[,which(group_list == "stable")]
data_progres. <- data[,which(group_list == "progres.")]
pvals <- apply(data, 1, function(x){t.test(as.numeric(x)~group_list)$p.value
})#这里注意要用原始数据，这样才是group_list对应的分组
p.adj = p.adjust(pvals, method = "BH")
avg_1 = rowMeans(data_stable)
avg_2 = rowMeans(data_progres.)
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)
# avg_1    avg_2     log2FC        pvals     p.adj
# 36129_at 8.791753 7.875615 -0.9161377 1.629755e-05 0.2057566
# 37676_at 7.965007 6.622749 -1.3422581 4.058944e-05 0.2436177
# 33791_at 5.786041 7.616197  1.8301554 6.965416e-05 0.2436177
# 39967_at 2.152471 4.456446  2.3039752 8.993339e-05 0.2436177
# 34594_at 7.058738 5.988866 -1.0698718 9.648226e-05 0.2436177
# 32198_at 3.407405 4.157971  0.7505660 2.454557e-04 0.3516678

19、使用limma包对表达矩阵及样本分组信息进行差异分析，得到差异分析表格

重点看logFC和P值，画个火山图(就是logFC和-log10(P值)的散点图。)

library(limma)
# 最重要是构建分组信息，构建好比较矩阵是关键
# 注意这里的表达矩阵信息行为gene，列为sample
#整理数据
pd <- pData(sCLLex)
group_list <- factor(pd$Disease,levels = c("stable","progres."))
data <- data.frame(exprs(sCLLex))
#构建分组信息
design <- model.matrix(~0+factor(group_list))
rownames(design) <- colnames(data)
colnames(design) <- c("stable","progres.")
design
contrast.matrix <- makeContrasts("progres.-stable",levels =design )
contrast.matrix
#这个矩阵我们要把progres.组跟stable进行差异分析比较
#拟合模型
fit <- lmFit(data,design)
fit2 <- contrasts.fit(fit,contrast.matrix)
fit2 <- eBayes(fit2)
tempOutput = topTable(fit2, coef=1, n=Inf)
DEG = na.omit(tempOutput)
head(DEG)
plot(DEG$logFC,-log10(DEG$P.Value))

image.png

20.对T检验结果的P值和limma包差异分析的P值画散点图，看看哪些基因相差很大。

DEG_t.test <- DEG_t.test[match(rownames(DEG),rownames(DEG_t.test)),]
plot(-log10(DEG_t.test$pvals),-log10(DEG$P.Value))#可以看到结果是差不多的
#整理行名为基因名
SYMBOL <- toTable(hgu95av2SYMBOL)
DEG_t.test2 <- DEG_t.test[!(rownames(DEG_t.test) %in% SYMBOL$probe_id),]
DEG_t.test <- DEG_t.test[(rownames(DEG_t.test) %in% SYMBOL$probe_id),]
DEG_t.test <- DEG_t.test[(rownames(DEG_t.test) %in% SYMBOL$probe_id),]
SYMBOL <- SYMBOL[match(rownames(DEG_t.test),SYMBOL$probe_id),]
SYMBOL$median <- apply(DEG_t.test, 1, median)
SYMBOL2 <- SYMBOL[order(SYMBOL$probe_id,SYMBOL$median,decreasing = T),]
SYMBOL2 <- SYMBOL2[!duplicated(SYMBOL2$symbol),]
DEG_t.test <- DEG_t.test[(rownames(DEG_t.test) %in% SYMBOL2$probe_id),]
rownames(DEG_t.test) <- SYMBOL2$symbol
##
SYMBOL <- toTable(hgu95av2SYMBOL)
DEG2 <- DEG[!(rownames(DEG) %in% SYMBOL$probe_id),]
DEG <- DEG[(rownames(DEG) %in% SYMBOL$probe_id),]
DEG <- DEG[(rownames(DEG) %in% SYMBOL$probe_id),]
SYMBOL <- SYMBOL[match(rownames(DEG),SYMBOL$probe_id),]
SYMBOL$median <- apply(DEG, 1, median)
SYMBOL2 <- SYMBOL[order(SYMBOL$probe_id,SYMBOL$median,decreasing = T),]
SYMBOL2 <- SYMBOL2[!duplicated(SYMBOL2$symbol),]
DEG <- DEG[(rownames(DEG) %in% SYMBOL2$probe_id),]
rownames(DEG) <- SYMBOL2$symbol
#
identical(rownames(DEG),rownames(DEG_t.test))
data_p <- cbind(DEG_t.test$pvals,DEG$P.Value)
rownames(data_p) <- rownames(DEG)
colnames(data_p) <- c("DEG_t.test$pvals","DEG$P.Value")
head(data_p)
plot(data_p)
#不明白题目的意思

image.png