TCGA数据库学习二：差异基因分析

多学点生信怎么了

已于 2024-03-04 23:18:09 修改

阅读量4.9k

点赞数 4

分类专栏： TCGA数据库学习文章标签：数据库学习 r语言数据挖掘

于 2023-03-31 21:40:16 首次发布

本文链接：https://blog.csdn.net/weixin_73362123/article/details/129791602

版权

TCGA数据库学习专栏收录该内容

4 篇文章 8 订阅

订阅专栏

在分析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.单个字符长度 相当于nchar
str_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 log2FoldChange     lfcSE      stat        pvalue         padj
ENSG00000236780.7    25.77041      -4.695302 0.2211499 -21.23130 4.908132e-100 1.957412e-95
ENSG00000106351.13 1263.63358      -2.366012 0.1199836 -19.71946  1.467863e-86 2.926993e-82
ENSG00000107159.13 1796.89319       5.853589 0.2971668  19.69799  2.243523e-86 2.982464e-82
ENSG00000244040.7    15.65783      -4.132633 0.2112587 -19.56196  3.262933e-85 3.253226e-81
ENSG00000070081.17 2879.11650      -2.142581 0.1104093 -19.40581  6.892827e-84 5.497857e-80
ENSG00000102547.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 columns
                    baseMean log2FoldChange     lfcSE
                   <numeric>      <numeric> <numeric>
ENSG00000236780.7    25.7704       -4.69530  0.221150
ENSG00000106351.13 1263.6336       -2.36601  0.119984
ENSG00000107159.13 1796.8932        5.85359  0.297167
ENSG00000244040.7    15.6578       -4.13263  0.211259
ENSG00000070081.17 2879.1165       -2.14258  0.110409
ENSG00000102547.19  431.5531       -2.21155  0.114979
                        stat       pvalue        padj
                   <numeric>    <numeric>   <numeric>
ENSG00000236780.7   -21.2313 4.90813e-100 1.95741e-95
ENSG00000106351.13  -19.7195  1.46786e-86 2.92699e-82
ENSG00000107159.13   19.6980  2.24352e-86 2.98246e-82
ENSG00000244040.7   -19.5620  3.26293e-85 3.25323e-81
ENSG00000070081.17  -19.4058  6.89283e-84 5.49786e-80
ENSG00000102547.19  -19.2344  1.90622e-82 1.26703e-78

#统计结果 有14382个基因上调 9214个基因下调
> summary(resordered)
out of 57747 with nonzero total read count
adjusted p-value < 0.1
LFC > 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] 10643     6
> 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 columns
                    baseMean log2FoldChange     lfcSE      stat       pvalue        padj
                   <numeric>      <numeric> <numeric> <numeric>    <numeric>   <numeric>
ENSG00000236780.7    25.7704       -4.69530  0.221150  -21.2313 4.90813e-100 1.95741e-95
ENSG00000106351.13 1263.6336       -2.36601  0.119984  -19.7195  1.46786e-86 2.92699e-82
ENSG00000107159.13 1796.8932        5.85359  0.297167   19.6980  2.24352e-86 2.98246e-82
ENSG00000244040.7    15.6578       -4.13263  0.211259  -19.5620  3.26293e-85 3.25323e-81
ENSG00000070081.17 2879.1165       -2.14258  0.110409  -19.4058  6.89283e-84 5.49786e-80
ENSG00000102547.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 log2FoldChange     lfcSE      stat        pvalue         padj change
ENSG00000236780.7    25.77041      -4.695302 0.2211499 -21.23130 4.908132e-100 1.957412e-95   DOWN
ENSG00000106351.13 1263.63358      -2.366012 0.1199836 -19.71946  1.467863e-86 2.926993e-82    NOT
ENSG00000107159.13 1796.89319       5.853589 0.2971668  19.69799  2.243523e-86 2.982464e-82     UP
ENSG00000244040.7    15.65783      -4.132633 0.2112587 -19.56196  3.262933e-85 3.253226e-81   DOWN
ENSG00000070081.17 2879.11650      -2.142581 0.1104093 -19.40581  6.892827e-84 5.497857e-80    NOT
ENSG00000102547.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)
                        logFC      logCPM       LR        PValue           FDR change
ENSG00000170369.4  -10.173381  4.57836960 1304.702 1.075140e-285 6.521799e-281   DOWN
ENSG00000162877.13  -6.593453 -0.06574913 1293.927 2.360343e-283 7.158920e-279   DOWN
ENSG00000131686.15 -10.027389  4.34086408 1199.671 7.190846e-263 1.453989e-258   DOWN
ENSG00000101441.4  -11.033503  5.55328702 1105.591 2.011992e-242 3.051185e-238   DOWN
ENSG00000160349.10 -11.086792  6.52774641 1067.857 3.198968e-234 3.880989e-230   DOWN
ENSG00000167748.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] 13076     6

4.4 limma-voom方法做差异分析

#生成design
library(limma)
design <- model.matrix(~0+group_list)
rownames(design) <- colnames(a)
colnames(design) <- levels(group_list)

#归一化处理
dge <- DGEList(counts=a)
dge <- calcNormFactors(dge)

#limma-trend
logCPM <- cpm(dge,log = TRUE,prior.count = 3)

fit <- lmFit(logCPM, design)
fit <- eBayes(fit, trend=TRUE)
topTable(fit, coef=constrasts)

#limma-voom
v <- voom(dge,design,normalize = "quantile")

#ImFit线性模型拟合
fit <- lmFit(v,design)

#levels(group_list)的结果是normal tumor 而差异分析里需要tumor normal即处理/对照 则用rev反转一下
#constrasts = tumor-normal
constrasts = 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.Value     adj.P.Val        B change
ENSG00000260976.2   3.032008 -3.3102384  27.37944 1.114750e-105 6.762072e-101 230.4398     UP
ENSG00000096006.12 -8.519398 -1.0377074 -27.11420 2.541990e-104 7.709856e-100 227.1797   DOWN
ENSG00000170178.7   3.032367 -3.6654207  26.55944 1.782659e-101  3.604537e-97 220.8719     UP
ENSG00000203740.4   3.824983 -2.6184232  26.14450  2.422231e-99  3.673314e-95 215.7615     UP
ENSG00000181355.21  4.202102 -2.4229131  26.02924  9.491357e-99  1.151491e-94 214.2791     UP
ENSG00000167588.13 -4.750523 -0.4712017 -24.49507  7.788336e-91  7.874008e-87 195.9838   DOWN

#统计上下调基因个数
> table(DEG3$change)

 DOWN   NOT    UP 
 2639 57864   157

4.5 保存三种方法得到的矩阵

DESeq2_DEG <- DEG
edgeR_DEG <- DEG2
limma_voom_DEG <- DEG3
save(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]=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)