RNA-seq分析（三）DESeq2

1、DESeq2差异表达分析

R语言数据类型：向量（vector）；列表（list）；矩阵（matrix）；数组（array）；因子（factor)；数据框（data.frame)
详细介绍可参考 R 数据类型 | 菜鸟教程 (runoob.com)
R语言数据的导入
read.table()，从带分隔符"\t"的文本文件中读入，返回的是数据框。
col.names 指定列名的向量，c() 是一个创造向量的函数。
在shell下写R语言脚本 vim DESeq2.R ；运行脚本 Rscript DESeq2.R。
或者进入R，分别执行每行的命令

#!bin/bash
library(DESeq2)
control1 <- read.table("SRR7059707.count", sep = "\t", col.names = c("gene_id","BY1"))
control2 <- read.table("SRR7059706.count", sep = "\t", col.names = c("gene_id","BY2"))
control3 <- read.table("SRR7059709.count", sep = "\t", col.names = c("gene_id","BY3"))
SY_1<- read.table("SRR7059708.count", sep = "\t", col.names = c("gene_id","SY1"))
SY_2<- read.table("SRR7059705.count", sep = "\t", col.names = c("gene_id","SY2"))
SY_3<- read.table("SRR7059704.count", sep = "\t", col.names = c("gene_id","SY3"))
raw_count1 <- merge(merge(control1, control2, by="gene_id"), control3, by="gene_id")
raw_count2 <- merge(merge(SY_1, SY_2, by="gene_id"), SY_3, by="gene_id")
raw_count <- merge(raw_count1, raw_count2, by="gene_id")
#merge 一次只能合并两个数据集
head(raw_count)
tail(raw_count)
raw_count_filt <- raw_count[-1:-5,]
#删除前五行比对信息
head(raw_count_filt)
tail(raw_count_filt)
#查看是否删除前五行
countData <- raw_count_filt
rownames(countData) <- countData[,1]
countData <- countData[,-1]
head(countData)
#去除gene_id，把第一列当做行名来处理，删除第一列列名
condition <- factor(c(rep("BY",3),rep("SY",3)))
colData <- data.frame(row.names=colnames(countData), condition)
dds <- DESeqDataSetFromMatrix(countData, colData, design= ~condition)
#构建dds对象，design= ~condition 只有一个变量
head(dds)
dds
dds$condition
#查看因子水平
dds <- DESeq(dds)
#使用DESeq()函数差异分析
res <- results(dds)
#用results()函数来提取分析结果
summary(res)
#summary看一下结果的概要信息
res
#查看treatment与control谁在前 ，这里是SY vs BY
table(res$padj<0.001)
DGE <- rep("NC", nrow(res))
head(DGE)
DGE[((res$padj) < 0.001) & (res$log2FoldChange >= 1)] = "Up"
DGE[((res$padj) < 0.001) & (res$log2FoldChange < -1)] = "Down"
rTable = data.frame(res[, c("baseMean","log2FoldChange", "pvalue", "padj")], DGE)
#新建一个数据框，取diffgene_R矩阵的所有行，"baseMean","log2FoldChange", "pvalue", "padj"列，以及DGE列
Up=rTable[which(rTable$DGE=="Up"),]
Down=rTable[which(rTable$DGE=="Down"),]
dim(Up)
dim(Down)
#查看Up、Down行数
save(rTable,Up,Down,dds,file="diff.RData")
#保存变量到diff.RData文件，使用时，load("diff.RData")即可，不用每次都保存.RData，否则拖慢R运行速度。
write.csv(rTable,file= "SY14_VSBY4742.csv")
write.csv(Up,file= "SY14_up.csv")
write.csv(Down,file= "SY14_down.csv")

导出SY14_VSBY4742.csv所有基因的表格，可用于GSEA差异分析
导出SY14_up.csv，可用于GO、KEGG通路分析。

2、PCA 主成分分析

DESeq2提供两种数据标准化方法：VST (variance stabilizing transformations)和rlog (regularized logarithm)。
两种方法都使用了log2缩放，并且已经进行了library size 或其他normalization factors的校正，可能使用VST标准化后数据的标准差比较稳定。

library(DESeq2)
library(ggplot2)
load("diff.RData")
vsd <- vst(dds, blind=TRUE)
head(assay(vsd), 6)
pcaData <- plotPCA(vsd, intgroup="condition", returnData=TRUE)
pcaData
#pcaData矩阵如下
          PC1        PC2 group condition name
BY1 -5.698124  0.3943525    BY        BY  BY1
BY2 -6.060169 -1.1671873    BY        BY  BY2
BY3 -6.200917  0.8404256    BY        BY  BY3
SY1  5.794393 -2.9737566    SY        SY  SY1
SY2  5.928922  1.6493463    SY        SY  SY2
SY3  6.235895  1.2568194    SY        SY  SY3
pdf("PCA_2.pdf",width = 5, height = 4) 
#输出pdf
ppi<- 300
png("PCA_2.png",width=5*ppi,height=4*ppi, res=ppi)
#或者直接输出png图片
percentVar <- round(100 * attr(pcaData, "percentVar"))
percentVar
xlab=paste0("PC1 (",percentVar[1],"%)") 
ylab=paste0("PC2 (",percentVar[2],"%)")
ggplot(pcaData, aes(PC1, PC2, color=condition)) +
geom_point(size=1.5) +theme_bw() +
labs(x= xlab, y= ylab,title= "PCA")+
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
theme(plot.title = element_text(hjust = 0.5)) + theme(legend.title=element_blank())
#去除网格线
dev.off()

从PC1主成分上看，样本聚类挺好的。