R语言学习笔记-Day07

1 GEO数据挖掘

数据类型

(1) 基因表达芯片

(2) 转录组

(3) 单细胞

(4) 突变、甲基化、拷贝数变异……

怎样筛选基因--分析方法

数据下载（DEO、TCGA）-差异分析（芯片与转录组不相同）-WGCNA（加权共表达网络）-富集分析（ORA、GSEA）-PPI网络-预后分析（影响生存的疾病）

1.1

1.1.1 热图

输入数值为数值型矩阵/数据框

以颜色变化代表数值大小

#聚类树：根据基因相似程度进行排序分类，与原表达矩阵基因顺序不同

1.1.2 散点图和箱线图

可以用箱线图代替散点图，显示整体差异

箱线图：

以连续型向量为纵坐标；有重复值的离散型向量为横坐标

箱线图的五条线

max - 75% - median#中位数 - 25% - min

最大值和最小值以外可能存在离群值#离群点

#用于单个基因在几组之间的表达差异

###多基因 --> 差异分析

1.1.3 火山图

两个数值：logFC、P.Value

logFC（横坐标）

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)

P.Value

会进行调整将其增大

-->

-log10(P.Value)

P.Value越小，-log10(P.Value)越大，差异越大的置信度越高

1.1.4 主成分分析

PCA样本聚类图降维

点与点之间的相对距离表示相似程度

横、纵坐标：Dimension（Dim1、2）——主成分（综合指标）

几个基因组合到一起成为一个主成分

例如：BMI

#括号内的数字越大越好，没有具体要求

#图中最大的点为聚类的中心点，不是样本点

#至少四个样本点才能在图中形成一簇

#将权重最高的两个主成分作为横、纵坐标，而非全部主成分

#用于简单查看组间是否存在差异

2 GEO背景知识及芯片表达分析思路

2.1 表达数据实验设计

实验目的：通过基因表达量数据的差异分析和富集分析来解释生物学现象

2.2 数据库介绍

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* 并不是所有基因都存在差异

2.3 表达矩阵

行名：探针id #需要转换为gene symbol

列名：GSM，样本编号 #需要分组信息

3 代码分析流程

芯片差异分析所需输入数据

表达矩阵

#数据分布范围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.