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

library(tidyverse)

library(DESeq2)

#导入数据

setwd("E:/2.Hitseq_counts/")

mycounts

head(mycounts)

#这里有个x，需要去除，先把第一列当作行名来处理

rownames(mycounts)

#把带X的列删除

mycounts

head(mycounts)

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

#载入colData文件

colData

head(colData)

3、构建矩阵

#构建数据矩阵

dds

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 1.5) ##或者diff_gene_deseq2 1 | log2FoldChange < -1))

dim(diff_gene_deseq2)

head(diff_gene_deseq2)

write.csv(diff_gene_deseq2,file= "DEG_treat_vs_control.csv")