TCGA_联合GTEx分析2_查看批次效应

在 TCGA_联合GTEx分析1_得到表达矩阵.tpm_老实人谢耳朵的博客-CSDN博客中，获取了TCGA和GTEx中样本的表达矩阵数据，数据格式均为tpm。本文对二者进行合并后，通过PCA分析、绘制内参箱线图等方法，查看是否存在批次效应。

关于批次效应的说明，可参看批次效应（Batch effect）解读

一、数据准备

1 合并后的表达矩阵

exp_tcga.tpm <- read.csv(file = "exp_tcga.tpm.csv", header=T, row.names=1,check.names=FALSE)
exp_gtex.tpm <- read.csv(file = "exp_gtex.tpm.csv", header=T, row.names=1,check.names=FALSE)t_index=intersect(rownames(exp_gtex.tpm),rownames(exp_tcga.tpm))
exp_m = cbind(exp_tcga.tpm[t_index,],exp_gtex.tpm[t_index,])
rm(exp_tcga.tpm,exp_gtex.tpm)

View(exp_tcga.tpm)

包括500个TP，52个NT

View(exp_gtex.tpm)

包括100个NT

View(exp_m)

2 分组信息

group = ifelse(colnames(exp_m[,1:552]) %in% t_dataSmTP,'tcgaTP','tcgaNT')
group = c(group,rep('gtexNT',ncol(exp_m)-length(group)))
barcode=colnames(exp_m)design = data.frame('Barcode'=barcode,'Group'=group)

View(design)

二、主成分分析（Principal Component Analysis）查看批次效应

R语言如何绘制PCA图（四）_心有灵犀啦的博客-CSDN博客_r语言绘制pca图

PCA分析是查看批次效应的最佳方式，如果样本明显按照批次聚类，说明存在批次效应。

library(ggplot2)
library(ggbiplot)

1 数据准备——转置矩阵

exp_m_t=as.data.frame(t(exp_m))

主成分分析，绘制的总是“行变量”的聚类图，因为想看的是barcode的聚类而不是基因的聚类，所以进行转置，使barcode转成行变量。

View(exp_m_t)

2 PCA分析

pca_result <- prcomp(exp_m_t,scale=T)  # 一个逻辑值，指示在进行分析之前是否应该将变量缩放到具有单位方差ggbiplot(pca_result, var.axes=F,            # 是否为变量画箭头obs.scale = 1,         # 横纵比例 groups = design$Group, # 添加分组信息，将按指定的分组信息上色ellipse = T,           # 是否围绕分组画椭圆circle = F) +
ggtitle('PCAplot_tcga&gtex') +
xlim(-100, 200) + ylim(-200, 100) #限制横纵轴范围

PCA图中，tcgaNT 和 gtexNT 明显分为两个亚群，表明存在较强的批次效应。

三、内参表达箱线图查看批次效应

library(ggplot2)
library(reshape2)

1 数据准备——构造 exp_Reshape 用于绘制箱线图

exp_R = melt(as.matrix(log2(exp_m+1)))     #melt()的输入必须为matrix,得到的exp_R为dataframe
colnames(exp_R) = c('Gene','Barcode','Value')
exp_R$Group = rep(group,each=nrow(exp_m))

View(exp_R)

2 ggplot2绘图

gene='All_gene'p_allgene = ggplot(exp_R,aes(x=Group, y=Value,fill=Group)) + #fill进行颜色填充geom_boxplot() +stat_summary(fun="mean",geom="point",shape=23, size=4,fill="white") + #添加均值点ggtitle(paste0('Expression of ',gene)) +xlab('Group') + ylab('log2(tpm+1)') + #x轴y轴标签ylim(0, 5)
p_allgene

所有样本 All_genes 表达量的平均值

gene='ACTB'p_actb = ggplot(exp_R[exp_R$Gene==gene,],aes(x=Group, y=Value,fill=Group)) + #fill进行颜色填充geom_boxplot() +stat_summary(fun="mean",geom="point",shape=23, size=4,fill="white") + #添加均值点ggtitle(paste0('Expression of ',gene)) +xlab('Group') + ylab('log2(tpm+1)') + #x轴y轴标签ylim(10, 15)
p_actb

所有样本 ACTB 表达量的平均值

gene='RPLP0'p_rplp0 = ggplot(exp_R[exp_R$Gene==gene,],aes(x=Group, y=Value,fill=Group)) + #fill进行颜色填充geom_boxplot() +stat_summary(fun="mean",geom="point",shape=23, size=4,fill="white") + #添加均值点ggtitle(paste0('Expression of ',gene)) +xlab('Group') + ylab('log2(tpm+1)') + #x轴y轴标签ylim(7, 12)
p_rplp0

所有样本 RPLP0 表达量的平均值