DEseq標準流程

step1 處理group_list

library(DESeq2)

colData <- data.frame(row.names=colnames(exprSet),group_list=group_list)

step2 處理數(shù)據(jù)構(gòu)建dds:需要exprSet colData 和group_list

# 按照DESeq2方法制作dds對象
dds <- DESeqDataSetFromMatrix(countData = exprSet,colData = colData,
                              design = ~ group_list)
dds_2 <- DESeq(dds)
save(dds_2,group_list, file = "4.DESeq_analysis.Rdata")

step3.處理結(jié)束，拿數(shù)據(jù)

rm(list = ls())
load(file = "4.DESeq_analysis.Rdata")
dds=dds_2
# 提取你想要的差異分析結(jié)果，我們這里是trt組對untrt組進行比較
res <- results(dds,contrast=c("group_list","trt","untrt"))
resOrdered <- res[order(res$padj),]

head(resOrdered)
DEG =as.data.frame(resOrdered)
DESeq2_DEG = na.omit(DEG)
## ttmpp=DEG[!(rownames(DEG) %in% rownames(DESeq2_DEG)),]
## 可以看出過濾na的全是padj=NA的行，一共是2302行

## 取出logFC和P值
nrDEG=DESeq2_DEG[,c(2,6)]
colnames(nrDEG)=c('log2FoldChange','pvalue')  
save(nrDEG, DESeq2_DEG, file = "4.DESeq.Rdata")

下面是我自己的數(shù)據(jù)實戰(zhàn)，但是我的表達矩陣是已經(jīng)被均一化的，所以dds的時候報錯ome values in assay are not integers。我的數(shù)字不是整數(shù)型。。最后我只做了PCA分析。

rm(list = ls())
group<- read.table(file = 'condition.txt', header = T)

### 第一步構(gòu)建dds對象
library(DESeq2)
#########step1 處理group_list  
colData <- as.data.frame(group[,-1])
row.names(colData)<- group[,1]
colnames(colData)<- c('group')
group_list<- group[,-1]

#########step2 處理數(shù)據(jù)

library(DESeq2)
library(dplyr)
expr<- read.table(file = 'ALL_summary_nrpm.txt', header = T)
exprselect<- select(expr,'Gene',c(rownames(colData)))  ###################只選擇了有分組信息的數(shù)據(jù)
row.names(exprselect)<- exprselect[,1] 
exprselect<- exprselect[,-1]
dds <- DESeqDataSetFromMatrix(countData = exprselect,
                              colData = colData,
                              design = ~ group_list)############################some values in assay are not integers

下面是簡單的PCA分析

需要準備成這樣的：需要有一個名為dat的表達矩陣，還需要一個group_list

image.png

image.png

# PCA 每次都要檢測數(shù)據(jù)
dat=exprselect
dat[1:4,1:4]
dim(dat)

################################step1:準備數(shù)據(jù)
## 下面是畫PCA的必須操作，需要看說明書。
dat=t(dat)
dat=as.data.frame(dat)
dat=cbind(dat,group_list)
# dat$group_list=group_list


################################step2:畫圖
library("FactoMineR")
library("factoextra") 
# The variable group_list (index = 54676) is removed
# before PCA analysis
ncol(dat)
dat[,ncol(dat)]
# 畫圖是用數(shù)據(jù)，因此去掉最后一列grouplist
dat.pca <- PCA(dat[,-ncol(dat)], graph = FALSE)
plot(dat.pca,choix="ind")
# Graph of individuals
fviz_pca_ind(dat.pca,
             geom.ind = "point", # show points only (nbut not "text")
             col.ind = dat$group_list, # color by groups
             palette = c("#00AFBB", "#E7B800", '#CC00FF', '#FF0099'),
             addEllipses = TRUE, # Concentration ellipses # 是否圈起來
             legend.title = "Groups"
)
ggsave('7.PCA_all_samples.png')

展示最后的結(jié)果：結(jié)果這么慘烈，我還發(fā)不發(fā)文章了