数据挖掘期末project-乳腺癌预后模型建立-2.1 筛选差异表达基因

1.讲到了在TCGA下载数据，并且得到了表达矩阵，现在开始2. 筛选基因

先是第一步：筛选差异表达基因

方法介绍：

参考链接：第四课：RNA-Seq数据分析——三种主流的差异基因分析（DESeq2，EdgeR，limma） - 知乎 (zhihu.com)

筛选差异表达基因有三种方法：limma、DESeq2、EdgeR

这三种方法都是用于筛选差异表达基因的常用工具，但它们在原理和实现上有所不同：

1. **limma（Linear Models for Microarray Analysis）**：
   - **原理简介：** limma最初是为微阵列数据设计的，但后来也被应用于RNA测序数据分析。它基于线性模型，使用的是贝叶斯推断的方法来鉴定基因的差异表达。
   - **工作方式：** limma通过拟合线性模型来考虑每个基因的表达量之间的变异，并使用经验贝叶斯方法来为每个基因的差异表达进行统计检验。
   - **优点：** 对于小样本RNA测序数据具有较高的灵敏度和准确性。

2. **DESeq2**：
   - **原理简介：** DESeq2是专门针对RNA测序数据设计的工具，采用负二项分布模型来估计基因表达水平的离散性，并鉴定差异表达基因。
   - **工作方式：** DESeq2使用负二项分布来建模RNA测序计数数据的离散性，并基于该模型进行差异表达分析。它考虑了基因表达量的变异性和数据的技术重复性。
   - **优点：** 在小样本情况下表现稳健，适用于处理RNA测序数据的离散计数。

3. **EdgeR**：
   - **原理简介：** EdgeR也是用于RNA测序数据的分析工具，它基于负二项分布模型来估计基因表达的离散性，并检测差异表达基因。
   - **工作方式：** EdgeR利用负二项分布来建模RNA测序数据的离散性，并使用经过修正的方法来进行差异表达分析。它特别关注低表达基因和小样本数据的特性。
   - **优点：** 对于小样本和低表达基因的RNA测序数据具有较好的表现。

RNA-seq数据分析 09：DESeq2差异表达分析 - 知乎 (zhihu.com)

DESeq2是一个为高维计量数据的归一化、可视化和差异表达分析而设计的一个R语言包。它通过经验贝叶斯方法(empirical Bayes techniques)来估计对数倍数变化(log2foldchange）和离差的先验值，并计算这些统计量的后验值。

这里选用的是DESeq2，R语言，参考的代码依然是上一篇的TCGA数据挖掘实战（一）——数据下载和差异分析

代码：

先剔除不表达/表达量为0的基因

countData <- exprSet[rowMeans(exprSet)>1,]        # 去除表达量过低的基因

剔除后：

使用DESeq2在疾病和正常样本之间筛选差异表达基因。

根据TCGA barcode中的sample对样本进行分类，其中01-09代表肿瘤类型，10-19代表正常类型。

rm(list = ls())

library(DESeq2)

exprSet <- read.table("data/COAD_expr.txt", header = T, check.names = F)

sample <- ifelse(substring(colnames(exprSet),14,15)<10, "tumor", "normal")
sample <- factor(sample, levels = c("tumor", "normal"))
group <- data.frame(barcode = colnames(exprSet), sample) 

dds <- DESeqDataSetFromMatrix(countData = round(exprSet),
                              colData = group,
                              design = ~sample)
dds <- DESeq(dds)

res <- results(dds, contrast = c("sample", "tumor", "normal"))
sel_res <- subset(res, padj < 0.05 & abs(log2FoldChange) > 1)
### sel_res是在res的基础上做了一个剔除，设置阙值为0.05和1

write.csv(sel_res, file = "data/DEGs.csv")
write.table(rownames(sel_res), file = "data/DEGs_gene_name.txt", quote = F,
             row.names = F, col.names = F)

小L生信学习日记-6丨你最关心的差异基因是怎么挑出来的？！ - 知乎 (zhihu.com)

前面几个列头代表的意思很容易理解，我们最关注的是log2FoldChange、P值和Padjust值

P-value：p值，是统计学检验变量，代表差异显著性，一般认为P < 0.05 为显著， P <0.01 为非常显著。其含义为：由抽样误差导致样本间差异的概率小于0.05 或0.01。

Padj：p.adjust。转录组测序的差异表达分析是对大量的基因表达值进行的独立统计假设检验，存在假阳性问题，因此引入Padj对显著性P值（P.adjust）进行校正。Padj是对P-value的再判断，筛选更为严格。

Fold Change：适用于两分组分析，表示两样本（组）间表达量的比值。

log2FoldChange：对Fold Change取log2，一般默认表达相差2倍以上是有意义的，可以根据情况适当放宽至1.5/1.2，但最好不要低于1.2倍（代码这里写的是>1，是否会存在问题？？？）

图来自： 小L生信学习日记-6丨你最关心的差异基因是怎么挑出来的？！ - 知乎 (zhihu.com)

得到的结果：

> table(substring(colnames(exprSet),14,15))

  01   06   11 
1097    7  113 


> table(group$sample)

 tumor normal 
  1104    113 


> head(dds)
class: DESeqDataSet 
dim: 6 1217 
metadata(1): version
assays(6): counts mu ... replaceCounts replaceCooks
rownames(6): MT-CO1 MT-ND4 ... MT-CO3 ACTB
rowData names(23): baseMean baseVar ... maxCooks replace
colnames(1217): TCGA-E9-A1NI-01A TCGA-A1-A0SP-01A ...
  TCGA-B6-A0RN-01A TCGA-BH-A203-01A
colData names(4): barcode sample sizeFactor replaceable


> summary(res)

out of 56486 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up)       : 18320, 32%
LFC < 0 (down)     : 11354, 20%
outliers [1]       : 0, 0%
low counts [2]     : 8746, 15%
(mean count < 0)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results


> head(sel_res)
log2 fold change (MLE): sample tumor vs normal 
Wald test p-value: sample tumor vs normal 
DataFrame with 6 rows and 6 columns
        baseMean log2FoldChange     lfcSE      stat       pvalue         padj
       <numeric>      <numeric> <numeric> <numeric>    <numeric>    <numeric>
COL1A1  395270.7        2.25520 0.1273349   17.7108  3.46219e-70  1.35702e-68
COL1A2  287915.9        1.57193 0.1243126   12.6450  1.19268e-36  1.59716e-35
FN1     273233.3        2.85107 0.1222278   23.3259 2.42064e-120 3.56799e-118
COL3A1  261230.3        1.38620 0.1264941   10.9586  6.04245e-28  5.65842e-27
GAPDH   129513.8        1.18825 0.0838662   14.1684  1.43809e-45  2.65587e-44
AHNAK    99416.3       -1.51978 0.0819419  -18.5471  8.60897e-77  4.19808e-75

> summary(sel_res)

out of 11298 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up)       : 7236, 64%
LFC < 0 (down)     : 4062, 36%
outliers [1]       : 0, 0%
low counts [2]     : 0, 0%
(mean count < 0)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

差异倍数大于1.5且p值小于0.05的为上调基因， 用红色点展示；差异倍数小于0.666667并且p值小于0.05的为下调基因， 用绿色点展示；其中非显著差异的基因用灰色点展示。

上面，我们提到了上调基因、下调基因。那么，什么是上调基因，下调基因？简单解释一下：

• 上调基因：up-regulated gene，相对于对照组，在实验组中该基因转录成mRNA时受到正向调控，表达量增加。
• 下调基因：down-regulatedgene，相对于对照组，在实验组中该基因转录成mRNA时受到抑制，表达量减少

绘制火山图：RNA-seq数据分析 09：DESeq2差异表达分析 - 知乎 (zhihu.com)

genes<- res1
# 根据上调、下调、不变为基因添加颜色信息
genes$color <- ifelse(genes$padj<0.05 & abs(genes$log2FoldChange)>= 1,ifelse(genes$log2FoldChange > 1,'red','blue'),'gray')
color <- c(red = "red",gray = "gray",blue = "blue")

p <- ggplot(
  # 指定数据、映射、颜色
  genes, aes(log2FoldChange, -log10(padj), col = color)) +  
  geom_point() +
  theme_bw() +
  scale_color_manual(values = color) +
  # 辅助线
  labs(x="log2 (fold change)",y="-log10 (q-value)") +
  geom_hline(yintercept = -log10(0.05), lty=4,col="grey",lwd=0.6) +
  geom_vline(xintercept = c(-1, 1), lty=4,col="grey",lwd=0.6) +
  # 图例
  theme(legend.position = "none",
        panel.grid=element_blank(),
        axis.title = element_text(size = 16),
        axis.text = element_text(size = 14)+
  # 注释
  geom_text_repel(
    data = subset(genes, padj < 1e-100 & abs(genes$log2FoldChange) >= 10),
    aes(label = rownames(genes)),
    size = 5,
    box.padding = unit(0.35, "lines"),
    point.padding = unit(0.3, "lines"))
        )
p