3.差异分析

上一步我们得到了基因在肝细胞癌样本中的表达矩阵，我这里下载了Counts值和TPM值的表达矩阵，现在使用Counts值进行差异分析：
一、读取表达矩阵，同时根据需求分组，这里分为两组—肝细胞癌组织和正常组织，这里表达矩阵的列名是全名（TCGA-G3-A25X-01A-11R-A16W-07），我们可以只保留前面的字符（TCGA-G3-A25X-01A），方便分组（“-01A”表示肿瘤，“-11A”表示正常组织 ），然后合并矩阵并构建分组。

rm(list = ls())

# 安装加载需要的R包
library(pacman)
p_load(data.table, tidyverse, edgeR, ggbeeswarm, ggsignif)
# 设置工作路径
setwd('D:/Bioinformatics/TCGA/差异分析/')
 
# 1、数据准备 
## 1.1 读取表达矩阵，将id列转化为行名。
exp_data = fread("LIHC_Portal_RNA_counts.txt", h = T, sep = "\t", check.names = F) %>% column_to_rownames("id")
exp_data[1:4,1:4] 
# 列名处理，保留前16个字符。 
colnames(exp_data) <- substr(colnames(exp_data), 1, 16)
# 这里“-01A”表示肿瘤，“-11A”表示正常组织 
exp_data_T = exp_data %>% dplyr::select(str_which(colnames(.), "-01A$")) # 匹配列名或用下示写法
nT = ncol(exp_data_T) # 369
exp_data_N = exp_data %>% dplyr::select(ends_with("-11A"))
nN = ncol(exp_data_N) # 50
exp_data_N[1:3,1:3]

## 1.2 合并表达数据
exp_data_final = cbind(exp_data_N, exp_data_T)
dim(exp_data_final) 
## 1.3 建立分组
SampleLabel = factor(c(rep("NT", nN), rep("HCC", nT)), levels = c("NT", "HCC"))

二、开始差异分析，这里有一个过滤基因的操作【留下的基因应该至少在50个样本（较少组的样本量）中的表达量大于0.5】，但并不是绝对的。

# 2、edgeR差异表达分析 
# 构建DGEList对象
dge = DGEList(counts = exp_data_final, group = SampleLabel)

# 过滤基因
keep = rowSums(cpm(dge) > 0.5) > nN # 至少在较少组样本量中的表达量cpm值大于0.5
dge_filter = dge[keep, keep.lib.sizes = FALSE]
dim(dge_filter$counts)

dge_norm = calcNormFactors(dge_filter)  # 默认 method = "TMM"

# 获取logCPM,用于可视化展示
logcpm = cpm(dge_norm, log = TRUE, prior.count = 1)

design = model.matrix(~SampleLabel)
y = estimateDisp(dge_norm, design, robust = TRUE)

# 得到差异分析结果 
fit = glmQLFit(y, design)
res = glmQLFTest(fit)
result = topTags(res, n = Inf)$table

# 保存结果 
write_csv(data.frame(Symbol = rownames(result), result, check.names = F), "RNA_DEA_result_info.csv")
# 筛选出显著差异的基因 并保存。 
sig_result = subset(result, abs(logFC) > 1 & FDR < 0.05)
write_csv(data.frame(Symbol = rownames(sig_result), sig_result, check.names = F), "DEA_result_logFC1_fdr0.05.csv")

三、绘制目标基因的差异表达图（箱线图、蜂群图、小提琴图等），这里绘制蜂群图（如下图1）。

#### beeswarm plot（蜂群图） ----
target = 'TK1'
# 绘图数据
DataForPlot = data.frame(Exp = logcpm[target,], Group = SampleLabel)
# 设置y轴范围
ymax = max(DataForPlot$Exp)
ymin = min(DataForPlot$Exp)
# 获取p.adj值
fdrvalue = result[target,5]
# 绘图
p = ggplot(DataForPlot, aes(x = Group, y = Exp, col = Group)) +
    geom_beeswarm(size = 1, cex = 1.5, alpha = 0.5) +  # cex设置分散的宽度，alpha设置点的透明度
    ylim(1.1*ymin, 1.3*ymax) + 
    labs(x = "", y = bquote(italic(.(target))~expression~(log[2]))) +
    scale_x_discrete(labels = paste(levels(SampleLabel), "\n", "(N=", table(SampleLabel), ")", sep = "")) +
    scale_colour_manual(values = c("#2700F6", "#C8181C")) + 
    stat_summary(fun = median, geom = "point", size = 1, col = "black") +  # median处画个点
    stat_summary(fun = median, fun.min = median, fun.max = median, geom = "crossbar", size = 0.3, width = 0.4, col = "black") +  # median处加条横线，注意三个参数都要为median!
    stat_summary(fun.data = function(x) median_hilow(x, 0.5), geom = "errorbar", width = 0.25, col = "black") +  # median_hilow(x, 0.5) 求中值和上下四分位，errorbar绘制"工"字型线
    theme_classic() + theme(legend.position = "none") +
    theme(plot.margin = unit(c(0.2, 0.2, -0.3, 0.2), "cm"))
# 添加显著性p值,注意修改分组名称
p + geom_signif(comparisons = list(c("NT", "HCC")), y_position = 1.15*ymax, textsize = 3, annotation = paste0("list(~italic(p.adj)==", signif(fdrvalue, 3), ")"), parse = T, color = "black") 
ggsave(paste0(targe, "_beeswarm.pdf"), width = 3, height =4)

差异分析.png