三种方法任选1个，或者三个全选
##deseq2处理（DEG为标注化后的exp，并且添加了change列，用来质控数据）

  1. rm(list = ls())
  2. load("TCGA-HNSC.Rdata")
  3. table(Group)
  4. #deseq2----
  5. library(DESeq2)
  6. colData <- data.frame(row.names =colnames(exp), 
  7.                       condition=Group)
  8. if(!file.exists(paste0(proj,"_dd.Rdata"))){
  9.   dds <- DESeqDataSetFromMatrix(
  10.     countData = exp,
  11.     colData = colData,
  12.     design = ~ condition)
  13.   dds <- DESeq(dds)
  14.   save(dds,file = paste0(proj,"_dd.Rdata"))
  15. }
  16. load(file = paste0(proj,"_dd.Rdata"))
  17. class(dds)
  18. res <- results(dds, contrast = c("condition",rev(levels(Group))))
  19. c("condition",rev(levels(Group)))
  20. class(res)
  21. DEG1 <- as.data.frame(res)
  22. DEG1 <- DEG1[order(DEG1$pvalue),] 
  23. DEG1 = na.omit(DEG1)
  24. head(DEG1)
  25. #添加change列标记基因上调下调
  26. logFC_t = 2
  27. pvalue_t = 0.05
  28. k1 = (DEG1$pvalue < pvalue_t)&(DEG1$log2FoldChange < -logFC_t);table(k1)
  29. k2 = (DEG1$pvalue < pvalue_t)&(DEG1$log2FoldChange > logFC_t);table(k2)
  30. DEG1$change = ifelse(k1,"DOWN",ifelse(k2,"UP","NOT"))
  31. table(DEG1$change)
  32. head(DEG1)

edgeR处理

  1. library(edgeR)
  2. dge <- DGEList(counts=exp,group=Group)
  3. dge$samples$lib.size <- colSums(dge$counts)
  4. dge <- calcNormFactors(dge) 
  6. design <- model.matrix(~Group)
  8. dge <- estimateGLMCommonDisp(dge, design)
  9. dge <- estimateGLMTrendedDisp(dge, design)
  10. dge <- estimateGLMTagwiseDisp(dge, design)
  12. fit <- glmFit(dge, design)
  13. fit <- glmLRT(fit) 
  15. DEG2=topTags(fit, n=Inf)
  16. class(DEG2)
  18. DEG2=as.data.frame(DEG2)
  19. head(DEG2)
  20. k1 = (DEG2$PValue < pvalue_t)&(DEG2$logFC < -logFC_t);table(k1)
  22. k2 = (DEG2$PValue < pvalue_t)&(DEG2$logFC > logFC_t);table(k2)
  23. DEG2$change = ifelse(k1,"DOWN",ifelse(k2,"UP","NOT"))
  25. head(DEG2)
  26. table(DEG2$change)

limma处理

  1. library(limma)
  2. dge <- DGEList(counts=exp)
  3. dge <- calcNormFactors(dge)
  4. v <- voom(dge,design, normalize="quantile")
  6. design <- model.matrix(~Group)
  7. fit <- lmFit(v, design)
  8. fit= eBayes(fit)
  10. DEG3 = topTable(fit, coef=2, n=Inf)
  11. DEG3 = na.omit(DEG3)
  13. k1 = (DEG3$P.Value < pvalue_t)&(DEG3$logFC < -logFC_t);table(k1)
  15. k2 = (DEG3$P.Value < pvalue_t)&(DEG3$logFC > logFC_t);table(k2)
  17. DEG3$change = ifelse(k1,"DOWN",ifelse(k2,"UP","NOT"))
  18. table(DEG3$change)
  20. head(DEG3)

组合三种数据并保存

  1. tj = data.frame(deseq2 = as.integer(table(DEG1$change)),
  2.                 edgeR = as.integer(table(DEG2$change)),
  3.                 limma_voom = as.integer(table(DEG3$change)),
  4.                 row.names = c("down","not","up")
  5. );tj
  7. save(DEG1,DEG2,DEG3,Group,tj,file = paste0(proj,"_DEG.Rdata"))