在分析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)