手动从GEO下载预处理Matrix文件并进行简单分析 | LiuMWei’s Recordings

手动从GEO下载预处理Matrix文件并进行简单分析

本篇将具体讲解和演示如何从GEO手动下载相关数据,并利用limma和ggplot2的R包进行简单的差异基因分析。

编者:龙飞、汪宣佚、柴子轩;补充校对:liumwei

1. 首先搭建R环境

Windows下的R环境,请参考这个页面进行。若想重现本文相关操作和信息输出结果,请参看该页面进行R版本安装与配置进行。同时,本文所涉及的GEOquery和limma包安装,均借助Bioconductor下的BiocManager来进行,因此,请务必确认您已进行了如下操作:

>if (!requireNamespace(“BiocManager”, quietly = TRUE))
    install.packages(“BiocManager”)

如果上述执行有错误,请回到这个页面查看解决办法。

2. limma和ggplot2包的安装

在RStudio Console界面的>提示符后输入如下命令,并回车执行:

>BiocManager::install(“limma”)
>BiocManager::install(“ggplot2”)

提示:命令输入后,记得回车执行;中文输入法状态下,可能会出现不能复制粘贴命令(最好采用shift+ctrl+v粘贴),或手动不能输入命令。此时切换英文输入法状态即可。如果出现包下载速度很慢的状态,请回到这个页面查看解决办法。同时,安装ggplot2,可能会出现“Installation path not writeable, unable to update packages: survival”的错误,此时需要使用管理员权限打开RStudio,然后安装survival:

>BiocManager::install()
>BiocManager::install(‘survival’)

如果安装limma也出现类似的错误,解决方式同上。

3. 加载相应的R包

>library(‘Biobase’)
>library(‘limma’)​​​​
>library(‘ggplot2’)​​​​

提示:* Biobase R包具有exprs()、pData()和fData()等函数,可利用它获取基因芯片上的表达数据、样品处理分组等信息以及芯片平台的设计注释信息;
          * limma包用于差异相关分析。
          * ggplot2用于绘制火山图。

4. 手动下载GEO matrixfile载并初始化数据(示例性数据GSE70213

首先打开GEO站点(https://www.ncbi.nlm.nih.gov/geo/),在搜索框中输入GSE70213,点击搜索,可获得GSE70213的记录页面。此时,找到页面底下的“Series Matrix File(s)”,点击下载,将数据压缩包保存到事先设置好的R用户目录下(如我设置的为C:UsersliumweiDocumentsR),然后解压可获得GSE70213_series_matrix.txt文件,如下图:

同时,在数据下载时,应对该数据的实验设计等信息充分了解(了解实验目的和设计很重要!),如下图:

可见,该数据集来自小鼠中的一种丝状蛋白Nebulin(与骨骼肌肌动蛋白丝共存的一种蛋白)研究。 尽管目前,先知Nebulin突变是导致Nemaline肌病(NEM)的主要原因,并在典型NEM成年患者中Nebulin低表达。不过,Nebulin在成年肌肉中的作用仍然知之甚少。 故,该研究目的是确定成年肌肉神经元蛋白的功能,模型上采用一种新型的条件性神经元蛋白KO(Neb cKO)雄性小鼠,实验设计上采用两组对照-疾病组设计,即quad-cKO vs quad-control 以及soleus-cKO vs soleus-control,且每个组分别6个重复。

在后面的表型数据分组时,cKo与control可作为表型标签进行提取。

5. 表达值矩阵获取生成

采用Excle、Notepad+++这样的工具打开GSE70213_series_matrix.txt,然后将文件中的首位中含有!的行全部删除,只留下从”ID_REF”开始的数据行(如下图),然后另存为GSE70213_matirx.

txt

将GSE70213_matirx.txt用Excel打开,数据形式如下:

然后另存为GSE70213_matrix.csv(也就是保存为csv格式。该操作务必进行!如果打开excle打开该文件,发现数据集中在一列里,还需进行数据分列操作。)。

6. 将上面生成的csv数据读入(RStudio Console中操作),然后进行log2转换

>GSE70213_expre <-read.csv(GSE70213_matrix.csv’, ,header=TRUE,row.names=1) #或者采用下面的命令读入亦可
>show(GSE70213_expre) #可以查看读入的数据结果究竟是怎么样的。后面的每个命令操作过程,均可用这个命令进行查验。这也是初学者必备的素质。
>GSE70213_expre<-read.table(‘GSE70213_matrix.csv’,header=TRUE,row.names=1) #需要注意
>GSE70213_log2 <- log2(GSE70213_expre) #将表达值进行对数转换

提示:对表达矩阵值进行对数转换,是否需要,其实可以先肉眼观察一下数据。一般数值中出现较大数据和过小数据,最好先进行转换。当然,来自GEO2R有一个自动判断脚本,经修改后如下:

qx <- as.numeric(quantile(GSE70213_expre, c(0., 0.25, 0.5, 0.75, 0.99, 1.0), na.rm=T))
LogC <- (qx[5] > 100) ||
               (qx[6]-qx[1] > 50 && qx[2] > 0) ||
              (qx[2] > 0 && qx[2] < 1 && qx[4] > 1 && qx[4] < 2)
if (LogC) { GSE70213_expre[which(GSE70213_expre <= 0)] <- NaN
exprSet <- log2(GSE70213_expre)
print(“log2 transform needed”)}
           else{print(“log2 transform not needed”)}

还需要注意的是,如果表达矩阵中出现负数,在进行log转换时,需要提前处理:要么舍弃,要么加上一个固定值(这个一般有商业公司的芯片数据处理手册提供指南)转化成正数。或者取负数的绝对值进行log转换,然后再把负号加上。相关的代码如下:

GSE70213_matrix_abs <- function(x){
          x[x==0] <-1
           si <- sign(x)
          si * log2(si*x)} 
GSE70213_expres_abs<-GSE70213_matrix_abs(GSE70213_expre)

7. 基于箱线图进行数据整体情况调查

>boxplot(data.frame(GSE70213_log2),col=”red”)

结果如下:

从箱线图可以看出,样本间的基因表达数据一致性还是很高。比如中位值均在5附近,最高值和最低值基本接近。不过,其他有些数据整齐可能很不好,甚至会出现下图的结果。尽管如此,我们看到数据整齐性还可以在优化一下。优化的最好办法是对数据进行标准化操作(至于为何要标准化操作及背后的原理,可查看这个页面的解释说明。),再进入下一步基因差异数据分析:

标准化是组学数据最基本的问题,有多种方法可以进行。在此,我们采用limma包中normalizeBetweenArrays函数来重新标准化,相关命令和结果如下:

GSE70213_expr_N <- normalizeBetweenArrays(GSE70213_expre)
GSE70213_N_log2 <- log2(GSE70213_expr_N)
boxplot(data.frame(GSE70213_N_log2),col=”red”)

可见,重新标准化的数据,其整齐度更好。

8. 基因差异分析

A. 数据分组准备

根据GSE70213的页面样本信息(表型信息)对数据进行手动分组准备(下图),保存为csv格式(GSE70213_type.csv)

准备好上面的表型分组数据后,在RStudio中将其载入到内存中:

>GSE70213_phenodata <- read.table(‘GSE70213_type.csv’, header=TRUE, row.names=1, check.names = FALSE)

提示:* 采用read.table()将表型分类数据表读入内存中。 header=TRUE表示要检查表头;row.names=1表示把第一列设为行名。
          * 注意,如果出现如下错误,这表明在”row.names=1″前输入了不能识别的符号(比如中文的空格)。为避免这个问题,请确认输入命令时,请切记将输入法切换到英文状态,并使用删除键,试着将”row.names=1″可能的空格删除。这是初学者常常会出现的错误和问题。在前面的命令输入时,也特别注意这个问题。

错误: unexpected input in "GSE70213_phenodata <- read.csv('GSE70213_type.csv', header=TRUE,"

B. 构建分组矩阵

>GSE70213_group <- GSE70213_phenodata [,1] #获取第二列的表型标签
>GSE70213_design <- model.matrix(~0+factor(GSE70213_group)) #将表型标签进行数字0和1的转换
>head(GSE70213_design) #查看表型标签进行数字转换后,整个表型表变换的结果,它和show()的结果是异曲同工的
>colnames(GSE70213_design) <- levels(factor(GSE70213_group)) #将表型标签作为列
>rownames(GSE70213_design) <- colnames(GSE70213_N_log2) #将对数表达值作为行

提示:colnames和rownames分别表示设置分组举证的列和行

C. 基因差异矩阵数据

>GSE70213_cont.matrix <- makeContrasts(paste0(unique(GSE70213_group),collapse=”-“),levels=GSE70213_design)

提示:makeContrasts()构建比较矩阵

D. 差异分析

>GSE70213_fit <- lmFit(GSE70213_N_log2,GSE70213_design) #给定系列数据,为每个基因进行线性模型拟合)
>GSE70213_fit2=contrasts.fit(GSE70213_fit,GSE70213_cont.matrix) #基于拟合线性模型进行系数估计和标准误差计算
>GSE70213_fit3 <- eBayes(GSE70213_fit2)#基于经验贝叶斯进行进行t统计量、F统计量和差异对数表达计算
>GSE70213_output = topTable(GSE70213_fit3,coef=1,n=Inf,adjust=”BH”) #从线性模型拟合结果中,进行排序,并提取排序靠前的基因列表
>GSE70213_nrDEG = na.omit(GSE70213_output) #删除基因排序列表中出现的缺失值,并将结果返回

提示R中通过na.fail和na.omit来处理样本中的缺失值。na.fail()用于先判断,后返回结果:如果样本中有缺失值,回返回错误,如果没有,则返回结果。而na.omit()则表示将缺失值删除后,把结果返回。attr(na.omit(),”na.action”)表示对缺失值进行检查,并返回缺失值的下标;is.na():判断向量内的元素是否有缺失(NA)

F.   差异分析结果文件

>write.table (GSE70213_nrDEG, paste(“GSE70213_DE.csv”)) #用下面的命令,也能将结果存到文档中
>write.table(GSE70213_nrDEG, “GSE70213_DE.csv”, quote=F, sep=”,”) #该命令加入了相关文档数字分隔的字符定义

提示:csv 文件一般是用逗号分隔的,故而 sep = “,”;tsv 文件是用制表符分隔的,故而 sep = “t”;常用的分隔符还有空格 sep = ” “

9. 差异基因火山图可视化结果

>GSE70213_alldiff<-read.csv(“GSE70213_nrDEG.csv”,header=TRUE,row.names=1,check.names=FALSE)
>GSE70213_alldiff$threshold[GSE70213_alldiff$P.Value < 0.05 & GSE70213_alldiff$logFC > 0.5] = “up”
>GSE70213_alldiff$threshold[GSE70213_alldiff$P.Value < 0.05 & GSE70213_alldiff$logFC < -0.5] = “down”
>GSE70213_alldiff$threshold[GSE70213_alldiff$P.Value > 0.05 | (GSE70213_alldiff$logFC >= -0.5 & GSE70213_alldiff$logFC <= 0.5)] = “non”
>ggplot(GSE70213_alldiff,aes(x=logFC,y=-log10(P.Value),colour=threshold))+ xlab(“log2 Fold Change”)+ ylab(“-log10 P.Value”)+ geom_point()+ scale_color_manual(values =c(“#6495ED”,”red”,”#FF4500″))+ geom_hline(yintercept=1.30103,linetype=3)+ geom_vline(xintercept=c(-0.5,0.5),linetype=3)+ theme(legend.title = element_blank(),panel.grid.major = element_blank(),panel.grid.minor = element_blank())

提示:GSE70213_alldiff$threshold用于设置取值“up”、“down”和  “non”;ggplot()绘图函数,其中P.Value和logFC的值需根据研究具体情况而定。这里,分别被设置为0.05和0.5,aes()为坐标轴取值设定,xlab和ylab为横纵坐标名称设定,scsle_color_manual为设定颜色值,geom_hline和vline为分割线设定,theme为主题设定。

Measure
Measure

点击数:0

发表评论

邮箱地址不会被公开。 必填项已用*标注