TCGA下载数据，并注释

    1. #下载数据
    2. proj = "TCGA-HNSC"
    3. download.file(url = paste0("https://gdc.xenahubs.net/download/",proj, ".htseq_counts.tsv.gz"),destfile = paste0(proj,".htseq_counts.tsv.gz"))
    4. #表达矩阵行名ID转换
    5. dat = read.table(paste0(proj,".htseq_counts.tsv.gz"),check.names = F,row.names = 1,header = T)
    6. ##逆转log
    7. dat = as.matrix(2^dat - 1)
    8. dat[1:4,1:4]
    9. ### 深坑一个
    10. dat[97,9]
    11. as.character(dat[97,9])
    12. ### 用apply转换为整数矩阵
    13. exp = apply(dat, 2, as.integer)
    14. exp[1:4,1:4] #行名消失
    15. rownames(exp) = rownames(dat)
    16. exp[1:4,1:4]
    18. #表达矩阵过滤
    19. nrow(exp)
    20. ##常用过滤标准1：
    21. ###仅去除在所有样本里表达量都为零的基因
    22. exp1 = exp[rowSums(exp)>0,]
    23. nrow(exp1)
    24. ##常用过滤标准2(推荐)：
    25. ###仅保留在一半以上样本里表达的基因
    26. exp = exp[apply(exp, 1, function(x) sum(x > 0) > 0.5*ncol(exp)), ]
    27. nrow(exp)
    28. #分组信息获取
    29. ##根据样本ID的第14-15位，给样本分组（tumor和normal）
    30. library(dplyr)
    31. library(tidyverse)
    32. table(str_sub(colnames(exp),14,15))
    33. ###分组
    34. Group = ifelse(as.numeric(str_sub(colnames(exp),14,15)) < 10,'tumor','normal')
    35. Group = factor(Group,levels = c("normal","tumor"))
    36. table(Group)
    37. save(exp,Group,proj,file = paste0(proj,".Rdata"))

    标准化数据

    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)
    35. #edgeR----
    36. library(edgeR)
    37. dge <- DGEList(counts=exp,group=Group)
    38. dge$samples$lib.size <- colSums(dge$counts)
    39. dge <- calcNormFactors(dge) 
    41. design <- model.matrix(~Group)
    43. dge <- estimateGLMCommonDisp(dge, design)
    44. dge <- estimateGLMTrendedDisp(dge, design)
    45. dge <- estimateGLMTagwiseDisp(dge, design)
    47. fit <- glmFit(dge, design)
    48. fit <- glmLRT(fit) 
    50. DEG2=topTags(fit, n=Inf)
    51. class(DEG2)
    53. DEG2=as.data.frame(DEG2)
    54. head(DEG2)
    55. k1 = (DEG2$PValue < pvalue_t)&(DEG2$logFC < -logFC_t);table(k1)
    57. k2 = (DEG2$PValue < pvalue_t)&(DEG2$logFC > logFC_t);table(k2)
    58. DEG2$change = ifelse(k1,"DOWN",ifelse(k2,"UP","NOT"))
    60. head(DEG2)
    61. table(DEG2$change)
    63. ###limma----
    64. library(limma)
    65. dge <- DGEList(counts=exp)
    66. dge <- calcNormFactors(dge)
    67. v <- voom(dge,design, normalize="quantile")
    69. design <- model.matrix(~Group)
    70. fit <- lmFit(v, design)
    71. fit= eBayes(fit)
    73. DEG3 = topTable(fit, coef=2, n=Inf)
    74. DEG3 = na.omit(DEG3)
    76. k1 = (DEG3$P.Value < pvalue_t)&(DEG3$logFC < -logFC_t);table(k1)
    78. k2 = (DEG3$P.Value < pvalue_t)&(DEG3$logFC > logFC_t);table(k2)
    80. DEG3$change = ifelse(k1,"DOWN",ifelse(k2,"UP","NOT"))
    81. table(DEG3$change)
    83. head(DEG3)
    85. tj = data.frame(deseq2 = as.integer(table(DEG1$change)),
    86.                 edgeR = as.integer(table(DEG2$change)),
    87.                 limma_voom = as.integer(table(DEG3$change)),
    88.                 row.names = c("down","not","up")
    89. );tj
    91. save(DEG1,DEG2,DEG3,Group,tj,file = paste0(proj,"_DEG.Rdata"))

    PCA图

    1. rm(list = ls()) 
    2. library(ggplot2)
    3. library(tinyarray)
    4. load("TCGA-HNSC.Rdata")
    5. load("TCGA-HNSC_DEG.Rdata")
    6. exp[1:4,1:4]
    7. ncol(exp)
    8. # cpm 去除文库大小的影响
    9. dat = log2(cpm(exp)+1)
    10. pca.plot = draw_pca(dat,Group);pca.plot
    11. save(pca.plot,file = paste0(proj,"_pcaplot.Rdata"))

    TCGA作图，拼图

    1. install.packages('VennDiagram')
    2. library('VennDiagram')
    3. UP=function(df){
    4.   rownames(df)[df$change=="UP"]
    5. }
    6. DOWN=function(df){
    7.   rownames(df)[df$change=="DOWN"]
    8. }
    10. up = intersect(intersect(UP(DEG1),UP(DEG2)),UP(DEG3))
    11. down = intersect(intersect(DOWN(DEG1),DOWN(DEG2)),DOWN(DEG3))
    12. dat = log2(cpm(exp)+1)
    13. hp = draw_heatmap(dat[c(up,down),],Group,n_cutoff = 2)
    15. #上调、下调基因分别画维恩图
    16. up_genes = list(Deseq2 = UP(DEG1),
    17.                 edgeR = UP(DEG2),
    18.                 limma = UP(DEG3))
    20. down_genes = list(Deseq2 = DOWN(DEG1),
    21.                   edgeR = DOWN(DEG2),
    22.                   limma = DOWN(DEG3))
    24. up.plot <- draw_venn(up_genes,"UPgene")
    25. down.plot <- draw_venn(down_genes,"DOWNgene")
    27. #维恩图拼图，终于搞定
    29. library(patchwork)
    30. #up.plot + down.plot
    31. # 拼图
    32. pca.plot + hp+up.plot +down.plot+ plot_layout(guides = "collect")+
    33.   plot_annotation(title = "HNSC",
    34.     tag_levels = "A")
    35. ggsave(paste0(proj,"_heat_ve_pca.png"),width = 15,height = 10)

    image.png