我们之前介绍了limma包，limma包是对基因芯片表达矩阵的分析，不能对逆转录RNAseq表达矩阵进行分析（因为数据特征不同），RNAseq需要用另一种方法：DESeq2

（注：基因芯片和RNAseq是测定表达量的两种方式，各有优劣，详细可自行百度）

一.读取数据

library(airway)
#Biocductor R包为三种：1.功能函数包2.数据包3.注释包（芯片基因之间的转换）
#此为中的一种，为数据包
data(airway)#加载数据
exprSet=assay(airway)#获取表达矩阵，默认airway获取表达矩阵就是assay，没有原因的
colnames(exprSet)#看表达矩阵的列名
dim(exprSet)#查看表达矩阵的维度

二.设定分组信息

group_list=colData(airway)[,3]#得出分组信息
tmp=data.frame(group_list)#把group_list向量变为数据框tmp
row.names(tmp)=colnames(exprSet)
#把tmp的行名改为exprSet的列名

三.筛选有意义的基因

筛选基因之后进行相关性计算，去掉基因中：有3个以上为表达量0的样本的基因

exprSet=exprSet[apply(exprSet,1,function(x)sum(x>1)>5),]
##分别对数据中每一行的数据进行一个什么运算，1代表行，2代表列

四.DESeq2

#DESeq2
library("DESeq2") 
colData<-data.frame(row.names = colnames(exprSet),group_list=group_list)
#给每个样本进行标签化，设定每个样本的性质特点，分组的前期准备

colData

dds <- DESeqDataSetFromMatrix(countData =exprSet, colData = colData, design = ~group_list)
#countData为表达矩阵，colData样本特点内涵分组信息，design 以某个种（group_list）样本特点设定分组
#keep <- rowSums(counts(dds)) >= 10#筛选表达量之和大于10的基因
#dds <- dds[keep,]
dds <- DESeq(dds)#四步矫正，处理数据
res <- results(dds, contrast=c("group_list","untrt","trt"))#按照"untrt","trt"比较得出差异分析结果
resOrdered<-res[order(res$padj),]#把res排序
head(resOrdered)
DEG=as.data.frame(resOrdered)#把差异分析结果变为数据框
DESeq2_DEG=na.omit(DEG)#删除差异分析中缺少值的结果

QQ截图20190401215947.jpg

得出结果后就可以进一步按照之前的方法富集分析等，原理是一样的

最后

感谢jimmy的生信技能树团队！

感谢导师岑洪老师！

感谢健明、孙小洁，慧美等生信技能树团队的老师一路以来的指导和鼓励！

特别注明:此文中编码来自生信技能树健明老师