之前的一个推文是从UCSC XENA获取TCGA的表达和表型数据，然后利用代码对表达数据进行了ID注释，以及mRNA、lncRNA和miRNA的区分筛选，最后将患者ID和临床信息进行配比，用于后续分析。详细内容见推文：《利用R代码从UCSC XENA下载mRNA, lncRNA, miRNA表达数据并匹配临床信息》。
本期我将继续上次的内容，从TCGA 546个头颈癌数据集（Tumor = 502，Normal = 44）中，提取出43对癌和癌旁样本，并做配对差异分析， 然后绘制某个基因的配对箱式图。实际上TCGA好多癌症比如头颈癌、肝癌等，都有癌和癌旁的数据，且癌和癌旁都是一一配对的关系。所以在分析癌vs.癌旁的过程中，可以选择普通的差异分析，例如头颈癌的502个癌vs.44癌旁；另外一种思路是从中挑选出配对的癌vs.癌旁进行配对差异分析，例如头颈癌的43个癌vs.43个癌旁。实际上，同一个数据的非配对分析和配对分析差异还是很大的，详见我之前写过的一个帖子《差异分析|DESeq2完成配对样本的差异分析》。
关于DEseq2配对差异分析，很少有帖子涉及 (生信宝典注：不是很少有帖子涉及，而是看有没有发现，配对的个体信息可以视作一个批次因素，在DESeq2差异基因分析和批次效应移除, 高通量数据中批次效应的鉴定和处理 - 系列总结和更新, 典型医学设计实验GEO数据分析 (step-by-step) - Limma差异分析、火山图、功能富集中都有类似的处理方式)。
总结来说，拿到配对设计的数据，如果不进行配对分析而用常规的差异分析，这样的结果可能会大不相同。因此，在分析数据的时候，一定要明白实验设计。
下面代码展示了如何从546个样本中挑选出一一配对的癌和癌旁数据，并进行DEseq2配对差异分析。
remove(list = ls())
##加载包
suppressMessages(library(DESeq2))
suppressMessages(library(dplyr))
library(ggplot2)
library(ggpubr)
library(ggthemes)1. 加载数据
首先，读入上期推文中处理好的表达矩阵（仅含mRNA）和表型数据。
# 1.1 表达矩阵
expr_data <- read.csv("./Rawdata/TCGA_mRNA_count_log2.csv",header = T,row.names = 1,check.names=F)
expr_data[1:3,1:4]
## TCGA.BB.4224.01 TCGA.H7.7774.01 TCGA.CV.6943.01 TCGA.CN.5374.01
## A4GNT 2.584963 0.00000 1.00000 2.321928
## AAAS 10.952741 10.43463 11.27671 10.750707
## AACS 11.098032 13.62822 12.52943 11.171802
# 1.2 表型数据
phenotype <- read.csv("./Rawdata/TCGA_HNSC_phenotype.csv",header = T,row.names = 1,check.names=F, stringsAsFactors = T)
phenotype$subject <- as.factor(substring(rownames(phenotype),1,12))
phenotype[1:3,c("group","subject")]
## group subject
## TCGA.BB.4224.01 Tumor TCGA.BB.4224
## TCGA.H7.7774.01 Tumor TCGA.H7.7774
## TCGA.CV.6943.01 Tumor TCGA.CV.6943
# 检查一下表达数据和表型数据是否是一一对应
table(colnames(expr_data) == rownames(phenotype)) #T
##
## TRUE
## 546表达数据和表型数据如果不是一一对应的话，运行下面这行代码:
# 利用match函数进行匹配
phenotype <- phenotype[match(colnames(expr_data),rownames(phenotype)),]
identical(colnames(expr_data),rownames(phenotype))
## [1] TRUE2. 筛选出配对的样本
2.1 将546个样本分为两部分，Tumor组和Normal组
# Tumor组
Tumor_group <- phenotype[which(phenotype$group == "Tumor"),] #TCGA的Tumor样本,502个
table(Tumor_group$group) #502个癌
##
## Nontumor Tumor
## 0 502
# Normal组
Normal_group <- phenotype[which(phenotype$group == "Nontumor"),] #TCGA的control样本,44
table(Normal_group$group) #44个癌旁
##
## Nontumor Tumor
## 44 0
# 找到癌和癌旁能够进行一一配对的样本的信息
Normal_group <- Normal_group[Normal_group$subject %in% Tumor_group$subject,]
#仅保留配对的对照组
table(Normal_group$group) # 43个存在配对的癌旁
##
## Nontumor Tumor
## 43 0上面这行代码就获取了43个配对的癌旁样本信息，根据这个信息可以获得所有配对的癌+癌旁，代码如下：
# 获取表达矩阵
paired_expr_data <- expr_data[,substring(colnames(expr_data),1,12) %in%
substring(rownames(Normal_group),1,12)]
# 实际上，这里的substring(rownames(Normal_group),1,12)和Normal_group$ID是一样的
paired_phe <- phenotype[match(colnames(paired_expr_data),row.names(phenotype)),]
paired_phe <- droplevels(paired_phe)
dim(paired_expr_data) # 获得了86个样本表达信息，其中43个癌，43个癌旁，是一一配对的关系
## [1] 18192 86
write.csv(paired_expr_data,"./Output/TCGA_HNSC_mRNA_count_paired_43vs43.csv")然后，进行配对的DEseq2差异分析。这里只保留了至少在75%的样本中都有表达的基因。
3. DEseq2配对差异分析
3.1 配对的表达矩阵
# UCSC Xena下载的数据是经过log2（count+1）处理的
data = (2^paired_expr_data-1) %>% apply(2, as.integer)# 还原log2的处理
row.names(data) = row.names(paired_expr_data)
# 保留至少在75%的样本中都有表达的基因
keep_data <- rowSums(data> 0) >= floor(0.75*ncol(data))
table(keep_data)
## keep_data
## FALSE TRUE
## 2599 15593
data <- data[keep_data,]
data[1:4,1:4]
## TCGA.CV.6943.01 TCGA.CV.6959.01 TCGA.CV.7438.11 TCGA.CV.7242.11
## AAAS 2479 3304 1900 2190
## AACS 5911 3881 4759 3070
## AADAC 7 5 449 2618
## AADAT 140 77 160 127
dim(data)
## [1] 15593 863.2 对应的表型数据和配对信息
table(paired_phe$subject,paired_phe$group)
##
## Nontumor Tumor
## TCGA.CV.6933 1 1
## TCGA.CV.6934 1 1
## TCGA.CV.6935 1 1
## TCGA.CV.6936 1 1
## TCGA.CV.6938 1 1
## TCGA.CV.6939 1 1
## TCGA.CV.6943 1 1
## TCGA.CV.6955 1 1
## TCGA.CV.6956 1 1
## TCGA.CV.6959 1 1
## TCGA.CV.6960 1 1
## TCGA.CV.6961 1 1
## TCGA.CV.6962 1 1
## TCGA.CV.7091 1 1
## TCGA.CV.7097 1 1
## TCGA.CV.7101 1 1
## TCGA.CV.7103 1 1
## TCGA.CV.7177 1 1
## TCGA.CV.7178 1 1
## TCGA.CV.7183 1 1
## TCGA.CV.7235 1 1
## TCGA.CV.7238 1 1
## TCGA.CV.7242 1 1
## TCGA.CV.7245 1 1
## TCGA.CV.7250 1 1
## TCGA.CV.7252 1 1
## TCGA.CV.7255 1 1
## TCGA.CV.7261 1 1
## TCGA.CV.7406 1 1
## TCGA.CV.7416 1 1
## TCGA.CV.7423 1 1
## TCGA.CV.7424 1 1
## TCGA.CV.7425 1 1
## TCGA.CV.7432 1 1
## TCGA.CV.7434 1 1
## TCGA.CV.7437 1 1
## TCGA.CV.7438 1 1
## TCGA.CV.7440 1 1
## TCGA.H7.A6C4 1 1
## TCGA.HD.8635 1 1
## TCGA.HD.A6HZ 1 1
## TCGA.HD.A6I0 1 1
## TCGA.WA.A7GZ 1 13.3 配对差异分析
# DEseq2的配对差异分析跑的有点慢，43对我的电脑需要约10分钟
# 所以，这里我为了加快运行速度，选了5对5.
data = data [,order(paired_phe$subject)]
data = as.matrix(data[,1:10])
data[1:3,1:4]
## TCGA.CV.6933.01 TCGA.CV.6933.11 TCGA.CV.6934.01 TCGA.CV.6934.11
## AAAS 2993 1328 2816 1977
## AACS 5124 4997 3522 4108
## AADAC 98 1 16 63
coldata_paired <- paired_phe[colnames(data),]
coldata_paired[,c('subject','group')]
## subject group
## TCGA.CV.6933.01 TCGA.CV.6933 Tumor
## TCGA.CV.6933.11 TCGA.CV.6933 Nontumor
## TCGA.CV.6934.01 TCGA.CV.6934 Tumor
## TCGA.CV.6934.11 TCGA.CV.6934 Nontumor
## TCGA.CV.6935.01 TCGA.CV.6935 Tumor
## TCGA.CV.6935.11 TCGA.CV.6935 Nontumor
## TCGA.CV.6936.01 TCGA.CV.6936 Tumor
## TCGA.CV.6936.11 TCGA.CV.6936 Nontumor
## TCGA.CV.6938.01 TCGA.CV.6938 Tumor
## TCGA.CV.6938.11 TCGA.CV.6938 Nontumor
# 3.4.1 DEseq2的配对差异分析
dds_paired <- DESeqDataSetFromMatrix(countData = data,
colData = coldata_paired,
design = ~ subject + group)
## factor levels were dropped which had no samples
dds_paired$group <- relevel(dds_paired$group, ref = "Nontumor") # 指定哪一组作为对照
dds_paired <- DESeq(dds_paired)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
nrDEG_paired <- as.data.frame(results(dds_paired))
nrDEG_paired = nrDEG_paired[order(nrDEG_paired$log2FoldChange),] # 按照logFC排序
nrDEG_paired[1:3,1:6]
## baseMean log2FoldChange lfcSE stat pvalue padj
## PRH2 544288.00 -15.63870 2.509593 -6.231567 4.617911e-10 7.423411e-08
## STATH 1385580.80 -15.33676 1.292428 -11.866626 1.764414e-32 9.170836e-29
## MUC21 26826.61 -12.83806 1.239818 -10.354790 3.980606e-25 7.758699e-22
write.csv(nrDEG_paired,file = "./Output/TCGA_HNSCpaired_DESeq2_nrDEG.csv")
# 定义差异基因，用于后续的分析，例如火山图，GO和KEGG等
logFC_cutoff <- 2
nrDEG_paired$change <- as.factor(ifelse(nrDEG_paired$padj < 0.05 & abs(nrDEG_paired$log2FoldChange) > logFC_cutoff,
ifelse(nrDEG_paired$log2FoldChange > logFC_cutoff ,'UP','DOWN'),'NOT'))
table(nrDEG_paired$change)
##
## DOWN NOT UP
## 753 14452 388
# 绘制火山图
Volcano_data <- na.omit(nrDEG_paired)
Volcano_data$logP <- -log10(Volcano_data$padj)
Volcano_paired <- ggscatter(Volcano_data, x = "log2FoldChange", y = "logP",
color = "change",
palette = c("#2f5688", "#BBBBBB", "#CC0000"),
size = 1,
font.label = 10,
repel = T,
xlab = "log2 FoldChange",
ylab = "-log10 (pvalue)",
title="Volcano Plot") +
theme_base() +
geom_hline(yintercept = 1.30, linetype="dashed") +
geom_vline(xintercept = c(-logFC_cutoff,logFC_cutoff), linetype="dashed")+
theme(legend.position = "right",plot.title = element_text(size = 14,color="black",hjust = 0.5))
Volcano_paired
ggsave("./Output/Volcano_paired.PDF",Volcano_paired,width=6,height=5)上述就是配对DEseq2差异分析的一个结果。接下来，我们看一下如果使用普通DEseq2差异分析的结果，将会有哪些区别？
# 3.5 DEseq2的普通差异分析
dds_2 <- DESeqDataSetFromMatrix(countData = data,
colData = coldata_paired,
design = ~ group)
dds_2$condition <- relevel(dds_2$group, ref = "Nontumor") # 指定哪一组作为对照
dds_2 <- DESeq(dds_2)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
nrDEG_2 <- as.data.frame(results(dds_2))
nrDEG_2 = nrDEG_2[order(nrDEG_2$log2FoldChange),] # 按照logFC排序
nrDEG_2[1:3,1:4]
## baseMean log2FoldChange lfcSE stat
## PRH2 544288.0 -20.27705 2.722897 -7.446865
## STATH 1385580.8 -18.75144 2.069789 -9.059590
## PRH1 104842.6 -15.65276 1.820363 -8.598703
write.csv(nrDEG_2,file = "./Output/TCGA_HNSC_DESeq2_nrDEG.csv")
# 这里我还提取了标准化后的表达矩阵，可以用于后续的热图箱式图绘制等等，当然也可以使用FPKM或TPM代替
rld_paired <- vst(dds_paired)
Nor_expr_data_paired <- assay(rld_paired)
write.csv(Nor_expr_data_paired,file = "./Output/TCGA_HNSCpaired_Norexpr_data_paired.csv")
# 3.5 定义差异基因，用于后续的分析，例如GO和KEGG等
logFC_cutoff <- 2
nrDEG_2$change <- as.factor(ifelse(nrDEG_2$padj < 0.05 & abs(nrDEG_2$log2FoldChange) > logFC_cutoff,
ifelse(nrDEG_2$log2FoldChange > logFC_cutoff ,'UP','DOWN'),'NOT'))
table(nrDEG_2$change)
##
## DOWN NOT UP
## 626 14093 375
# 3.6 绘制火山图
Volcano_data_2 = na.omit(nrDEG_2)
Volcano_data_2$logP <- -log10(Volcano_data_2$padj)
Volcano_2 <- ggscatter(Volcano_data_2, x = "log2FoldChange", y = "logP",
color = "change",
palette = c("#2f5688", "#BBBBBB", "#CC0000"),
size = 1,
font.label = 10,
repel = T,
xlab = "log2 FoldChange",
ylab = "-log10 (pvalue)",
title="Volcano Plot") +
theme_base() +
geom_hline(yintercept = 1.30, linetype="dashed") +
geom_vline(xintercept = c(-logFC_cutoff,logFC_cutoff), linetype="dashed")+
theme(legend.position = "right",plot.title = element_text(size = 14,color="black",hjust = 0.5))
Volcano_2
ggsave("./Output/Volcano_2.PDF",Volcano_2 ,width=6,height=5)可以看到配对的DEseq2差异分析有753基因下调，388基因上调；而常规的差异分析有626个基因下调，375基因上调，这个结果看似差不多 (生信宝典注：这里应该做个Venn图比较下的，除了数量，具体基因是否一致。可参考高通量数据中批次效应的鉴定和处理 - 系列总结和更新)。
实际上，笔者曾处理过一个项目，用配对和不用配对，得到的结果相差甚远，主要还是体现在火山图、KEGG和GO分析等需要输入差异基因的分析中。但好像做GSEA这样纳入全部基因的分析，并不会受到太大的影响，这点有机会再探讨。
生信宝典注 (具体见一文掌握GSEA，超详细教程！)：
- 如果用log2FC做GSEA结果是肯定会有差别的，影响多大不好说。
- 如果用标准化后的矩阵做GSEA结果没差别，因为配对分析没有影响到表达标准化，而只是在差异分析步骤起作用。
绘制单基因配对t检验箱式图
4.1 提取PRH2基因的表达量
## 4. PRH2表达量的配对箱式图绘制
PRH2_expression <- t(paired_expr_data["PRH2",] ) %>% data.frame(PRH2 = .)
PRH2_expression$group <- paired_phe$group
PRH2_expression$paired <- paired_phe$subject
PRH2_expression = PRH2_expression[order(as.numeric(paired_phe$subject)),] #这里需要按照配对信息排序好，以做后需的配对差异分析4.2 开始绘图
# 这个是我自己写的一个ggplot2的主题，可以自定义修改其中的参数
if(T){mytheme <- theme(plot.title = element_text(size = 12,color="black",hjust = 0.5),
axis.title = element_text(size = 12,color ="black"),
axis.text = element_text(size= 12,color = "black"),
panel.grid.minor.y = element_blank(),
panel.grid.minor.x = element_blank(),
panel.grid=element_blank(),
legend.position = "none",
legend.text = element_text(size= 12),
legend.title= element_text(size= 12)
)}
# 配对t检验箱式图图
box_paired <- ggplot(PRH2_expression, aes(x = group,y = PRH2,fill= group)) +
geom_boxplot(aes(fill = group),width= 1) +
geom_point(size= 1, alpha=0.6) +
geom_line(aes(group = paired), colour="black", linetype="11",size=0.5)+
scale_fill_manual(values = c("#1CB4B8", "#EB7369"))+
scale_y_continuous(limits=c(0, 25),breaks =seq(0,20,5),expand = c(0, 0) )+
labs(y="log2 Expression",x=NULL,title = "TCGA")+
theme_bw()+ mytheme +
stat_compare_means(
label = "p.format",
method = "t.test",
paired = T,
size = 4,
label.x = 1.5,
label.y = 22
)
box_paired
推荐一款高颜值免费在线SCI绘图工具~~~ 可以直接做配对箱线图
ggsave("./Output/TCGA_PRH2_paired_TumorvsNormal.PDF",
box_paired ,height = 8 ,width = 6,units = "cm")
sessionInfo()
## R version 4.0.2 (2020-06-22)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 14393)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=Chinese (Simplified)_China.936 LC_CTYPE=Chinese (Simplified)_China.936
## [3] LC_MONETARY=Chinese (Simplified)_China.936 LC_NUMERIC=C
## [5] LC_TIME=Chinese (Simplified)_China.936
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] ggthemes_4.2.0 ggpubr_0.4.0 ggplot2_3.3.2
## [4] dplyr_1.0.0 DESeq2_1.28.1 SummarizedExperiment_1.18.2
## [7] DelayedArray_0.14.0 matrixStats_0.56.0 Biobase_2.48.0
## [10] GenomicRanges_1.40.0 GenomeInfoDb_1.24.2 IRanges_2.22.2
## [13] S4Vectors_0.26.1 BiocGenerics_0.34.0
##
## loaded via a namespace (and not attached):
## [1] tidyr_1.1.0 bit64_0.9-7 splines_4.0.2 carData_3.0-4
## [5] blob_1.2.1 cellranger_1.1.0 GenomeInfoDbData_1.2.3 yaml_2.2.1
## [9] pillar_1.4.6 RSQLite_2.2.0 backports_1.1.7 lattice_0.20-41
## [13] glue_1.4.1 digest_0.6.25 RColorBrewer_1.1-2 ggsignif_0.6.0
## [17] XVector_0.28.0 colorspace_1.4-1 htmltools_0.5.0 Matrix_1.2-18
## [21] XML_3.99-0.4 pkgconfig_2.0.3 broom_0.7.0 haven_2.3.1
## [25] genefilter_1.70.0 zlibbioc_1.34.0 purrr_0.3.4 xtable_1.8-4
## [29] scales_1.1.1 openxlsx_4.1.5 rio_0.5.16 BiocParallel_1.22.0
## [33] tibble_3.0.2 annotate_1.66.0 farver_2.0.3 generics_0.1.0
## [37] car_3.0-8 ellipsis_0.3.1 withr_2.2.0 readxl_1.3.1
## [41] survival_3.1-12 magrittr_1.5 crayon_1.3.4 memoise_1.1.0
## [45] evaluate_0.14 rstatix_0.6.0 forcats_0.5.0 foreign_0.8-80
## [49] tools_4.0.2 data.table_1.12.8 hms_0.5.3 lifecycle_0.2.0
## [53] stringr_1.4.0 munsell_0.5.0 locfit_1.5-9.4 zip_2.1.1
## [57] AnnotationDbi_1.50.1 compiler_4.0.2 rlang_0.4.6 grid_4.0.2
## [61] RCurl_1.98-1.2 labeling_0.3 bitops_1.0-6 rmarkdown_2.3
## [65] gtable_0.3.0 abind_1.4-5 DBI_1.1.0 curl_4.3
## [69] R6_2.4.1 knitr_1.29 bit_1.1-15.2 stringi_1.4.6
## [73] Rcpp_1.0.5 vctrs_0.3.1 geneplotter_1.66.0 tidyselect_1.1.0
## [77] xfun_0.15参考资料
- 《利用R代码从UCSC XENA下载mRNA, lncRNA, miRNA表达数据并匹配临床信息》
- 《差异分析|DESeq2完成配对样本的差异分析》
- 《R数据科学》
- 数据和代码下载：
https://gitee.com/ct5869/shengxin-baodian/tree/master/TCGA - 作者：赵法明
编辑：生信宝典
