在分析TCGA数据库里的RNA-seq数据之前，先了解TCGA样本的id名称。
TCGA样本id中第14-15位代表分组信息，01-09是tumor(肿瘤)，10-29是normal。
1.准备R包
#如果安装不成功 可以换用BiocManager::install()方法继续安装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))install.packages("DESeq2")if(!require(edgeR))install.packages("edgeR")if(!require(limma))install.packages("limma")
2.载入数据
##rm(list = ls())清空当前空间中所有对象rm(list = ls())a = read.csv("TCGAbiolinks_HNSC_counts.csv")dim(a)#查看表达矩阵 确保第一列不是基因名 也就是如果第一列是基因名 将其设置为行名 然后删除第一列rownames(a) <- a[,1]a = a[,-1]
3.将样品分类
补充知识点1：字符串操作
字符串的基本操作类型：
查找和替换大小写转换字符数统计字符串链接和拆分
stringr使用指南
stringr函数主要分为四类：
字符操作：操作字符向量中的单个字符str_length，str_sub，str_dup；添加、移除和操作空白符：str_pad，str_trim，str_wrap；大小写转换处理：str_to_lower，str_to_upper，str_to_title；模式匹配函数：str_detect, str_subset, str_count, str_locate, str_locate_all, str_match, str_match_all, str_replace, str_replace_all, str_split_fix, str_split, str_extract, str_extract_all 。
单个字符的处理
#1.单个字符长度 相当于ncharstr_length("abc")#2.根据位置信息提取或替换字符 类似于substr()#substr(x,start,end) x代表要操作的字符串 start和end表明起始和结束位点 如果两者一致代表只取这一个位置的> x <- c("abcdef","ghijk")#第三个> str_sub(x,3,3)[1] "c" "i"#第二个到倒数第二个> str_sub(x,2, -2)[1] "bcde" "hij" #第三个替换为7> str_sub(x,3,3) <- 7> x[1] "ab7def" "gh7jk" #第三个到第四个> str_sub(x,3,4)[1] "7d" "7j"#3.重复字符串 不同于rep#将x中第一个字符串重复两次 第二个字符串重复三次str_dup(x,c(2,3))
剩余操作参考：R语言的字符操作
该部分的代码：
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)
4.差异基因分析
4.1理论基础：线性模型 设计矩阵和比较矩阵
参考文章：转录组入门（7）：差异表达分析，仍然需要去学习一下统计学知识和数据挖掘第六章：线性模型和广义线性模型。
统计课上会介绍如何使用t检验比较两个样本之间的差异，然后在样本比较多的时候适用方差分析确定样本间是否有差异。样本来自于正态分布的群体，或者随机独立的大量抽样。
对于基因芯片的差异分析而言，由于普遍认为其数据是服从正态分布的，因此差异表达分析无非就是用t检验或方差分析应用到每一个基因上。高通量一次找的基因多，于是需要对多重试验进行矫正，控制假阳性，目前用的最多的就是limma。
但是，高通量测序(HTS)的read count普遍认为是服从泊松分布，不可能直接用正态分布的t检验和方差分析。 当然我们可以简单粗暴的使用对于的非参数检验的方法，但是统计力不够，结果的p值矫正之估计一个差异基因都找不到。因此，还是得要用参数检验的方法，于是就要说到方差分析和线性模型之间的关系了。
线性回归 一般是用于量化的预测变量来预测量化的响应变量。
负二向分布有两个参数，均值(mean)和离散值(dispersion)，离散值描述方差偏离均值的程度。泊松分布可以认为是负二向分布的离散值为1，也就是均值等于方差(mean=variance)的情况。
4.2 DESeq方法做差异分析
知识补充见：RNAseq分析全过程的DEseq2筛选差异表达基因。
分析后的几个参数解释
log2FoldChange：样本质检表达量的差异倍数，log2FoldChange的意思是取log2，这样可以让差异特别大的和差异特别小的数值之间缩小之间的差距。
Qvalue：是p-value的校正值，P值是统计差异的显著性的，Q值比P值更严格的一种统计。P-value是t-test来的。
padj：矫正过后的P值。
分析结果
#查看结果 contrast中是两两比较 一般是处理/对照> res <- results(dds,contrast = c("condition",rev(levels(group_list))))> #按pvalue值排序> resordered <- res[order(res$pvalue),]#转化为数据框> DEG <- as.data.frame((resordered))> head(DEG)baseMean log2FoldChangelfcSEstat pvalue padjENSG00000236780.7 25.77041-4.695302 0.2211499 -21.23130 4.908132e-100 1.957412e-95ENSG00000106351.13 1263.63358-2.366012 0.1199836 -19.71946 1.467863e-86 2.926993e-82ENSG00000107159.13 1796.89319 5.853589 0.2971668 19.69799 2.243523e-86 2.982464e-82ENSG00000244040.7 15.65783-4.132633 0.2112587 -19.56196 3.262933e-85 3.253226e-81ENSG00000070081.17 2879.11650-2.142581 0.1104093 -19.40581 6.892827e-84 5.497857e-80ENSG00000102547.19 431.55310-2.211553 0.1149789 -19.23443 1.906221e-82 1.267034e-78> head(resordered)log2 fold change (MLE): condition tumor vs normal Wald test p-value: condition tumor vs normal DataFrame with 6 rows and 6 columnsbaseMean log2FoldChangelfcSE<numeric><numeric> <numeric>ENSG00000236780.7 25.7704 -4.69530 0.221150ENSG00000106351.13 1263.6336 -2.36601 0.119984ENSG00000107159.13 1796.8932 5.85359 0.297167ENSG00000244040.7 15.6578 -4.13263 0.211259ENSG00000070081.17 2879.1165 -2.14258 0.110409ENSG00000102547.19 431.5531 -2.21155 0.114979stat pvalue padj<numeric> <numeric> <numeric>ENSG00000236780.7 -21.2313 4.90813e-100 1.95741e-95ENSG00000106351.13 -19.7195 1.46786e-86 2.92699e-82ENSG00000107159.13 19.6980 2.24352e-86 2.98246e-82ENSG00000244040.7 -19.5620 3.26293e-85 3.25323e-81ENSG00000070081.17 -19.4058 6.89283e-84 5.49786e-80ENSG00000102547.19 -19.2344 1.90622e-82 1.26703e-78#统计结果 有14382个基因上调 9214个基因下调> summary(resordered)out of 57747 with nonzero total read countadjusted p-value < 0.1LFC > 0 (up) : 14382, 25%LFC < 0 (down): 9214, 16%outliers [1] : 0, 0%low counts [2]: 17867, 31%(mean count < 0)[1] see 'cooksCutoff' argument of ?results[2] see 'independentFiltering' argument of #统计p值小于0.05的个数> table(res$padj<0.05)FALSE TRUE 18665 21216 > table(resordered$padj<0.05)FALSE TRUE 18665 21216
提取结果中的差异表达基因
#对DESeq2分析后具有显著性差异的结果进行提取 保存#获取padj(p值经过多重校验校正后的值)小于0.05 表达倍数取以2为对数后大于1或小于-1的差异表达基因#subset()函数过滤需要的结果到新变量diff_gene_Group2中#垚垚爸爱学习的方法> diff_gene_Group <- subset(resordered,padj < 0.05 & abs(log2FoldChange)>1)> dim(diff_gene_Group)[1] 106436> head(diff_gene_Group)log2 fold change (MLE): condition tumor vs normal Wald test p-value: condition tumor vs normal DataFrame with 6 rows and 6 columnsbaseMean log2FoldChangelfcSEstat pvalue padj<numeric><numeric> <numeric> <numeric> <numeric> <numeric>ENSG00000236780.7 25.7704 -4.69530 0.221150 -21.2313 4.90813e-100 1.95741e-95ENSG00000106351.13 1263.6336 -2.36601 0.119984 -19.7195 1.46786e-86 2.92699e-82ENSG00000107159.13 1796.8932 5.85359 0.297167 19.6980 2.24352e-86 2.98246e-82ENSG00000244040.7 15.6578 -4.13263 0.211259 -19.5620 3.26293e-85 3.25323e-81ENSG00000070081.17 2879.1165 -2.14258 0.110409 -19.4058 6.89283e-84 5.49786e-80ENSG00000102547.19 431.5531 -2.21155 0.114979 -19.2344 1.90622e-82 1.26703e-78> write.csv(diff_gene_Group,file = "DESeq2_diff_gene_Group")#生信start_site的过滤方法> DEG <-na.omit((DEG))> logFC_cutoff <- with(DEG,mean(abs(log2FoldChange)) + 2*sd(abs(log2FoldChange)))> 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 log2FoldChangelfcSEstat pvalue padj changeENSG00000236780.7 25.77041-4.695302 0.2211499 -21.23130 4.908132e-100 1.957412e-95 DOWNENSG00000106351.13 1263.63358-2.366012 0.1199836 -19.71946 1.467863e-86 2.926993e-82 NOTENSG00000107159.13 1796.89319 5.853589 0.2971668 19.69799 2.243523e-86 2.982464e-82UPENSG00000244040.7 15.65783-4.132633 0.2112587 -19.56196 3.262933e-85 3.253226e-81 DOWNENSG00000070081.17 2879.11650-2.142581 0.1104093 -19.40581 6.892827e-84 5.497857e-80 NOTENSG00000102547.19 431.55310-2.211553 0.1149789 -19.23443 1.906221e-82 1.267034e-78 NOT> table(DEG$change)DOWN NOT UP 920 37905 1056
4.3 edgeR方法做差异分析
利用DGEList函数读取count matrix数据，构建DGEList对象，构建该对象需要提供两类信息：表达量矩阵和分组信息。
1.构建DGEList对象
library(edgeR)dge <- DGEList(counts=a,group=group_list)
2.过滤counts数据
与DESeq2的预过滤不同，DESeq2的预过滤只是为了改善后续运算性能，在运行过程中依旧会自动处理low count数据，edgeR需要在分析前就要排除那些low count数据，而且非常严格。从生物学角度，有生物学意义的基因的表达量必须高于某一个阈值。从统计学角度上， low count的数据不太可能有显著性差异，而且在多重试验矫正阶段还会拖后腿。
dge$samples$lib.size <- colSums(dge$counts)
3.根据组成偏好(composition bias)标准化
edgeR的calcNormFactors函数使用TMM算法对DGEList标准化。
dge <- calcNormFactors(dge)
4.实验设计矩阵(Design matrix)
类似于DESeq2中的design参数，edgeR的线性模型和差异表达分析需要定义一个实验设计矩阵。
design <- model.matrix(~0+group_list)rownames(design) <- colnames(dge)colnames(design) <- levels(group_list)
5.估计离散值
负二项分布需要均值和离散值两个参数，edgeR对每个基因都估测一个经验贝叶斯稳健离散值(mpirical Bayes moderated dispersion)，还有一个公共离散值(common dispersion，所有基因的经验贝叶斯稳健离散值的均值)以及一个趋势离散值。
估计离散值这个步骤其实有许多estimate*Disp函数。当不存在实验设计矩阵(design matrix)的时候，estimateDisp 等价于 estimateCommonDisp 和 estimateTagwiseDisp 。而当给定实验设计矩阵(design matrix)时， estimateDisp 等价于 estimateGLMCommonDisp, estimateGLMTrendedDisp 和 estimateGLMTagwiseDisp。 其中tag与gene同义。
dge <- estimateGLMCommonDisp(dge,design)dge <- estimateGLMTrendedDisp(dge,design)dge <- estimateGLMTagwiseDisp(dge,design)#进一步通过quasi-likelihood(QL)拟合NB模型 用于解释生物学和技术性导致的基因特异性变异fit <- glmFit(dge,design)
6.差异表达检验
主要构建比较矩阵，类似于DESeq2中的results函数的contrast参数。
这里用的是glmQLFTest而不是glmLRT是因为前面用了glmQLTFit进行拟合，所以需要用QL F-test进行检验。如果前面用的是glmFit，那么对应的就是glmLRT. 作者称QL F-test更加严格。多重试验矫正用的也是BH方法。
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)))#生信start_site方法> DEG2$change = as.factor(+ ifelse(DEG2$PValue < 0.05 & abs(DEG2$logFC) > logFC_cutoff2,+ifelse(DEG2$logFC > logFC_cutoff2,'UP','DOWN'),'NOT'))> head(DEG2)logFClogCPM LR PValue FDR changeENSG00000170369.4 -10.173381 4.57836960 1304.702 1.075140e-285 6.521799e-281 DOWNENSG00000162877.13 -6.593453 -0.06574913 1293.927 2.360343e-283 7.158920e-279 DOWNENSG00000131686.15 -10.027389 4.34086408 1199.671 7.190846e-263 1.453989e-258 DOWNENSG00000101441.4 -11.033503 5.55328702 1105.591 2.011992e-242 3.051185e-238 DOWNENSG00000160349.10 -11.086792 6.52774641 1067.857 3.198968e-234 3.880989e-230 DOWNENSG00000167748.11 -5.976151 4.62530984 1039.784 4.044347e-228 4.088834e-224 DOWN> table(DEG2$change)DOWN NOT UP 1281 57855 1524 > edgeR_diff_gene_Group <- subset(DEG2,PValue < 0.05 & abs(logFC)>1)> dim(edgeR_diff_gene_Group)[1] 130766
4.4 limma-voom方法做差异分析
#生成designlibrary(limma)design <- model.matrix(~0+group_list)rownames(design) <- colnames(a)colnames(design) <- levels(group_list)#归一化处理dge <- DGEList(counts=a)dge <- calcNormFactors(dge)#limma-trendlogCPM <- cpm(dge,log = TRUE,prior.count = 3)fit <- lmFit(logCPM, design)fit <- eBayes(fit, trend=TRUE)topTable(fit, coef=constrasts)#limma-voomv <- voom(dge,design,normalize = "quantile")#ImFit线性模型拟合fit <- lmFit(v,design)#levels(group_list)的结果是normal tumor 而差异分析里需要tumor normal即处理/对照 则用rev反转一下#constrasts = tumor-normalconstrasts = paste(rev(levels(group_list)),collapse = "-")#比较矩阵(tumor/normal)cont.matrix <-makeContrasts(contrasts = constrasts,levels=design)#根据比对模型进行差值计算fit3 = contrasts.fit(fit,cont.matrix)#eBayes贝叶斯检验fit3 = eBayes((fit3))#生成所有基因的检测报告DEG3 = topTable(fit3,coef=constrasts,n=Inf)DEG3 = na.omit(DEG3)#计算logFC阈值logFC_cutoff3 <- with(DEG2,mean(abs(logFC))+2*sd(abs(logFC)))#筛选符合阈值要求的基因 并添加一列change 表明是基因是上调还是下调DEG3$change = as.factor(ifelse(DEG3$P.Value < 0.05 & abs(DEG2$logFC) > logFC_cutoff3,ifelse(DEG3$logFC > logFC_cutoff3,'UP','DOWN'),'NOT'))> head(DEG3)logFC AveExpr t P.Valueadj.P.Val B changeENSG00000260976.2 3.03 -3.3102384 27.37944 1.114750e-105 6.762072e-101 230.4398UPENSG00000096006.12 -8.519398 -1.0377074 -27.11420 2.541990e-104 7.709856e-100 227.1797 DOWNENSG00000170178.7 3.032367 -3.6654207 26.55944 1.782659e-101 3.604537e-97 220.8719UPENSG00000203740.4 3.824983 -2.6184232 26.14450 2.422231e-99 3.673314e-95 215.7615UPENSG00000181355.21 4.02 -2.4229131 26.02924 9.491357e-99 1.151491e-94 214.2791UPENSG00000167588.13 -4.750523 -0.471 -24.49507 7.788336e-91 7.874008e-87 195.9838 DOWN#统计上下调基因个数> table(DEG3$change)DOWN NOT UP 2639 57864 157
4.5 保存三种方法得到的矩阵
DESeq2_DEG <- DEGedgeR_DEG <- DEG2limma_voom_DEG <- DEG3save(DESeq2_DEG,edgeR_DEG,limma_voom_DEG,group_list,file = "DEG.Rdata")
4.6 热图
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)#DESeq2_DEG矩阵的热图#mat是每一个基因ID对应的第一个样品的表达量值mat = a[cg1,1]n=t(scale(t(mat)))n[n>1]=1n[n<-1]=-1ac=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)
