在分析TCGA數(shù)據(jù)庫里的RNA-seq數(shù)據(jù)之前涯穷，先了解一下TCGA樣本的id名稱里的小秘密：參考文章：TCGA的樣本id里藏著分組信息

根據(jù)文章里的內(nèi)容趋观，我查看了前一篇文章里我下載的count文件（利用R包TCGAbiolinks進(jìn)行各種數(shù)據(jù)下載）粱侣，打開是這樣的：

列是樣品名稱，格式舉例：TCGA.BA.A4IF.01A.11R.A266.07平痰，這里面包含分組信息艺骂，比如說這個(gè)樣品是腫瘤樣品還是正常組織诸老。分組信息是在這個(gè)id的第14-15位，01-09是tumor钳恕，10-29是normal别伏。比如第一個(gè)樣品，是01A忧额，說明這個(gè)樣品是個(gè)腫瘤樣品厘肮。

OK,知道了哪里是分組信息，就可以開始進(jìn)行操作了宙址。

（一）準(zhǔn)備R包

if(!require(stringr))install.packages('stringr')
if(!require(ggplotify))install.packages("ggplotify")
if(!require(patchwork))install.packages("patchwork")
if(!require(cowplot))install.packages("cowplot")
if(!require(DESeq2))BiocManager::install('DESeq2')
if(!require(edgeR))BiocManager::install('edgeR')
if(!require(limma))BiocManager::install('limma')

（二）載入數(shù)據(jù)

這里的數(shù)據(jù)是我上一篇文章下載好的count文件轴脐，不知道怎么下的請(qǐng)移步：利用R包TCGAbiolinks進(jìn)行各種數(shù)據(jù)下載

#加載數(shù)據(jù)
> rm(list = ls())
> a= read.csv("TCGAbiolinks_HNSC_counts.csv")
> dim(a)
[1] 56512   547
#另外要關(guān)注一下，你的這個(gè)表達(dá)矩陣的第一列是樣品還是基因名抡砂，如果第一列是基因名，要把第一列設(shè)置成為行名恬涧，不然只有的差異分析會(huì)出錯(cuò)
#將第一列換成行名
> row.names(a) <- a[, 1]
> a <- a[, -1]

NOTE：如果你沒有check你的矩陣注益，導(dǎo)致了你的矩陣第一列是基因名，后面DESeq2運(yùn)行會(huì)顯示：“Error in DESeqDataSet(se, design = design, ignoreRank) : some values in assay are negative”這樣的報(bào)錯(cuò)溯捆。

（三）將樣品分組

> group_list <- ifelse(as.numeric(str_sub(colnames(a),14,15))<10,"tumor","normal")
> group_list <- factor(group_list,levels = c("normal","tumor"))
> table(group_list)
group_list
normal  tumor 
    44    502 

（四）差異分析

參考文章：TCGA-6.（轉(zhuǎn)錄組）差異分析三大R包及其結(jié)果對(duì)比

DESeq2方法做差異分析

> library(DESeq2)
> colData <- data.frame(row.names =colnames(a), 
                       condition=group_list) #condition是你的分組信息
> dds <- DESeqDataSetFromMatrix(
   countData = a,
   colData = colData,
   design = ~ condition)
> dds <- DESeq(dds)
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
-- replacing outliers and refitting for 3594 genes
-- DESeq argument 'minReplicatesForReplace' = 7 
-- original counts are preserved in counts(dds)
estimating dispersions
fitting model and testing

# 兩兩比較
> res <- results(dds, contrast = c("condition",rev(levels(group_list))))
> resOrdered <- res[order(res$pvalue),] # 按照P值排序
> DEG <- as.data.frame(resOrdered)
> head(DEG)
                 baseMean log2FoldChange     lfcSE      stat        pvalue          padj
ENSG00000231887  233.5976      -7.414328 0.3208937 -23.10525 4.100424e-118 1.594778e-113
ENSG00000107159 1962.6861       5.849713 0.2969441  19.69971  2.168630e-86  4.217225e-82
ENSG00000106351 1273.1916      -2.371015 0.1206022 -19.65980  4.765762e-86  6.178493e-82
ENSG00000070081 2922.2020      -2.145545 0.1094426 -19.60430  1.420977e-85  1.381652e-81
ENSG00000102547  434.7033      -2.210238 0.1154164 -19.15013  9.655057e-82  7.510283e-78
ENSG00000151882  882.7880      -4.178352 0.2216513 -18.85102  2.882343e-79  1.868383e-75

去除NA值:

> DEG <- na.omit(DEG)

在DEG里添加一列名為change列丑搔，標(biāo)記基因上調(diào)下調(diào):

> logFC_cutoff <- with(DEG,mean(abs(log2FoldChange)) + 2*sd(abs(log2FoldChange)) )
> logFC_cutoff
[1] 2.610411
> DEG$change = as.factor(
+   ifelse(DEG$pvalue < 0.05 & abs(DEG$log2FoldChange) > logFC_cutoff,
+          ifelse(DEG$log2FoldChange > logFC_cutoff ,'UP','DOWN'),'NOT')
+ )
> head(DEG)
                 baseMean log2FoldChange     lfcSE      stat        pvalue          padj change
ENSG00000231887  233.5976      -7.414328 0.3208937 -23.10525 4.100424e-118 1.594778e-113   DOWN
ENSG00000107159 1962.6861       5.849713 0.2969441  19.69971  2.168630e-86  4.217225e-82     UP
ENSG00000106351 1273.1916      -2.371015 0.1206022 -19.65980  4.765762e-86  6.178493e-82    NOT
ENSG00000070081 2922.2020      -2.145545 0.1094426 -19.60430  1.420977e-85  1.381652e-81    NOT
ENSG00000102547  434.7033      -2.210238 0.1154164 -19.15013  9.655057e-82  7.510283e-78    NOT
ENSG00000151882  882.7880      -4.178352 0.2216513 -18.85102  2.882343e-79  1.868383e-75   DOWN
> View(DEG)
> DESeq2_DEG <- DEG

最后看一下用這個(gè)方法做差異分析，上調(diào)和下調(diào)的基因有多少：

> table(DEG$change)
 DOWN   NOT    UP 
  929 36955  1009

edgeR方法做差異分析

> library(edgeR)
> dge <- DGEList(counts=a,group=group_list)
> dge$samples$lib.size <- colSums(dge$counts)
> dge <- calcNormFactors(dge)
> design <- model.matrix(~0+group_list)
> rownames(design)<-colnames(dge)
> colnames(design)<-levels(group_list)
> dge <- estimateGLMCommonDisp(dge,design)
> dge <- estimateGLMTrendedDisp(dge, design)
> dge <- estimateGLMTagwiseDisp(dge, design)
> fit <- glmFit(dge, design)
> fit2 <- glmLRT(fit, contrast=c(-1,1)) 
> DEG2=topTags(fit2, n=nrow(a))
> DEG2=as.data.frame(DEG2)
> logFC_cutoff2 <- with(DEG2,mean(abs(logFC)) + 2*sd(abs(logFC)) )
> DEG2$change = as.factor(
+   ifelse(DEG2$PValue < 0.05 & abs(DEG2$logFC) > logFC_cutoff2,
+          ifelse(DEG2$logFC > logFC_cutoff2 ,'UP','DOWN'),'NOT')
+ )
> head(DEG2)
                     logFC     logCPM       LR        PValue           FDR change
ENSG00000170369  -9.842798  4.4837048 1227.908 5.248913e-269 2.966266e-264   DOWN
ENSG00000162877  -6.470581 -0.1480235 1215.091 3.203213e-266 9.050998e-262   DOWN
ENSG00000231887  -9.129113  4.9924958 1140.670 4.781111e-250 9.006338e-246   DOWN
ENSG00000131686  -9.939181  4.2888099 1135.359 6.822362e-249 9.638633e-245   DOWN
ENSG00000101441 -10.903794  5.4644612 1059.704 1.893018e-232 2.139565e-228   DOWN
ENSG00000160349 -10.982631  6.4426001 1030.047 5.287725e-226 4.980332e-222   DOWN

看一下用edgeR方法提揍，上調(diào)和下調(diào)基因的數(shù)量：

> table(DEG2$change)
 DOWN   NOT    UP 
 1188 53914  1410 
> edgeR_DEG <- DEG2

limma-voom方法做差異分析

> library(limma)
> design <- model.matrix(~0+group_list)
> colnames(design)=levels(group_list)
> rownames(design)=colnames(a)
> dge <- DGEList(counts=a)
> dge <- calcNormFactors(dge)
> logCPM <- cpm(dge, log=TRUE, prior.count=3)
> v <- voom(dge,design, normalize="quantile")
> fit <- lmFit(v, design)
> constrasts = paste(rev(levels(group_list)),collapse = "-")
> cont.matrix <- makeContrasts(contrasts=constrasts,levels = design) 
> fit3=contrasts.fit(fit,cont.matrix)
> fit3=eBayes(fit3)
> DEG3 = topTable(fit3, coef=constrasts, n=Inf)
> DEG3 = na.omit(DEG3)
> logFC_cutoff3 <- with(DEG3,mean(abs(logFC)) + 2*sd(abs(logFC)) )
> DEG3$change = as.factor(
+   ifelse(DEG3$P.Value < 0.05 & abs(DEG3$logFC) > logFC_cutoff3,
+          ifelse(DEG3$logFC > logFC_cutoff3 ,'UP','DOWN'),'NOT')
+ )
> head(DEG3)
                    logFC   AveExpr         t       P.Value     adj.P.Val        B change
ENSG00000203740  3.394627 -3.404611  28.73533 3.090631e-111 1.746578e-106 243.1880     UP
ENSG00000096006 -8.588036 -1.083975 -27.14114 2.835110e-103  8.010887e-99 224.8165   DOWN
ENSG00000260976  3.045123 -3.341040  27.00751 1.329154e-102  2.503772e-98 223.4002     UP
ENSG00000181092 -7.489810 -4.776791 -24.73134  4.130357e-91  5.835368e-87 196.2190   DOWN
ENSG00000198478 -4.043231  3.231204 -24.35482  3.353473e-89  3.790229e-85 192.7091   DOWN
ENSG00000181355  4.163352 -2.440333  24.29635  6.639596e-89  6.253614e-85 191.7566     UP
> limma_voom_DEG <- DEG3

> table(DEG3$change)

 DOWN   NOT    UP 
 1593 53742  1177 

（五）保存三種方法得到的矩陣

> save(DESeq2_DEG,edgeR_DEG,limma_voom_DEG,group_list,file = "DEG.Rdata")

（六）熱圖

下面得到了三個(gè)方法得來了矩陣啤月，就可以可視化了。教程里用的包draw_heatmap在我的R版本上安裝不了劳跃，所以我用的是pheatmap畫的熱圖谎仲。

#這一步是定義后面畫的熱圖是三個(gè)矩陣?yán)颿hange那一列上調(diào)和下調(diào)的基因，不包括沒變化的基因
> cg1 = rownames(DESeq2_DEG)[DESeq2_DEG$change !="NOT"]
> cg2 = rownames(edgeR_DEG)[edgeR_DEG$change !="NOT"]
> cg3 = rownames(limma_voom_DEG)[limma_voom_DEG$change !="NOT"]

> library(pheatmap)
> library(RColorBrewer)
#定義熱圖的顏色
> color<- colorRampPalette(c('#436eee','white','#EE0000'))(100) 
#下面畫第一個(gè)DESeq2矩陣的熱圖
> mat=a[cg1,]
> n=t(scale(t(mat)))
> n[n>1]=1
> n[n< -1]= -1
> ac=data.frame(group=group_list)
> rownames(ac)=colnames(mat)
> ht1 <- pheatmap(n,show_rownames = F,show_colnames = F, 
         cluster_rows = F,cluster_cols = T,
         annotation_col = ac,color=color)

#下面畫edgeR矩陣的熱圖
> mat2=a[cg2,]
> n2=t(scale(t(mat2)))
> n2[n2 > 1]=1
> n2[n2< -1]= -1
> ht2 <- pheatmap(n2,show_rownames = F,show_colnames = F, 
         cluster_rows = F,cluster_cols = T,
         annotation_col = ac,color=color)

#下面畫limma矩陣的熱圖
> mat3=a[cg3,]
> n3=t(scale(t(mat3)))
> n3[n3 > 1]=1
> n3[n3< -1]= -1
> ht3 <- pheatmap(n3,show_rownames = F,show_colnames = F, 
         cluster_rows = F,cluster_cols = T,
         annotation_col = ac,color=color)

DESeq2矩陣熱圖

edgeR矩陣熱圖

limma矩陣熱圖

（七）火山圖

> library(EnhancedVolcano)
> library(airway)

> v1=EnhancedVolcano(DESeq2_DEG,
                lab = rownames(DESeq2_DEG),
                x = 'log2FoldChange',
                y = 'pvalue',
                xlim = c(-8, 8),
                title = 'DESeq2_DEG',
                pCutoff = 10e-17,
                FCcutoff = 2.5,
                transcriptPointSize = 1.0,
                transcriptLabSize = 3.0,
                col=c('black', 'blue', 'green', 'red1'),
                colAlpha = 1,
                legend=c('NS','Log2 fold-change','P value',
                         'P value & Log2 fold-change'),
                legendPosition = 'right',
                legendLabSize = 10,
                legendIconSize = 3.0,
)

> v2=EnhancedVolcano(edgeR_DEG,
                lab = rownames(edgeR_DEG),
                x = 'logFC',
                y = 'PValue',
                xlim = c(-8, 8),
                title = 'edgeR_DEG',
                pCutoff = 10e-17,
                FCcutoff = 2.5,
                transcriptPointSize = 1.0,
                transcriptLabSize = 3.0,
                col=c('black', 'blue', 'green', 'red1'),
                colAlpha = 1,
                legend=c('NS','LogFC','P value',
                         'P value & LogFC'),
                legendPosition = 'right',
                legendLabSize = 10,
                legendIconSize = 3.0,
)

> v3=EnhancedVolcano(limma_voom_DEG,
                lab = rownames(limma_voom_DEG),
                x = 'logFC',
                y = 'P.Value',
                xlim = c(-8, 8),
                title = 'limma_voom_DEG',
                pCutoff = 10e-17,
                FCcutoff = 2.5,
                transcriptPointSize = 1.0,
                transcriptLabSize = 3.0,
                col=c('black', 'blue', 'green', 'red1'),
                colAlpha = 1,
                legend=c('NS','LogFC','P value',
                         'P value & LogFC'),
                legendPosition = 'right',
                legendLabSize = 10,
                legendIconSize = 3.0,
)

把三個(gè)火山圖組合起來：

#multiplot function
> multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  require(grid)
  
  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)
  
  numPlots = length(plots)
  
  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                     ncol = cols, nrow = ceiling(numPlots/cols))
  }
  
  if (numPlots==1) {
    print(plots[[1]])
    
  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
      
      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}

> multiplot(v1,v2,v3,cols=3)

（八）三個(gè)不同方法得到的矩陣的交集可視化

上面的可視化是用三個(gè)矩陣分別畫出的熱圖和火山圖刨仑，下面要把這三個(gè)矩陣?yán)锷险{(diào)和下調(diào)的基因的交集提出來郑诺，然后再進(jìn)行可視化：

> UP=function(df){
  rownames(df)[df$change=="UP"]
}
> DOWN=function(df){
  rownames(df)[df$change=="DOWN"]
}

> up = intersect(intersect(UP(DESeq2_DEG),UP(edgeR_DEG)),UP(limma_voom_DEG))
> down = intersect(intersect(DOWN(DESeq2_DEG),DOWN(edgeR_DEG)),DOWN(limma_voom_DEG))

> mat_total=a[c(up,down),]
> n4=t(scale(t(mat_total)))
> n4[n4 >1]=1
> n4[n4< -1]= -1
> ac=data.frame(group=group_list)
> rownames(ac)=colnames(mat_total)
> ht_combine <- pheatmap(n4,show_rownames = F,show_colnames = F, 
                cluster_rows = F,cluster_cols = T,
                annotation_col = ac,color=color)

三個(gè)矩陣的交集畫韋恩圖：

#上調(diào)、下調(diào)基因分別畫維恩圖
> up.plot <- venn(list(UP(DESeq2_DEG),UP(edgeR_DEG),UP(limma_voom_DEG)))
> down.plot <- venn(list(DOWN(DESeq2_DEG),DOWN(edgeR_DEG),DOWN(limma_voom_DEG)))

三個(gè)矩陣?yán)锷险{(diào)基因的交集杉武，一共533個(gè)

三個(gè)矩陣?yán)锵抡{(diào)基因的交集辙诞，一共630個(gè)