DESeq2分析RNAseq数据

生信喵

背景
       DESeq2 主要需要两个数据，一个是表达矩阵，该表达矩阵必须为计数数据，而不能是归一化之后的数据，也就是必须为 reads 数，而不能是 FPKM，TPM 等值。另外一个是样品信息，主要提供分组信息。
       可以直接将 salmon 或者 hisat2 的比对结果读入 DESeq2 中。

一、读取表达矩阵

1. library(DESeq2)
   
2. countdata <- read.csv("CountMatrix.csv",header = T,row.names = 1)
   
3. head(countdata, 10)
   
4. coldata <- read.csv("sample.csv",header = T)
   
5. #关键步骤，生成Deseq对象
   
6. dds <- DESeqDataSetFromMatrix(countData = countdata,colData = coldata,design = ~ cell + dex)

复制代码

      上述代码中的count和分组，我提供一个压缩包，可以关注公众号回复"20221009"获取。
      

二、SummarizedExperiment 类
       SummarizedExperiment 类用于存储测序或者芯片试验的结果。每个对象中包含一个或者多个样本的观察值，以及描述观察值（observations），即特征（features） 和样本（samples），即表型（phenotypes）的元数据（meta-data）。
       SummarizedExperiment 类的一个重要特征是元数据与实验数据的协调，例如想要剔除某个样本，那么单个操作即可完成对元数据和试验的处理，确保元数据和观察值保持同步。
       SummarizedExperiment 在其行信息更加灵活，可以是基于 GRanges 或者任意的 DataFrame。
2.1 解析 SummarizedExperiment
       SummarizedExperiment 包中包含两个类：SummarizedExperiment 和RangedSummarizedExperiment。
       SummarizedExperiment：SummarizedExperiment 是一个矩阵样的容器，每一行代表感兴趣的特征（feature），如基因，转录本，外显子等，每一列则代表样本（sample）。每个对象可以包含一个或者多个实验（assay），而它们每一个都用类矩阵的对象表示。
       特征相关的信息存储在一个数据框（DataFrame）对象中，可以通过函数 rowData() 来获取。数据框的每一行提供了 SummarizedExperiment 对象
中对应每一行特征的信息，数据框中的每一列则提供了感兴趣的特征属性，例如基因或者转录本的 ID 等。
       RangedSummarizedExperiment：RangedSummarizedExperiment 是SummarizedExperiment 的子类，所以后者所有的方法在RangedSummarizedExperiment 中也适用。
       两者最主要的区别在于 RangedSummarizedExperiment 中的每一行均为 genomic ranges，而不是包含特征的数据框。这些 genomic ranges 通过 GRanges 或GRangesList 对象描述，可通过 rowRanges()函数来获取。

1. http://www.bioconductor.org/packages/release/bioc/html/SummarizedExperiment.html

复制代码

     

2.2 SummarizedExperiment 类与 DESeqDataSet 类
       1 后者使用 assay 函数获取矩阵数据
       2 后者数据必须为非负数
       3 后者包含分组比较公式
      
       SE 类

2.3 处理 SE 类常用函数

1. dim(dds)
   
2. assay(dds)
   
3. assayNames(dds)
   
4. colSums(assay(dds))
   
5. rowRanges(dds)
   
6. colData(dds)
   
7. 
   
8. #过滤没有reads比对上的基因，所有reads数为零
   
9. nrow(dds)
   
10. dds <- dds[rowSums(counts(dds)) > 1,]
    
11. nrow(dds)
    
12. colSums(assay(dds))
    
13. 
    
14. par('mar')
    
15. par(mar=c(7.1,4.1,4.1,2.1))
    
16. barplot(colSums(assay(dds)),las=3,col=rep(rainbow(4),each=2))

复制代码

三、多维数据探索
       将数据导入 R 中并不能直接用于基因差异表达的计算。首先要对数据进行质控，看组内差别是否小于组间差别。一个好的 RNAseq 实验，要求组内差别越小越好，组间差别越大越好。如果组内差别已经大于组间差别，则该实验结果没有意义。多维数据探索通过聚类分析，PCA，MDS 等方法，可以检测样品之间的相互关系。
3.1 数据转换

1. # 将数据通过rlog方法与vst方法转换，这样可以用于后面计算距离矩阵
   
2. ## ----rlog方法----------------------------------------------------------------
   
3. rld <- rlog(dds, blind = FALSE)
   
4. head(assay(rld), 3)
   
5. 
   
6. ## ----vst方法-----------------------------------------------------------------
   
7. vsd <- vst(dds, blind = FALSE)
   
8. head(assay(vsd), 3)

复制代码

3.2 相关性热图
      
       欧式距离相关性热图

1. #利用转换后的结果计算样品之间距离关系
   
2. #方法1--欧氏距离--------------------------
   
3. sampleDists <- dist(t(assay(rld)))
   
4. sampleDists
   
5. 
   
6. library("pheatmap")
   
7. library("RColorBrewer")
   
8. 
   
9. ## ----distheatmap, fig.width = 6.1, fig.height = 4.5----------------------
   
10. sampleDistMatrix <- as.matrix( sampleDists )
    
11. rownames(sampleDistMatrix) <- paste( rld$dex, rld$cell, sep = " - " )
    
12. colnames(sampleDistMatrix) <- NULL
    
13. pheatmap(sampleDistMatrix,clustering_distance_rows = sampleDists,
    
14.          clustering_distance_cols = sampleDists)
    
15. 
    
16. ## 方法2，使用Poisson Distance方法------------------------------------------
    
17. library("PoiClaClu")
    
18. poisd <- PoissonDistance(t(counts(dds)))
    
19. samplePoisDistMatrix <- as.matrix( poisd$dd )
    
20. rownames(samplePoisDistMatrix) <- paste( rld$dex, rld$cell, sep=" - " )
    
21. colnames(samplePoisDistMatrix) <- NULL
    
22. pheatmap(samplePoisDistMatrix,clustering_distance_rows = poisd$dd,
    
23.          clustering_distance_cols = poisd$dd)

复制代码

3.3 PCA
      
        PCA 图

1. plotPCA(rld, intgroup = c("dex", "cell"))
   
2. 
   
3. ## ------------------------------------------------------------------------
   
4. pcaData <- plotPCA(rld, intgroup = c( "dex", "cell"), returnData = TRUE)
   
5. pcaData
   
6. percentVar <- round(100 * attr(pcaData, "percentVar"))
   
7. library(ggplot2)
   
8. ggplot(pcaData, aes(x = PC1, y = PC2, color = dex, shape = cell)) +
   
9.   geom_point(size =3) +  xlab(paste0("PC1: ", percentVar[1], "% variance")) +
   
10.   ylab(paste0("PC2: ", percentVar[2], "% variance")) +  coord_fixed()

复制代码

3.4 MDS
      

1. library(dplyr)
   
2. mds <- as.data.frame(colData(rld))  %>%   cbind(cmdscale(sampleDistMatrix))
   
3. ggplot(mds, aes(x = `1`, y = `2`, color = dex, shape = cell)) +  
   
4.   geom_point(size = 3) + coord_fixed()
   
5. mdsPois <- as.data.frame(colData(dds)) %>%
   
6.   cbind(cmdscale(samplePoisDistMatrix))
   
7. ggplot(mdsPois, aes(x = `1`, y = `2`, color = dex, shape = cell)) +
   
8.   geom_point(size = 3) + coord_fixed()

复制代码

四、差异表达基因筛选
       FoldChange：就是两样品中同一个基因表达水平的变化倍数。一般取 log2 FoldChange（lfc），这样差异倍数小于 1 倍的，也就是低表达量的值可以为负值，便于区分。
       pvalue： Probability，Pr，统计学根据显著性检验方法所得到的 P 值，一般以 P < 0.05 为显著， P <0.01 为非常显著，其含义是样本间的差异由抽样误差所致的概率小于 0.05 或 0.01。
       FDR：是 False Discovery Rate control 的简称。 就是多重检验校正中一种校正 p-value 的统计学方法。在实际中，FDR 是错误发现率的期望。
       Qvalue：衡量错误发现率的指标（False discovery rate，简称 FDR）, 由于Qvalue 需要利用公式从 P value 校正计算后得到，所以 Q value 通常又被称为 adjusted p value。一般来说，Q value = FDR = adjusted p value。

1. dep <- DESeq(dds)
   
2. res <- results(dep)
   
3. res
   
4. write.csv(x = res,file = "des.csv")
   
5. 
   
6. #筛选出差异表达基因，log2foldchange >=1或者<=-1,并且q值小于0.05的基因
   
7. res <- na.omit(res) #17264
   
8. res0.05 <- res[(abs(res$log2FoldChange)>=1 & res$padj <=0.05), ]
   
9. nrow(res0.05) #1035
   
10. 
    
11. #筛选出差异表达基因，log2foldchange >=2,并且q值小于0.01的基因
    
12. res0.01 <- res[(abs(res$log2FoldChange)>=1 & res$padj <=0.01), ]
    
13. nrow(res0.01) #878
    
14. 
    
15. write.csv(res0.05,file = "res005.csv")
    
16. #找出差异表达最显著的基因，按log2差异倍数排序
    
17. head(res0.05[order(res0.05$log2FoldChange,decreasing = TRUE), ])
    
18. head(res0.05[order(res0.05$log2FoldChange,decreasing = FALSE), ])

复制代码

五、结果可视化
       火山图用来可视化基因差异表达的结果，输入数据为基因表达差异分析的结果，利用log2foldchange 与 p 值进行绘图。
      
       差异表达基因火山图

1. #MA图
   
2. plotMA(res, ylim = c(-10,10))
   
3. 
   
4. #火山图
   
5. m <- res
   
6. head(m)
   
7. m <- na.omit(m)
   
8. plot(m$log2FoldChange,m$padj)
   
9. plot(m$log2FoldChange,-1*log10(m$padj))
   
10. plot(m$log2FoldChange,-1*log10(m$padj),xlim = c(-10,10),ylim=c(0,100))
    
11. m <- transform(m,padj=-1*log10(m$padj))
    
12. down <- m[m$log2FoldChange<=-1,]
    
13. up <- m[m$log2FoldChange>=1,]
    
14. no <- m[m$log2FoldChange>-1 & m$log2FoldChange <1,]
    
15. 
    
16. plot(no$log2FoldChange,no$padj,xlim = c(-10,10),ylim=c(0,100),col="blue",pch=16,
    
17.      cex=0.8,main = "Gene Expression",
    
18.      xlab = "log2FoldChange",ylab="-log10(Qvalue)")
    
19. 
    
20. points(up$log2FoldChange,up$padj,col="red",pch=16,cex=0.8)
    
21. points(down$log2FoldChange,down$padj,col="green",pch=16,cex=0.8)
    
22. abline(v = c(-1,1),h = 2)

复制代码

六、处理批次效应
       所谓批次效应（batch effect）主要是因为 RNAseq 鉴定的是基因表达动态变化的一个快照。与 DNA 测序不同，RNAseq 是实时变化的。因此，不同批次的样本之间往往会有很大的差别。这种差别甚至可以大到大于组间差别，直接导致整个实验和最终的结论失败。
       批次效应主要是因为不同平台的数据，同一平台的不同时期的数据，同一个样品不同试剂的数据，以及同一个样品不同时间的数据等等原因。
       如何检测是否存在这种效应呢？
       最简单的就是记录实验中时间这个变量，然后对差异表达的基因进行聚类，看是否都和时间相关，如果相关就证明存在 batch effect。

生信喵

1. library("sva")
   
2. dat <- counts(dds, normalized = FALSE)#TRUE
   
3. idx <- rowMeans(dat) > 1
   
4. dat <- dat[idx, ]
   
5. mod <- model.matrix(~ dex, colData(dds))
   
6. mod0 <- model.matrix(~ 1, colData(dds))
   
7. svseq <- svaseq(dat, mod, mod0, n.sv = 2)
   
8. svseq$sv
   
9. par(mfrow = c(2, 1), mar = c(3,5,3,1))
   
10. for (i in 1:2) {
    
11. stripchart(svseq$sv[, i] ~ dds$cell, vertical = TRUE, main = paste0("SV",
    
12. i))
    
13. abline(h = 0)
    
14. }

复制代码

1. ddssva <- dds
   
2. ddssva$SV1 <- svseq$sv[,1]
   
3. ddssva$SV2 <- svseq$sv[,2]
   
4. design(ddssva) <- ~ SV1 + SV2 + dex

复制代码

      目前 RNA-seq 分析去除批次效应的工具主要有 R 包 sva 和 RUVSeq。去除批次效应主要针对两种情况：（1）已知批次变量的，（2）未知批次变量的。第一种情况通常用 sva 的 combat 和 limma 的 removeBatchEffect 功能进行去除。第二种情况可以用 sva 的 svaseq 功能和 RUVSeq 的 RUV 功能进行矫正；这并不是真正的去除，只是在差异基因筛选的时候，根据模型考虑了批次效应的影响。