三阴性乳腺癌表达矩阵探索笔记之差异性分析

生信技能树

发布于 2020-10-26 10:51:36

7560

发布于 2020-10-26 10:51:36

文章被收录于专栏：生信技能树生信技能树

学徒和学员已经陆续出师，是时候把生信技能树的舞台交给后辈了！

下面是学徒写的《GEO数据挖掘课程》的配套笔记（第3篇）

以第一个基因为例，根据group_list来检验在分组之间是否存在差异

load(file='step1-output.RData')
dat[1:4,1:4]
table(group_list) #分组信息
boxplot(dat[1,]~group_list)
t.test(dat[1,]~group_list) #p=0.00998有显著性差异

library(ggpubr)
df <- data.frame(gene=dat[1,], stage=group_list) #比较下一个基因可以改为dat[2,]
p <- ggboxplot(df, x = "stage", y = "gene",
               color = "stage", palette = "jco",
               add = "jitter")

p+stat_compare_means() #p值和之前不一样，因为换了一种统计学检验方法

以第一个基因为例进行表达差异性分析.Rplot

==Note== : 第一个基因是随机挑选的，虽然在两个类群中有差异性，但是从图上可以看出，noTNBC 有一部分是被包含在TNBC中的，并不是完全独立分离的关系，统计学显著性也不好说。

使用limma来进行批量的全部的基因的差异分析

#将绘制箱图的函数包装成函数便于使用
pb <- function(g){
  library(ggpubr)
  df <- data.frame(gene=g, stage=group_list) #比较下一个基因可以改为dat[2,]
  p <- ggboxplot(df, x = "stage", y = "gene",
                 color = "stage", palette = "jco",
                 add = "jitter")
  
  p+stat_compare_means() #p值和之前不一样，因为换了一种统计学检验方法
  
}
pb(dat[1,])

library(limma)
design <- model.matrix(~factor(group_list)) #design就是limma做差异分析需要的一个输入值，了解它可以有group_list构建
fit <- lmFit(dat, design)
fit <- eBayes(fit)
options(digits = 4)
topTable(fit,coef = 2,adjust="BH")
pb(dat["241662_x_at",]) #
deg <- topTable(fit, coef=2, adjust="BH", number = Inf) #Inf就是无穷大,把所有的数值都拿出来

limma筛选出来的在两个分组中差异表达比较显著的基因, logFC为负表示下调，以noTNBC作为对照

limma识别到的差异表达基因

以上面的第一个基因241662_x_at为例绘制箱线图，这个基因在两个分组之间的表达差异非常显著，而且没有重叠部分，TNBC和noTNBC完全分开了。

241662_x_at.Rplot

下面我们需要将这些差异基因的基因名找出来，也就是将探针ID转换到基因symbol

##进行注释
rm(list=ls())
load(file = "deg.Rdata")
deg$probe_id = rownames(deg)

library(hgu133plus2.db)
?hgu133plus2.db #学习包的内容
ids <- toTable(hgu133plus2ENSEMBL) #取出探针和基因名对应的数据集，并将其转化为表格

#merge函数，根据两个数据框的相同的列名进行合并
deg_1 <- merge(deg,ids,by="probe_id")

绘制火山图

将上调和下调的基因清楚地显示出来

nrDEG = deg_1
head(nrDEG)
attach(nrDEG)
plot(logFC,-log10(P.Value))
library(ggpubr)
df = nrDEG
df$v = -log10(P.Value)
ggscatter(df,x="logFC", y ='v', size = 0.5)
df$g = ifelse(df$P.Value > 0.01, "stable",
                ifelse(df$logFC>1.5,"up",
                       ifelse(df$logFC < -1.5,"down","stable")))
table(df$g)
df$name = rownames(df)
ggscatter(df,x = "logFC", y = "v", size=0.5, color = "g")
ggscatter(df,x ="logFC",y="v",color = "g",size=0.5,
            label = "symbol", repel = T,
            label.select = c("PROM1","GABRP","AGR3","SCGB2A2")) #标记比较显著的基因

火山图1

绘制热图

火山图不需要表达矩阵，只要差异分析结果的表格就可以

##绘制热图
load(file='step1-output.RData')
dat[1:4,1:4]
x = deg$logFC
names(x) = deg$probe_id
cg= c(names(head(sort(x),100)), names(tail(sort(x),100))) #分别挑出上调和下调表达最多的100个基因
library(pheatmap)
n <- t(scale(t(dat[cg,])))
n[n>2]=2
n[n< -2,]=-2
pheatmap(n,show_rownames = F,show_colnames = F)
ac=data.frame(g=group_list) #添加分组信息
rownames(ac)=colnames(n)
pheatmap(n, show_colnames = F,
         show_rownames = F,
         cluster_cols = F,
         cluster_rows = F,
         annotation_col = ac)

使用前表达差异性特别显著的基因绘制热图就可以很明显地看出基因表达的模式差异