学习笔记Day8：GEO数据挖掘-基因表达芯片

GEO数据挖掘

• 数据库：GEO、NHANCE、TCGA、ICGC、CCLE、SEER等

• 数据类型：基因表达芯片、转录组、单细胞、突变、甲基化、拷贝数变异等等

• 

常见图表

表达矩阵

• 一行为一个基因，一列为一个样本，内容是基因表达量。

热图

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

颜色变化表示数值大小

• Complexheatmap：自行探索~可以定义某些需要标注的基因。

散点图和箱线图

1. 数据要求：xy轴分别为一个连续型向量和一个有重复值的离散型向量（分类向量）

2. 数据意义：箱线图会删除离群点

   

3. 意义：单个基因在两组之间的表达量差异

   矩阵/数据框中无法直接添加分组信息，因此需要一个单独向量来说明样本的分组信息。

   

火山图

1. 意义：展示多个基因在两组间的表达量差异

2. 横坐标：FC(Foldchange) 处理组平均值/对照组平均值；logFC（log2Foldchange）

   • 芯片差异分析的起点是一个取过log的表达矩阵，如果未取log需要自行log

     

3. 理解logFC：log2(X/Y)=log2(x)-log2(y)

   • log后的表达矩阵：表达量在0-24之间。
   • 未log的表达矩阵，表达量在0，10，100，1000…
   • logFC的正常范围：个位数居多。

4. logFC>0，treat>control，基因表达量上升；logFC<0，treat<control，基因表达量下降。通常说的上调和下调基因是指表达量显著上升下降的基因。

5. logFC常见的阈值：1/2/1.2/1.5/2.2/0.585[log2(1.5)]；需要根据情况取值。

   

6. 纵坐标：P.Value→log10(P.Value)→-log10(P.Value)

   

主成分分析

1. PCA：主成分，把多个指标转换成少数综合指标（主成分）来代表样本。代表样本的点在坐标轴上距离越远，说明样本差异越大。

2. PCA样本聚类图：dim1/dim2中数字不重要（尽量大）。

• 图上点代表样本（中心点除外），点与点之间距离代表样本差异。

3. 可以用于“预实验”，简单看组间是否有差别。

表达芯片分析思路

表达数据实验设计

1. 实验目的：通过基因表达量数据的差异分析和富集分析来解释生物学现象。
   • 设计实验组和对照组
2. 有差异的材料→差异基因→找功能→解释差异，缩小差异范围。

基因表达芯片

探针的表达量代表基因的表达量

探针根据要测量的基因设计，是一段与基因互补杂交的短核苷酸序列。探针和序列绑定，不和基因绑定。

数据库介绍

1. GEO数据库

   • GEO工具：GEO2R，可以导入代码进入R修改。

2. Series：用户提交给数据库一个完整的研究，包括其样本数据(GSM)，包含提供研究描述，包括对数据描述并总结分析(GSE)。

3. GEO数据集筛选

   表达芯片数据：Expression profiling by array

   单细胞/普通转录组（高通量测序）：Expression profiling by high throughput sequencing

4. GSE界面：GPL（平台）中看ID和Symbol Gene；GSM（样本）中看表达量是否正常、是否需要取log等。

分析思路

1. 找数据，找到GSE编号
   • GEO数据库中检索
   • 文献中查找GSE编号
2. 下载数据：表达矩阵、临床信息（分组信息）、GPL编号（探针注释）
   • 网页中点选下载
   • 代码下载（推荐）
3. 数据探索：分组之间是否有差异、PCA、热图（方差排名靠前的1k个基因）
4. 差异分析和可视化：P值、logFC；火山图、热图
5. 富集分析：KEGG、GO

表达矩阵

以分组为单位说问题，而不是以样本为单位

代码分析流程

安装R包

下载数据

library(GEOquery)
eSet = getGEO("GSE7305", destdir = '.', getGPL = F)##探索eSet
class(eSet)
##[1] "list"
length(eSet)
##[1] 1
eSet = eSet[[1]]     ##将list解开
class(eSet)          ##是一种特殊的数据类型，可以从帮助文档中找到详细说明。出自Biobase包
##[1] "ExpressionSet"
##attr(,"package")
##[1] "Biobase"

提取表达矩阵

exp <- exprs(eSet)    ##提取表达矩阵
dim(exp)              ##查看数据属性
range(exp)            ##查看数据范围，决定是否需要取log，是否有负值
##[1]     5.020951 22011.934000       ##这样的数据就需要取log
exp = log2(exp+1)     ##取log
boxplot(exp,las = 2)  ##检查数据情况

• 有异常样本：1. 删掉异常样本；2. limma包的标准化函数拉齐。
• 负值：log后少量负值可以接受，log前有负值或一半都是负值（标准化）弃用。
• 如果logFC在2-4之间，有可能取了2次log

提取临床信息

pd <- pData(eSet)       ##包含分组信息

让表达矩阵和临床信息顺序匹配

p = identical(rownames(pd),colnames(exp));p    ##判断是否对应if(!p) {                                       ##如果不对应，运行s = intersect(rownames(pd),colnames(exp))   exp = exp[,s]pd = pd[s,]
}

提取芯片平台编号

gpl_number <- eSet@annotation;gpl_number        ##等同于网页查找
save(pd,exp,gpl_number,file = "step1output.Rdata")  ##保存数据

• 原始数据处理方法（有时间再整理吧~）：https://mp.weixin.qq.com/s/0g8XkhXM3PndtPd-BUiVgw

引用自生信技能树课程~ 给小洁老师比心~