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.核心環(huán)節(jié)：構(gòu)建dds對(duì)象
###必備條件：countData, colData（分組信息）席里。
### 不是配對(duì)：design = ~grouplist
### 第一列如果有基因名稱，tidy=TRUE
dds <- DESeqDataSetFromMatrix(countData = rt,
colData = metadata,
design = ~sample)
dds <- dds[rowSums(counts(dds))>1,] # 過濾原則：有原則是某基因在所有樣本中的read counts之和小于樣本總數(shù)（樣本均read counts<1)則去除拢驾。
nrow(dds)
vsd <- vst(dds,blind = F) #數(shù)據(jù)標(biāo)準(zhǔn)化,去除批次
exprSet_vst <- as.data.frame(assay(vsd)) #assay提取vst標(biāo)準(zhǔn)化后的數(shù)據(jù)奖磁，用于表達(dá)量作圖、熱圖等
plotPCA(vsd,"sample")
View(dds)
#2. 計(jì)算差異倍數(shù)及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轉(zhuǎn)換
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矯正還是很有作用的嘛
果子洲更老師強(qiáng)烈推薦的lfcShrink
聲明：本文僅供學(xué)習(xí)交流用繁疤，禁止用于商業(yè)用途署穗，如有侵權(quán)，請(qǐng)聯(lián)系刪除嵌洼。
參考鏈接：
wx:喜大普奔，Deseq2重大升級(jí)封恰，速度提升了60倍麻养！(果子學(xué)生信）
wx:學(xué)好統(tǒng)計(jì)，我10分鐘就完成了別人服務(wù)器通宵才完成的分析诺舔。(果子學(xué)生信）
