GEO数据挖掘-DAY1

常用数据库：

1、不分方向——GEO、NHANES（临床资料）
2、肿瘤专属——TCGA、ICGC、CCLE、SEER

一般可挖掘的数据：基因表达芯片、转录组学、单细胞、甲基化、拷贝数变异...

（基因表达芯片和转录组学虽然都是测定表达水平，但是分析方法和数据格式并不相同）

如何去筛选基因：

1. 数据下载
2. 差异分析（表达芯片array和转录组分析方法不一样）
3. WGCNA（加权共表达网络）
4. 富集分析（ORA、GSEA）
5. PPI网络（非R语言工作）
6. 预后分析（疾病存在一个影响生存/疾病生化复发的结局）

图表的简介：（具体代码和图表解读待补充）

1. 热图（heatmap/pheatmap）- 必须输入数值型矩阵/数据框的格式，如果热图中基因较多，需要使用complicated heatmap去标注后文研究相关的特征基因。
2. 散点图
3. 箱线图- 输入连续型向量（一般是纵轴）和一个有重复值的离散型向量（一般是横轴），箱线图的五条线从上至下分别为max，75%，median，25%，min。但是max和min并非真正意义上的最大值和最小值，因为其考虑到了离群值，可能会有有些点落在max和min之外，可以表示单个基因在两组之间的表达量差异（散点图看不清整体趋势）。在表达矩阵之外（行为gene，列为sample），还需要额外添加分组信息，因此需要注意对应关系，哪些是实验组哪些是对照组。
4. 火山图（多基因差异分析，logFC-P.Value/adj.P.Val）

 

• Foldchange（FC）=实验组平均值/对照组平均值，一般需要取log2FC来展示，一般log2FC范围在0-25之间，也就是log2(实验组) - log2(对照组)，顺序错了上调下调就全反了！！！
• 横坐标是logFC，纵坐标是-log10（P.Value / adj.P.Val）
• 表达芯片进行差异分析的出发点就是一个取过log的表达矩阵（有些已经取过了，转录组的数据分析没有这个要求）
• log2FC＞0，即 log2实验组 - log2对照组＞0，即基因表达上升，反之则为下降。但是一般我们在将表达量上升or下降的基因定义为上调or下调基因时，会自定义一个logFC和P.Value，比如常将上调基因和下调基因阈值分别设置为log2FC＞1，P＜0.01；log2FC＜-1，P＜0.01。如果数据不太争气，差异不太明显，整成log2FC＞0.585（1.5倍差异以上）。
• 火山图的解读：横坐标越偏向两侧，上调or下调越显著；纵坐标越往上，该数据越可信。logFC=4，说明该基因实验组表达量是对照组的2^4倍。
• Add：图做出来之后看一下横坐标的logFC，看看有没有取log，还是重复取log了（如果在阈值调到很低的情况下还是找不出什么差异，很有可能是在别人取了log的情况下又取了一遍）。

主成分分析PCA样本聚类图（具体可见另一篇独立笔记）

        总而言之，使用降维的算法，自定义坐标轴——把多个指标转化为少数几个综合指标（即主成分），根据这些主成分进行样本聚类，进而实现样本组内/组间的差异可视化。同组内点聚集→组内重复好，中心点距离远→组间差异明显。

GEO数据库使用之基因表达芯片数据

1. browse content → series → Expression profiling by array（基因表达芯片）
2. 设定filter条件→series type/organism（homo or mouse）/samples/...
3. 下载数据集，点击matrix下载 → GSExxxx series matrix.txt.gz文件格式，转录组和单细胞是在supplementary file下面

     4.点进样本名，然后去看一下表达矩阵中的探针注释/value/行数

 

1. 找到合适数据，下载数据（表达矩阵、临床信息/分组信息、GPL编号、探针注释）
2. 初步数据探索，查看分组之间是否有差异（PCA、热图）
3. 差异分析及可视化 → P值、LogFC → 火山图、热图

4. 富集分析KEGG、GO

 

 

#传统下载方式
library(GEOquery)
eSet = getGEO("GSE7305", destdir = '.', getGPL = F)
#网速太慢，下不下来怎么办
#1.从网页上下载后放在工作目录里
#2.试试geoChina（下面这个包）,只能下载2019年前的表达芯片数据
#library(AnnoProbe)  #几百个样本的时候这个包下的很快
#eSet = geoChina("GSE7305")  #选择性代替eSet = getGEO()

#研究一下这个eSet
class(eSet)
length(eSet)

eSet = eSet[[1]] #eset只有一列，直接提取出去
class(eSet)

#(1)提取表达矩阵exp
exp <- exprs(eSet)
dim(exp)      #看行和列数，如果是0行说明没导入表达矩阵，也可能下的这个文件里没有。
range(exp)   #看数据范围决定是否需要自己去手动log，是否有负值，异常值
exp = log2(exp+1)     #需要log才log
boxplot(exp,las = 2)         #看是否有异常样本（跟其他样本截然不同），还要看纵坐标有无log       （纵坐标范围在2-4之间大概率是取重复了）
exp=limma::normalizeBetweenArrays(exp)          #处理异常样本的函数，还有一种办法就是直接删除异常样本

#关于表达矩阵里的负值
##1、取过log，有少量的负值——正常
##2、没取log，存在负值——错误数据，弃用or手动处理原始数据
##3、一半负值，中位数为0——原作者使用zscore进行了标准化，只能拿来做热图，一般弃用
# 原始数据处理的代码，按需学习
# https://mp.weixin.qq.com/s/0g8XkhXM3PndtPd-BUiVgw

#(2)提取临床信息
pd <- pData(eSet)
view(pd)      #一般看title这一列（也可能是在别的列），哪些是实验组，哪些是对照组，并找到提取的关键词
#(3)让exp列名与pd的行名顺序完全一致，
p = identical(rownames(pd),colnames(exp));p #先判断是否一致，再使用 if 来决定后续一致化代码是否运行
if(!p) {
  s = intersect(rownames(pd),colnames(exp))
  exp = exp[,s]
  pd = pd[s,]
}

#(4)提取芯片平台编号，后面要根据它来找探针注释
gpl_number <- eSet@annotation;gpl_number
save(pd,exp,gpl_number,file = "step1output.Rdata")