生信入门（二）——使用limma、Glimma和edgeR,RNA-seq数据分析

生信入门（二）——使用limma、Glimma和edgeR,RNA-seq数据分析

一、简介

简单且高效地分析RNA 测序数据的能力是生信的核心优势之一。通常在获得RNA-seq基因表达矩阵后，通常需要对数据进行数据预处理、探索性数据分析、差异表达检验以及通路分析，以得到可以帮助进一步试验和验证研究的结果。未来几篇将介绍使用edgeR包，limma包，Glimma包对RNA-seq数据展开分析，从原始计数（counts）中挖掘生物学意义

二、数据背景

在本篇文章中，描述了一个用于分析RNA-seq数据的edgeR - limma工作流程，使用基因水平的计数（gene-level counts）作为输入，经过预处理和探索性数据分析，然后得到差异表达（DE）基因和基因表达特征（gene signatures）的列表。Glimma包(Su et al. 2017)提供的交互式图表可以同时呈现整体样本层面与单个基因层面的数据，相对静态的R图表而言，更便于我们探索更多的细节。
此工作流程中我们分析的数据来自Sheridan等人的实验（2015）(Sheridan et al. 2015)，它包含三个细胞群，即基底（basal）、管腔祖细胞（liminal progenitor, LP）和成熟管腔（mature luminal, ML）。细胞群皆分选自雌性处女小鼠的乳腺，每种都设三个生物学重复。RNA样品分三个批次使用Illumina HiSeq 2000进行测序，得到长为100碱基对的单端序列片段。

三、初始配置

library(limma)
library(Glimma)
library(edgeR)
library(Mus.musculus)

四、数据整合

1、数据下载

数据来源：https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE63310&format=file 下载GSE63310_RAW.tar放在一定的文件夹

2、解压文件

setwd("D:\\RData\\prac002")#设置当前路径，下载数据在该文件夹下
utils::untar("GSE63310_RAW.tar", exdir = ".")
files<-c("GSM1545535_10_6_5_11.txt", "GSM1545536_9_6_5_11.txt", "GSM1545538_purep53.txt",
         "GSM1545539_JMS8-2.txt", "GSM1545540_JMS8-3.txt", "GSM1545541_JMS8-4.txt",
         "GSM1545542_JMS8-5.txt", "GSM1545544_JMS9-P7c.txt", "GSM1545545_JMS9-P8c.txt")
for(i in paste(files,".gz",sep=""))
  R.utils::gunzip(i,overwrite=TURE)

每一个文本文件均为对应样品的原始基因水平计数矩阵。需要注意我们的这次分析仅包含了此实验中的basal、LP和ML样品（可见下方所示文件名）。

files<-c("GSM1545535_10_6_5_11.txt", "GSM1545536_9_6_5_11.txt", 
         "GSM1545538_purep53.txt", "GSM1545539_JMS8-2.txt", 
         "GSM1545540_JMS8-3.txt", "GSM1545541_JMS8-4.txt", 
         "GSM1545542_JMS8-5.txt", "GSM1545544_JMS9-P7c.txt", 
         "GSM1545545_JMS9-P8c.txt")
read.delim(files[1],nrow=5)#查看导入数据

运行结果

相比于分别读入这九个文本文件然后合并为一个计数矩阵，edgeR提供了更方便的途径，使用readDGE函数即可一步完成。得到的DGEList对象中包含一个计数矩阵，它的27179行分别对应每个基因不重复的Entrez基因ID，九列分别对应此实验中的每个样品。

class(x)
dim(x)

运行结果

注：如果数据不是每个样本一个文件的形式，而是包含所有样本的技术的文件，则可以先将文件读入R，再使用DGEList函数转换为一个DGEList对象

3、组织样本信息

为进行下游分析，需要将有关实验设计的样品信息与计数矩阵的列关联起来。这里需要包括各种对表达水平有影响的实验变量，无论是生物变量还是技术变量。例如，细胞类型（在这个实验中是basal、LP和ML）、基因型（野生型、敲除）、表型（疾病状态、性别、年龄）、样品处理（用药、对照）和批次信息（如果样品是在不同时间点进行收集和分析的，需要记录进行实验的时间）等。
我们的DGEList对象中包含的samples数据框同时存储了细胞类型（group）和批次（测序泳道lane）信息，每种信息都包含三个不同的水平。在x$samples中，程序会自动计算每个样品的文库大小（即样品的总序列计数），归一化系数会被预先设置为1。 为了方便阅读，我们从DGEList对象x的列名中删去了GEO样品ID（GSM*）。

samplenames<-substring(colnames(x),12,nchar(colnames(x)))
colnames(x)<-samplenames
group<-as.factor(c("LP", "ML", "Basal", "Basal",