常用数据库:
1、不分方向——GEO、NHANES(临床资料)
2、肿瘤专属——TCGA、ICGC、CCLE、SEER
一般可挖掘的数据:基因表达芯片、转录组学、单细胞、甲基化、拷贝数变异...
(基因表达芯片和转录组学虽然都是测定表达水平,但是分析方法和数据格式并不相同)
如何去筛选基因:
- 数据下载
- 差异分析(表达芯片array和转录组分析方法不一样)
- WGCNA(加权共表达网络)
- 富集分析(ORA、GSEA)
- PPI网络(非R语言工作)
- 预后分析(疾病存在一个影响生存/疾病生化复发的结局)
图表的简介:(具体代码和图表解读待补充)
- 热图(heatmap/pheatmap)- 必须输入数值型矩阵/数据框的格式,如果热图中基因较多,需要使用complicated heatmap去标注后文研究相关的特征基因。
- 散点图
- 箱线图- 输入连续型向量(一般是纵轴)和一个有重复值的离散型向量(一般是横轴),箱线图的五条线从上至下分别为max,75%,median,25%,min。但是max和min并非真正意义上的最大值和最小值,因为其考虑到了离群值,可能会有有些点落在max和min之外,可以表示单个基因在两组之间的表达量差异(散点图看不清整体趋势)。在表达矩阵之外(行为gene,列为sample),还需要额外添加分组信息,因此需要注意对应关系,哪些是实验组哪些是对照组。
- 火山图(多基因差异分析,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数据库使用之基因表达芯片数据
- browse content → series → Expression profiling by array(基因表达芯片)
- 设定filter条件→series type/organism(homo or mouse)/samples/...
- 下载数据集,点击matrix下载 → GSExxxx series matrix.txt.gz文件格式,转录组和单细胞是在supplementary file下面
4.点进样本名,然后去看一下表达矩阵中的探针注释/value/行数
- 找到合适数据,下载数据(表达矩阵、临床信息/分组信息、GPL编号、探针注释)
- 初步数据探索,查看分组之间是否有差异(PCA、热图)
- 差异分析及可视化 → P值、LogFC → 火山图、热图
-
富集分析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")