TCGA差异分析|2.DESeq2

rm(list = ls())
library(DESeq2)
library(tidyverse)
library(pheatmap)
library(vioplot)
library(BiocParallel)
rt <- data.table::fread("F:/TCGA-counts/HTSeq_Counts_PRAD.txt",data.table = F) %>% 
  column_to_rownames("Tags")
metadata <- data.frame(TCGA_id=colnames(rt),
                       sample=factor(ifelse(str_sub(colnames(rt),14,15)<10,"tumor","normal"),levels = c("normal","tumor")))  #分组信息 

################################################################################
#1.核心环节：构建dds对象
###必备条件：countData, colData（分组信息）。
### 不是配对：design = ~grouplist
### 第一列如果有基因名称，tidy=TRUE
dds <- DESeqDataSetFromMatrix(countData = rt,
                              colData = metadata,
                              design = ~sample)
dds <- dds[rowSums(counts(dds))>1,] # 过滤原则：有原则是某基因在所有样本中的read counts之和小于样本总数（样本均read counts<1)则去除。
nrow(dds)
vsd <- vst(dds,blind = F) #数据标准化,去除批次
exprSet_vst <- as.data.frame(assay(vsd)) #assay提取vst标准化后的数据，用于表达量作图、热图等
plotPCA(vsd,"sample")
View(dds)
#2. 计算差异倍数及p值
dds <- DESeq(dds,parallel = T) 

#3.logFC矫正

contrast <- c("sample","tumor","normal")
dd1 <- results(dds,contrast,alpha = 0.05)
plotMA(dd1,ylim=c(-5,5))

dd2 <- lfcShrink(dds,contrast = contrast,res=dd1,type = "ashr")
plotMA(dd2,ylim=c(-5,5))

res <- dd2 %>% 
  as.data.frame() %>% 
  rownames_to_column("gene_id")

# 4.id转换
anno <- data.table::fread("F:/TCGA-anno/gencode.v22.annotation.gene.probeMap",data.table = F) %>% 
  select(id,gene)
res1 <- anno %>% 
  inner_join(res,by=c("id"="gene_id"))
write.table(res1,file = "diff-tcga-prad.txt",sep  = "\t",row.names = F,col.names = T,quote = F)
save(exprSet_vst,file = "vst-tcga-prad.Rdata")

logFC矫正还是很有作用的嘛

果子洲更老师强烈推荐的lfcShrink

声明：本文仅供学习交流用，禁止用于商业用途，如有侵权，请联系删除。
参考链接：
wx:喜大普奔，Deseq2重大升级，速度提升了60倍！(果子学生信）
wx:学好统计，我10分钟就完成了别人服务器通宵才完成的分析。(果子学生信）