miRNA分析流程学习(二)/TCGAmiRNA数据三大R包整合差异分析再学习

获得了miRNA之后，我们可以尝试做一下差异分析，那么这种差异分析本质上是于mRNAseq的流程一样的。曾老师/小洁老师也已经在多个推文中展示了mRNAseq的整合差异分析方法。

分析流程：

1、导入

rm(list = ls())
proj = "TCGA-HNSC"
load(paste0(proj,"_miRNA_count.Rdata"))

2、读取和整理数据

2.1 表达矩阵(见流程学习一)

exp <- miRNA_count

2.2 临床信息

load("hnsc_clinical.Rdata")
clinical = hnsc_clinical

3、过滤

# exp1 = exp[rowSums(exp)>0,]
# nrow(exp1)

exp = exp[apply(exp, 1, function(x) sum(x > 0) > 0.5*ncol(exp)), ]
nrow(exp)

4、分组信息获取

# TCGA数据处理方式1
library(tinyarray) # 小洁老师的R包可以快速划分肿瘤/非肿瘤组织
Group = make_tcga_group(exp)
table(Group)
# Group
# normal  tumor 
#     44    525 

5、保存数据

save(exp,Group,proj,clinical,file = paste0(proj,"_miRNApre.Rdata"))
setwd("..")

6、三大R包差异分析

rm(list = ls())

proj = "TCGA-HNSC"
load(paste0(proj,"_miRNApre.Rdata"))
table(Group)

#deseq2----
library(DESeq2)
colData <- data.frame(row.names =colnames(exp), # 创建colData数据框
                      condition=Group)
#注意事项:如果多次运行，表达矩阵改动过的话，需要从工作目录下删除下面if括号里的文件
if(!file.exists(paste0(proj,"_dd.Rdata"))){
  dds <- DESeqDataSetFromMatrix(
  countData = exp,
  colData = colData,
  design = ~ condition)
  dds <- DESeq(dds)
  save(dds,file = paste0(proj,"_dd.Rdata"))
}  
##因为上述步骤中的DESeq是计算差异的函数，数据量大了之后容易出现限速；
##因此为了避免每次加载都要等待，小洁老师用了一个if函数
##如果存在了R.data文件就不需要再加载了
load(file = paste0(proj,"_dd.Rdata"))
class(dds)
res <- results(dds, contrast = c("condition",rev(levels(Group))))
#constrast
c("condition",rev(levels(Group)))
class(res)
DEG1 <- as.data.frame(res)
library(dplyr)
DEG1 <- arrange(DEG1,pvalue)
DEG1 = na.omit(DEG1)
head(DEG1)

#添加change列标记基因上调下调
logFC_t = 1
pvalue_t = 0.05

k1 = (DEG1$pvalue < pvalue_t)&(DEG1$log2FoldChange < -logFC_t);table(k1)
k2 = (DEG1$pvalue < pvalue_t)&(DEG1$log2FoldChange > logFC_t);table(k2)
DEG1$change = ifelse(k1,"DOWN",ifelse(k2,"UP","NOT"))
table(DEG1$change)
head(DEG1)

1. DESeqDataSetFromMatrix：创建 DESeqDataSet 对象。
2. countData = exp：使用基因表达计数矩阵 exp。
3. colData = colData：指定样本的元数据。
4. design = ~ condition：使用 condition（组别）作为差异分析的设计公式。
5. DESeq(dds)：进行差异表达分析。
6. contrast = c("condition", rev(levels(Group)))：指定比较的组别。
7. rev(levels(Group))：获取组别的反向顺序，这样可以设置哪个组为对照组，哪个为处理组。
8. results(): 是 DESeq2 包中的一个函数，用于从 dds 对象中提取差异表达分析的结果。contrast 参数: 用来指定要比较的组别。格式是 c("factor_name", "level1", "level2")，表示 level1 相对于 level2 的变化情况。"condition": 这是用于分组的因子（通常是实验设计中的某个分组变量）。在 colData 中已经定义了这个变量。rev(levels(Group)):levels(Group): 获取分组变量 Group 中的所有级别（例如，Group 的级别可能是 "Control" 和 "Treatment"）。rev(): 反转分组变量的顺序。假设 levels(Group) 的级别是 ("Control", "Treatment")，那么 rev(levels(Group)) 将返回 ("Treatment", "Control")。效果: 使用 rev(levels(Group)) 反转顺序，意味着 contrast 会比较 Group 的最后一个级别相对于第一个级别的差异。比如，如果 levels(Group) 是 ("Control", "Treatment")，那么 rev(levels(Group)) 会变成 ("Treatment", "Control")，表示你想要比较 "Treatment 相对于 Control"。

注意最初输入的Group因子数据中level差异是normal在前，tumor在后面。

#edgeR----
library(edgeR)

# edgeR 中的一个数据结构，用于存储基因表达数据和相关的样本信息。
dge <- DGEList(counts=exp,group=Group)
dge$samples$lib.size <- colSums(dge$counts)
dge <- calcNormFactors(dge) 

design <- model.matrix(~Group)

dge <- estimateGLMCommonDisp(dge, design)
dge <- estimateGLMTrendedDisp(dge, design)
dge <- estimateGLMTagwiseDisp(dge, design)

fit <- glmFit(dge, design)
fit <- glmLRT(fit) 

DEG2=topTags(fit, n=Inf)
class(DEG2)
DEG2=as.data.frame(DEG2)
head(DEG2)

k1 = (DEG2$PValue < pvalue_t)&(DEG2$logFC < -logFC_t)
k2 = (DEG2$PValue < pvalue_t)&(DEG2$logFC > logFC_t)
DEG2$change = ifelse(k1,"DOWN",ifelse(k2,"UP","NOT"))

head(DEG2)
table(DEG2$change)

1. colSums(dge$counts)：计算每个样本的总计数，作为库大小。这一步是为了确保库大小被正确计算和存储在 dge 对象中。
2. calcNormFactors：计算样本之间的标准化因子，用于校正测序深度和样本间差异，确保可以进行跨样本比较。
3. model.matrix：创建一个用于线性模型拟合的设计矩阵。~Group：表示使用 Group 因子作为模型中的解释变量。这个矩阵用于表示不同样本组的比较关系。
4. estimateGLMCommonDisp：估计通用离散度，所有基因使用相同的离散度值。
5. estimateGLMTrendedDisp：估计有趋势的离散度，根据基因的平均表达量来调整离散度。
6. estimateGLMTagwiseDisp：估计每个基因的特异离散度。
7. glmFit：使用设计矩阵拟合广义线性模型，用于计算每个基因的差异表达。
8. glmLRT：基于拟合的 GLM，进行似然比检验（Likelihood Ratio Test, LRT）来检测差异表达基因。

#limma----
library(limma)
dge <- edgeR::DGEList(counts=exp)
dge <- edgeR::calcNormFactors(dge)
design <- model.matrix(~Group)
v <- voom(dge,design, normalize="quantile")

fit <- lmFit(v, design)
fit= eBayes(fit)

DEG3 = topTable(fit, coef=2, n=Inf)
DEG3 = na.omit(DEG3)

k1 = (DEG3$P.Value < pvalue_t)&(DEG3$logFC < -logFC_t)
k2 = (DEG3$P.Value < pvalue_t)&(DEG3$logFC > logFC_t)
DEG3$change = ifelse(k1,"DOWN",ifelse(k2,"UP","NOT"))
table(DEG3$change)
head(DEG3)

1. dge <- edgeR::DGEList(counts=exp)：使用 edgeR 包创建一个 DGEList 对象，exp 是表达矩阵，其中行是基因，列是样本。
2. dge <- edgeR::calcNormFactors(dge)：使用 edgeR 的 calcNormFactors 函数对数据进行归一化，调整不同样本间测序深度的差异。
3. design <- model.matrix(~Group)：创建设计矩阵，Group 是实验设计中的分组信息。这个矩阵用来指示不同样本属于哪些实验组。
4. v <- voom(dge, design, normalize="quantile")：使用 voom 函数对数据进行转换，以适应线性模型分析。voom 会对计数数据进行方差稳定化转换，并计算权重，同时使用量化标准化方法对数据进行归一化。
5. fit <- lmFit(v, design)：使用 lmFit 函数对转换后的表达数据 v 进行线性模型拟合。设计矩阵 design 确定了模型中的因变量和自变量。
6. fit = eBayes(fit)：应用贝叶斯修正，生成统计显著性的估计。eBayes 使用贝叶斯方法来调整方差估计，提高小样本数据的稳定性。

tj = data.frame(deseq2 = as.integer(table(DEG1$change)),
           edgeR = as.integer(table(DEG2$change)),
           limma_voom = as.integer(table(DEG3$change)),
           row.names = c("down","not","up")
          );tj
save(DEG1,DEG2,DEG3,Group,tj,file = paste0(proj,"_DEG.Rdata"))

7、可视化

library(ggplot2)
library(tinyarray)
exp[1:4,1:4]
# cpm 去除文库大小的影响
dat = log2(cpm(exp)+1)
pca.plot = draw_pca(dat,Group, addEllipses = T);pca.plot ##加了addellipase参数取消了椭圆
save(pca.plot,file = paste0(proj,"_pcaplot.Rdata"))

cg1 = rownames(DEG1)[DEG1$change !="NOT"]
cg2 = rownames(DEG2)[DEG2$change !="NOT"]
cg3 = rownames(DEG3)[DEG3$change !="NOT"]

h1 = draw_heatmap(dat[cg1,],Group)
h2 = draw_heatmap(dat[cg2,],Group)
h3 = draw_heatmap(dat[cg3,],Group)

v1 = draw_volcano(DEG1,pkg = 1,logFC_cutoff = logFC_t)
v2 = draw_volcano(DEG2,pkg = 2,logFC_cutoff = logFC_t)
v3 = draw_volcano(DEG3,pkg = 3,logFC_cutoff = logFC_t)

library(patchwork)
(h1 + h2 + h3) / (v1 + v2 + v3) +plot_layout(guides = 'collect') &theme(legend.position = "none")

ggsave(paste0(proj,"_heat_vo.png"),width = 15,height = 10)

UP=function(df){
  rownames(df)[df$change=="UP"]
}
DOWN=function(df){
  rownames(df)[df$change=="DOWN"]
}

up = intersect(intersect(UP(DEG1),UP(DEG2)),UP(DEG3))
down = intersect(intersect(DOWN(DEG1),DOWN(DEG2)),DOWN(DEG3))
dat = log2(cpm(exp)+1)
hp = draw_heatmap(dat[c(up,down),],Group)
##导出上下调基因


#上调、下调基因分别画维恩图
up_genes = list(Deseq2 = UP(DEG1),
          edgeR = UP(DEG2),
          limma = UP(DEG3))

down_genes = list(Deseq2 = DOWN(DEG1),
          edgeR = DOWN(DEG2),
          limma = DOWN(DEG3))

up.plot <- draw_venn(up_genes,"UPgene")
down.plot <- draw_venn(down_genes,"DOWNgene")

#维恩图拼图，终于搞定

library(patchwork)
#up.plot + down.plot
# 拼图
pca.plot + hp+up.plot +down.plot+ plot_layout(guides = "collect")
ggsave(paste0(proj,"_heat_ve_pca.pdf"),width = 15,height = 10)

既往推文：

miRNA分析流程学习(一)/TCGAmiRNA数据下载

参考资料：

1. 生信技能树时间线：https://mp.weixin.qq.com/mp/appmsgalbum?action=getalbum&__biz=MzAxMDkxODM1Ng==&scene=24&album_id=2201138830328528899&count=3&uin=&key=&devicetype=iMac+Mac14%2C7+OSX+OSX+14.6.1+build(23G93)&version=13080810&lang=zh_CN&nettype=WIFI&ascene=0&fontScale=100
2. 生信技能树B站视频：https://www.bilibili.com/video/BV1zK411n7qr/?vd_source=3a13860df939bc922ad1fd6099e42c1d
3. 生信技能树/小洁老师：https://mp.weixin.qq.com/s/zTN32UA6Nu6PMmqXYTOnnQ

致谢：感谢曾老师/小洁老师以及生信技能树团队全体成员。

注：若对内容有疑惑或者有发现明确错误的朋友，请联系后台(欢迎交流)。更多内容可关注公众号：生信方舟