专栏首页生信技能树三阴性乳腺癌表达矩阵探索笔记之差异性分析

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

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

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

  1. B站课程《三阴性乳腺癌表达矩阵探索》笔记之文献解读
  2. 三阴性乳腺癌表达矩阵探索之数据下载及理解

以第一个基因为例,根据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)

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

top100_rail100热图.Rplot01

视频观看方式

我把3年前的收费视频课程:3年前的GEO数据挖掘课程你可以听3小时或者3天甚至3个月,免费到B站:

  • 这个课程超级棒,B站免费学习咯:https://m.bilibili.com/video/BV1dy4y1C7jz
  • 配套代码在GitHub哈:https://github.com/jmzeng1314/GSE76275-TNBC
  • TCGA数据库挖掘,代码在:https://github.com/jmzeng1314/TCGA_BRCA
  • GTEx数据库挖掘,代码在:https://github.com/jmzeng1314/gtex_BRCA
  • METABRIC数据库挖掘,代码在:https://github.com/jmzeng1314/METABRIC

然后马上就有了3千多学习量,而且有学员给出来了图文并茂版本万字笔记,让我非常感动!

本文分享自微信公众号 - 生信技能树(biotrainee),作者:生信技能树

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2020-10-06

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 三阴性乳腺癌表达矩阵探索笔记之GSEA

    基于超几何分布检验的富集分析做KEGG数据库的时候,它总共只有七千多个基因,人类总的背景基因有两万多个,被KEGG记住的只有6500个(一直在增加),假设一条通...

    生信技能树
  • 这个WGCNA作业终于有学徒完成了!

    共画了3张热图,最后一张热图展示如下图,与原文对比'Ligamentocyte'和'Chondrocyte'相比较其他组是高表达的。

    生信技能树
  • 有趣的基因命名

    还有,如果你看到HS.开头的基因,它是unigene的ID了,已经不再是symbol啦。

    生信技能树
  • 将数据文件(csv,Tsv)导入Hbase的三种方法

    (1)使用HBase的API中的Put是最直接的方法,但是它并非都是最高效的方式(2)Bulk load是通过一个MapReduce Job来实现的,通过Job...

    黑白格
  • python爬取豆瓣电影Top250的信息

    拓荒者
  • JVM之SerialOld收集器

    WindWant
  • 【设计模式】几需体验三欢钟,里造会干我一样理解Facade和Mediator模式

    其实,现实生活中,我们可能远不止跑这三个部门,如果Client来回穿梭于N个部门间,办事效率是不是很垃圾了!

    行百里er
  • 使用Kafka订阅Binlog之字段值获取防坑指南(阿里云DTS)

    在《如果可以,我想并行消费Kafka拉取的数据库Binlog》这篇文章中,笔者介绍如何实现并行消费Binlog,具体实现就是让同一张表的Binlog放到同一个线...

    Java艺术
  • NodeController 源码分析

    在早期的版本中 NodeController 只有一种,v1.16 版本中 NodeController 已经分为了 NodeIpamController 与 ...

    田飞雨
  • Oracle参数解析(nls_sort)

    前面介绍了Oracle的基本参数,从这节开始讲其他的参数,参数从v$parameter中提取

    bsbforever

扫码关注云+社区

领取腾讯云代金券