RNA-Seq数据使用DESeq2、edgeR和limma-voom进行差异分析

###DESeq2  DEseq2针对有生物学重复的样本。
###数据来源https://4va.github.io/biodatasci/r-rnaseq-airway.html#data_needed  

metadata <-  read.csv("airway_metadata.csv")
library(DESeq2)
dds <- DESeqDataSetFromMatrix(countData=counts, 
                              colData=metadata, 
                              design=~dex, 
                              tidy=TRUE)
###其中countData存放counts数据，colData存放样本信息的数据，design就是实验设计。
#####数据预处理
dds1 <- dds [ rowSums(counts(dds)) > 1, ]

##方差齐性转换，rlog更适用于样本量小于30的数据，相应的vst负责的是样本量大于30的数据
rld <- rlog(dds, blind = FALSE)
head(assay(rld), 3)
vsd <- vst(dds, blind = FALSE)
head(assay(vsd), 3)

sampleDists <- dist(t(assay(rld)))##样本距离，查看样本之间的相似度和差异程度
library("pheatmap")
library("RColorBrewer")#自动生成一组颜色
sampleDistMatrix <- as.matrix( sampleDists )

colors <- colorRampPalette( rev(brewer.pal(8, "Blues")) )(255)
pheatmap(sampleDistMatrix,
         clustering_distance_rows = sampleDists,
         clustering_distance_cols = sampleDists,
         col = colors)#聚类的作用是判断样本之间的相似度

##PCA分析
#还有一种可视化样本-样本距离的方法就是主成分分析。
##DESeq2提供了专门的方法用于作图
plotPCA(rld,intgroup=c("dex"))
 

########差异分析
dds <- DESeq(dds)#包含了以下三个函数在内
#dds <- estimateSizeFactors(dds) 
#计算归一化系数sizeFactor
##dds <- estimateDispersions(dds) 
##估计基因的离散程度
##dds <- nbinomWaldTest(dds) 
##统计检验，差异分析
##获得分析结果
res1 <- results(dds,tidy=TRUE)#tidy为TRUE是数据框，否则为S4对象

metadata$dex=c("a","a","a","b","b","c","c","c")#重新定义分组,3个组

dds <- DESeqDataSetFromMatrix(countData=counts,