DESeq2处理TCGA数据库Seq-count数据

1、DESeq2需要导入两个数据集：mycounts, colData。先说mycounts，这就是处理完的TCGA数据RNAmatrix.txt，直接读入即可。

library(tidyverse)
library(DESeq2)
#导入数据
setwd("E:/2.Hitseq_counts/")
mycounts<-read.table("RNAmatrix.txt", header=TRUE)
head(mycounts)
#这里有个x，需要去除，先把第一列当作行名来处理
rownames(mycounts)<-mycounts[,1]
#把带X的列删除
mycounts<-mycounts[,-1]
head(mycounts)

2、colData就是对每个样本的一个情况说明。这个可以生成，也可以自己写一个保存为csv格式。我一般自己写。

#载入colData文件
colData <-read.csv("RNAmatrix_condition.csv")
head(colData)

3、构建矩阵

#构建数据矩阵
dds <- DESeqDataSetFromMatrix(mycounts, colData, design= ~ condition)
dds <- DESeq(dds)
#查看dds的内容
dds
#接下来，我们要查看treat VS control的总体结果，并根据p-value进行重新排序。利用summary命令统计显示一共多少个genes上调和下调（FDR0.1）
res = results(dds, contrast=c("condition", "control", "treat")) ##或者res= results(dds)
res = res[order(res$pvalue),]
head(res)
summary(res)

4、输出结果

#所有结果先进行输出
write.csv(res,file="All_results.csv")
table(res$padj<0.05)
#获取padj（p值经过多重校验校正后的值）小于0.05，表达倍数取以2为对数后大于1或者小于-1的差异表达基因。代码如下
diff_gene_deseq2 <-subset(res, padj < 0.05 & abs(log2FoldChange) > 1.5) ##或者diff_gene_deseq2 <-subset(res,padj < 0.05 & (log2FoldChange > 1 | log2FoldChange < -1))
dim(diff_gene_deseq2)
head(diff_gene_deseq2)
write.csv(diff_gene_deseq2,file= "DEG_treat_vs_control.csv")