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分钟就完成了别人服务器通宵才完成的分析。(果子学生信）
