(1) 基因表达芯片
(2) 转录组
(3) 单细胞
(4) 突变、甲基化、拷贝数变异……
数据下载(DEO、TCGA)-差异分析(芯片与转录组不相同)-WGCNA(加权共表达网络)-富集分析(ORA、GSEA)-PPI网络-预后分析(影响生存的疾病)
输入数值为数值型矩阵/数据框
以颜色变化代表数值大小
#聚类树:根据基因相似程度进行排序分类,与原表达矩阵基因顺序不同
可以用箱线图代替散点图,显示整体差异
箱线图:
以连续型向量为纵坐标;有重复值的离散型向量为横坐标
箱线图的五条线
max - 75% - median#中位数 - 25% - min
最大值和最小值以外可能存在离群值#离群点
#用于单个基因在几组之间的表达差异
###多基因 --> 差异分析
两个数值:logFC、P.Value
Foldchange(FC):处理组均值/对照组均值
log2Foldchange(logFC):Foldchange取log2
#实际运算中先取log再相减
#logFC表示处理组和对照组相比的基因表达差异倍数
#存在负值,表示表达降低
#基因的上调/下调,指基因表达量显著上升/下降
--> P.Value
芯片差异分析的起点是一个取过log的表达矩阵(0-20);
若未进行该操作,数值将非常大,需要先取log
通常设置阈值,例如:
上调基因:logFC > 1, p < 0.01
下调基因:logFC < -1, p < 0.01
#logFC常见阈值:1, 2, 1.2, 1.5, 2.2, 0.585 = log2(1.5)
会进行调整将其增大
-->
-log10(P.Value)
P.Value越小,-log10(P.Value)越大,差异越大的置信度越高
点与点之间的相对距离表示相似程度
横、纵坐标:Dimension(Dim1、2)——主成分(综合指标)
几个基因组合到一起成为一个主成分
例如:BMI
#括号内的数字越大越好,没有具体要求
#图中最大的点为聚类的中心点,不是样本点
#至少四个样本点才能在图中形成一簇
#将权重最高的两个主成分作为横、纵坐标,而非全部主成分
#用于简单查看组间是否存在差异
实验目的:通过基因表达量数据的差异分析和富集分析来解释生物学现象
Gene Expression Omibus
#Tools-Analyze a Study with GEO2R-代码,可以辅助完成一些操作
Sample:用户提交给GEO的样本数据(GSM)
Series:一个完整的研究,提供了整个研究的描述,包括对数据的描述、总结、分析(GSE)
Platform:用户测定表达量使用的芯片/平台(GPL)
基因表达芯片的原理:探针的表达量代表基因的表达量
#分析思路
找数据,找到GSE编号 -->
下载并读取数据 -->
#表达矩阵 #临床(分组)信息 #GPL编号(探针注释)
数据探索 -->
#分组间是否存在差异,PCA、热图
差异分析并可视化 -->
#P.Value, logFC #火山图、热图
富集分析
#KEGG #GO
为什么不画全部基因的热图
1* 数据太大
2* 并不是所有基因都存在差异
行名:探针id #需要转换为gene symbol
列名:GSM,样本编号 #需要分组信息
表达矩阵
#数据分布范围0-20
#无异常值,如NA、INF、负值
#无异常样本
分组信息
#同一分组对应同一关键词
#顺序与表达矩阵的列一一对应
#因子,对照组的levels在前
探针注释
#根据GPL编号查找
#探针与基因之间的对应关系
#只能有两列,且均为字符型
#列名必须是probe_id和symbol
options("repos"="https://mirrors.ustc.edu.cn/CRAN/") if(!require("BiocManager")) install.packages("BiocManager",update = F,ask = F) options(BioC_mirror="https://mirrors.westlake.edu.cn/bioconductor") cran_packages <- c('tidyr', 'tibble', 'dplyr', 'stringr', 'ggplot2', 'ggpubr', 'factoextra', 'FactoMineR', 'devtools', 'cowplot', 'patchwork', 'basetheme', 'paletteer', 'AnnoProbe', 'ggthemes', 'VennDiagram', 'survminer', "tinyarray") Biocductor_packages <- c('GEOquery', 'GO.db', 'hgu133plus2.db', 'ggnewscale', "limma", "impute", "GSEABase", "GSVA", "clusterProfiler", "org.Hs.eg.db", "preprocessCore", "enrichplot")
for (pkg in cran_packages){ if (! require(pkg,character.only=T,quietly = T) ) { install.packages(pkg,ask = F,update = F) require(pkg,character.only=T) } }
for (pkg in Biocductor_packages){ if (! require(pkg,character.only=T,quietly = T) ) { BiocManager::install(pkg,ask = F,update = F) require(pkg,character.only=T) } }
#前面的所有提示和报错都先不要管。主要看这里
for (pkg in c(Biocductor_packages,cran_packages)){ require(pkg,character.only=T) } #没有任何提示就是成功了,如果有warning xx包不存在,用library检查一下。 #library报错,就单独安装。
rm(list = ls()) #打破下载时间的限制,改前60秒,改后10w秒 options(timeout = 100000) options(scipen = 20)#不要以科学计数法表示
#传统下载方式
library(GEOquery) eSet = getGEO("GSE7305", destdir = '.', getGPL = F) #网速太慢,下不下来怎么办 #1.从网页上下载/发链接让别人帮忙下,放在工作目录里 #2.试试geoChina(),只能下载2019年前的表达芯片数据 library(AnnoProbe) eSet = geoChina("GSE7305") #选择性代替getGEO() #研究一下这个eSet class(eSet) length(eSet)
eSet = eSet[1] class(eSet) 1 "ExpressionSet" attr(,"package") 1 "Biobase" #“ExpressionSet”为来自R包“Biobase”中的一个对象
#(1)提取表达矩阵exp
exp <- exprs(eSet)
#⭐第一个要检查的地方👇,表达矩阵行列数,正常是几万行,列数=样本数,
#如果0行说明不是表达芯片或者是遇到特殊情况,不能用此流程分析
dim(exp)
#⭐二个要检查的地方👇
range(exp)#看数据范围决定是否需要log,是否有负值,异常值,如有负值,结合箱线图进一步判断 #数据范围应为0-20之间 #0-4可能取了两次log2,其它情况也有可能取成log10
#⭐可能要修改的地方👇
exp = log2(exp+1) #需要log才log,不需要log要注释掉这一句
#⭐第三个要检查的地方👇
boxplot(exp,las = 2) #箱线图看是否有异常样本 #应当在大概相等的范围内 #处理异常样本 第一个办法:删掉异常样本 第二个办法:exp = limma::normalizeBetweenArrays(exp) #中位数在0附近,是不正常的标准化数据#做过不可逆操作,无法继续分析 #取过log,存在少量负值,4<中位数<15——正常 #没取log,有负值——错误数据
#(2)提取临床信息
pd <- pData(eSet) #临床信息表格中的行为表达矩阵的列 #⭐多分组中提取两分组的代码示例,二分组不需要 if(F){ #因为现在这个例子不是多分组,所以编造一列做示例。 pd$fake = paste0(rep(c("a","b","c","d"),each = 5),1:5)#用于编造一列,实际操作中无需此句 k1 = str_detect(pd$fake,"b");table(k1) k2 = str_detect(pd$fake,"c");table(k2) pd = pdk1|k2, }
#(3)让exp列名与pd的行名顺序完全一致
p = identical(rownames(pd),colnames(exp));p if(!p) { s = intersect(rownames(pd),colnames(exp)) exp = exp,s pd = pds, }
#(4)提取芯片平台编号,后面要根据它来找探针注释
gpl_number <- eSet@annotation;gpl_number save(pd,exp,gpl_number,file = "step1output.Rdata") #"@"为对象提取子集的编号
1* Bioconductor的注释包
2* GPL页面的表格文件解析
3* 官网下载对应产品的注释表格
4* 自主注释
#不是所有GPL都能找到注释
代码及部分相关注释源自课程脚本
引用自生信技能树
This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
原创声明:本文系作者授权腾讯云开发者社区发表,未经许可,不得转载。
如有侵权,请联系 cloudcommunity@tencent.com 删除。
原创声明:本文系作者授权腾讯云开发者社区发表,未经许可,不得转载。
如有侵权,请联系 cloudcommunity@tencent.com 删除。