R语言转录组分析-GEO数据库下载及Illumina数据处理

一、简单认识一下GEO数据库

入口：在NCBI主页 http://www.ncbi.nlm.nih.gov/ 的搜索入口输入关键字，比如某个基因或某个疾病，选择数据类型为GEO DataSets。如果想直接看某个基因在一些不同的GEO数据集中表达量，可以直接选择数据类型为GEO Profiles。

详细的页面介绍，搜索方法，参数设置，数据介绍可以看这里https://www.omicsclass.com/article/1101

GEO Platform (GPL) 芯片平台样品记录描述了处理单个样品的条件，它经历的操作以及从中得到的每个元素的丰度测量。每个样品记录都分配有唯一且稳定的GEO登录号（GSMxxx）。Sample实体必须仅引用一个Platform，可以包含在多个Series中。

GEO Sample (GSM) 样本ID号样品记录描述了处理单个样品的条件，它经历的操作以及从中得到的每个元素的丰度测量。每个样品记录都分配有唯一且稳定的GEO登录号（GSMxxx）。Sample实体必须仅引用一个Platform，可以包含在多个Series中。

GEO Series (GSE) study的ID号系列记录定义了一组被认为是组的一部分的相关样本，样本如何相关，以及它们是否以及如何排序。A系列提供了整个实验的焦点和描述。系列记录还可能包含描述提取数据，摘要结论或分析的表格。每个系列记录都分配有唯一且稳定的GEO登记号（GSExxx）。系列记录有几种格式，由GEOquery独立处理。较小的和新的GSEMatrix文件解析起来非常快; GEOquery使用一个简单的标志来选择使用GSEMatrix文件

GEO Dataset (GDS) 数据集的ID号 GEO DataSet（GDSxxx）是GEO样本数据的精选集合。GDS记录代表了一系列具有生物学和统计学意义的GEO样本，构成了GEO数据显示和分析工具套件的基础。GDS中的样本指的是同一个平台，也就是说，它们共享一组共同的探测元素。假设GDS中每个样本的值测量值以等效方式计算，即背景处理和标准化等考虑因素在整个数据集中是一致的。通过GDS子集提供反映实验设计的信息。

二、下载GEO数据

1. 用GEOquery包下载

#安装并加载GEOquery包**

**BiocManager::install('devtools')**

**library(GEOquery)**

**#在平台找好数据库，找到GSE编号，destdir设置当前目录，getGPL 和AnnotGPL都设置TRUE，可以下载和获得平台的注释文件**

**gse=getGEO("GSE29250",GSEMatrix = TRUE,destdir = ".",getGPL = T,AnnotGPL = T)**

下载完以后要查看数据集里的包含的实验信息、分组、数据处理等可以看这个翻译过的GEOquery包的说明书。https://www.jianshu.com/p/0d32fd410bcf

2. 网页直接下载

在NCBI上找到自己想要的数据集，可以在页面上看一下摘要、实验设计、数据分组。如果不是从文献中看到找的数据集，还可以看一下应用这个数据集的文献里做了哪些图，有哪些发现。

image

GSE提供的下载文件，包括SOFT文件，MINiML文件，Series Matrix File(s)以及原始数据文件，数据分析时，我们会用到Series Matrix File(s)文件或者原始数据。

image

具体的数据文件介绍先参考这个学习。

https://www.jianshu.com/p/74d570cb8c29

三、Illumina数据初步处理

这次先试了用R包直接下serise matrix中的数据集文件，打开看数据的时候发现数据里面有负数，还有很多1000多的极大值，应该是经过了各种数据预处理还有标准化之类的操作，但是看不明白作者处理数据用的方法。并且发现用的探针平台是Illumina的，Illumina出厂的表达芯片跟affymetrix系列芯片的分析方法不太一样。不懂! 看Jimmy大神对Illumina平台的介绍。www.biotrainee.com/thread-899-1-1.html

Illumina的原始数据有3个层级

GEO中给的原始数据一般来说是表中的第一个summary-level(c)数据，non bacground corrected, non normalized。然后给的是一个叫做sample probe profile 的数据，我们这找的是类似non-normalized.txt.gz名字的文件。

image

于是重新手动从网页下载好了non-normalized数据文件，按照Jimmy大神的Illumina数据处理教程，http://www.bio-info-trainee.com/1944.html 进行数据处理。

首先看一下我下载的数据，probeid ID, summarized expression level(AVG_Signal) 就是几个样本的表达量, a detection p-value(Dectection Pval) ，这些内容的含义在下图。

image

“对数据进行背景矫正以及标准化如分位数标准化的时候，一般还会有第二个文件control probe profile，一般这个文件GEO里边不会提供，因此文档中说也不是必须的。这个时候就需要用到这个dectection pval这一列来进行背景矫正以及标准化处理”

我的这个数据集有4个样本，但是下载得到的数据文件中每个样本的Dectection Pval的表头都是一样的，导致后面读取文件的一直报错，把这几个Dectection Pval表头分别加上其样品的后缀就能读取了。但是后面进行背景矫正以及标准化处理的时候也没给它指定哪个Dectection Pval对应哪个样本，不知道是不是直接按顺序来算的，反正是跑通了。

开始在R上加载文件分析数据

#安装并加载limma包，安装失败了好多次才装上，心累！

library(limma)

BiocManager::install('limma')

 #装不上包的时候多搜索一下，多试一下加一些参数，虽然不懂啥意思

BiocManager::install("scran", dependencies=TRUE, INSTALL_opts = c('--no-lock'))

#读取illumina的芯片结果文件，用R的lumi包来获取表达矩阵。

x <- read.ilmn("GSE40486-non-normalized.txt",expr="MUS",probeid="ID_REF",other.columns="Detection Pval")

y <- neqc(x,detection.p="Detection Pval") #背景矫正并且分位数标准化

y <- nec(x,detection.p="Detection Pval")#只有背景校正

expdata <- y$E

expdata

boxplot(data.frame(expdata),col="blue") ## 画箱式图，比较数据分布情况

exprSet=log2(expdata+1) #数据转化

exprSet

boxplot(data.frame(exprSet),col="blue") ## 画箱式图，比较数据分布情况

exprSet=normalizeBetweenArrays(exprSet)#数据标准化

boxplot(data.frame(exprSet),col="blue")## 画箱式图，比较数据分布情况

image

Jimmy大神说芯片数据质量控制结合了,N,T,B,Q(normalization,transformation,backgroud correction,qulity control)四个步骤，其中Q这个步骤又包括8种统计学方法。

这里用到了normalization ,transformation, backgroud correction, 这三个，还没有qulity control，不知道这几个步骤的顺序有没有要求。这里我用的顺序是试着backgroud correction，transformation，normalization。试着跑了一下backgroud correction，normalization，transformation的顺序，发现这样得到的箱线图就不在一条水平线上了，需要在重新一次normalization才能在一条水平线上。需要学习数据质量控制的统计学知识。

最后看到箱式图看起来几个样品的数据总体在一条水平线上。不知道有没有其他的判断数据是否处理好的指标。