感谢老俊俊的大力支持。我们会每日跟新，欢迎您关注老俊俊的生信笔记。

Part1RNA-seq：Salmon 定量结果差异分析

上篇文章讲到如何使用 salmon 快速进行定量分析，这回讲解如何将定量结果导入到 DESeq2 软件中进行差异分析。
Salmon 下游分析主要有：

differential transcript usage (DTU) 差异转录本的使用
differential gene expression (DGE) 差异基因表达分析

这里主要介绍差异基因表达分析，DTU 分析可参考：Swimming downstream: statistical analysis of differential transcript usage following Salmon quantification^[1] 或者 https://f1000research.com/articles/7-952/v3^[2]

1安装并加载 R 包

BiocManager::install("tximport")
BiocManager::install("GenomicFeatures")
library(tximport)
library(GenomicFeatures)

2构造输入文件

# 设置定量结果存放目录
setwd('D:/rnaseq/salmon/quants')
# 获取当前工作路径
dir <- getwd()
# 保存当前文件夹名给file
file <- c('test1_rep1_quant','test1_rep2_quant',
          'test2_rep1_quant','test2_rep2_quant')
# 或者
file <- list.files(pattern = '*.quant')

# 构建定量结果文件路径
file_dir <- file.path(dir, file, 'quant.sf')
file_dir
[1] "D:/rnaseq/salmon/quants/test1_rep1_quant/quant.sf"
[2] "D:/rnaseq/salmon/quants/test1_rep2_quant/quant.sf"
[3] "D:/rnaseq/salmon/quants/test2_rep1_quant/quant.sf"
[4] "D:/rnaseq/salmon/quants/test2_rep2_quant/quant.sf"
# 给文件命名
names(file_dir) <- c('test1_rep1','test1_rep2','test2_rep1','test2_rep2')
file_dir
                                         test1_rep1                                          test1_rep2
"D:/rnaseq/salmon/quants/test1_rep1_quant/quant.sf" "D:/rnaseq/salmon/quants/test1_rep2_quant/quant.sf"
                                         test2_rep1                                          test2_rep2
"D:/rnaseq/salmon/quants/test2_rep1_quant/quant.sf" "D:/rnaseq/salmon/quants/test2_rep2_quant/quant.sf"

3准备转录本 id 和基因 id 对应关系文件

可以从 TxDb 包中获得，也可以自己从 gtf 注释文件中获得，推荐后者，尽量和上游分析注释文件保持一致，注释文件 R 包都是 UCSC 的注释类型，可能注释完整性每个数据库也不一样。

# 从自己gtf文件获取，把gtf文件放到当前工作目录
# 读取gtf文件
txdb <- makeTxDbFromGFF(file = 'gencode.vM27.annotation.gtf',format = 'gtf',organism = 'Mus musculus')

# 查看keys类型
keytypes(txdb)
[1] "CDSID"    "CDSNAME"  "EXONID"   "EXONNAME" "GENEID"   "TXID"     "TXNAME"

# 选择GENEID和TXNAME
k <- keys(txdb, keytype = "TXNAME")
tx2gene <- select(txdb, k, "GENEID", "TXNAME")

# 查看内容
head(tx2gene,4)
                TXNAME               GENEID
1 ENSMUST00000193812.2 ENSMUSG00000102693.2
2 ENSMUST00000082908.3 ENSMUSG00000064842.3
3 ENSMUST00000192857.2 ENSMUSG00000102851.2
4 ENSMUST00000161581.2 ENSMUSG00000089699.2

# 从注释文件R包中获取
# 人hg19注释信息
BiocManager::install("TxDb.Hsapiens.UCSC.hg19.knownGene")
# 人hg38注释信息
BiocManager::install("TxDb.Hsapiens.UCSC.hg38.knownGene")
# 小鼠mm9注释信息
BiocManager::install("TxDb.Mmusculus.UCSC.mm9.knownGene")
# 小鼠mm10注释信息
BiocManager::install("TxDb.Mmusculus.UCSC.mm10.knownGene")
##########################################################
library(TxDb.Mmusculus.UCSC.mm10.knownGene)
txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene
keytypes(txdb)
k <- keys(txdb, keytype = "TXNAME")
tx2gene <- select(txdb, k, "GENEID", "TXNAME")
head(tx2gene,4)

4使用 tximport^[3] 导入

tximport 包能导入以下软件的定量结果：

Salmon (Patro et al. 2017)
Alevin (Srivastava et al. 2019)
Sailfish (Patro, Mount, and Kingsford 2014)
kallisto (Bray et al. 2016)
RSEM (Li and Dewey 2011)
StringTie (Pertea et al. 2015)

具体可参考 tximport 的 bioconductor 帮助文档。

# 读取定量结果,保存为一个list
txi <- tximport(file_dir, type = "salmon", tx2gene = tx2gene)
reading in files with read_tsv
1 2 3 4
summarizing abundance
summarizing counts
summarizing length

# 查看list
names(txi)
[1] "abundance"   "counts"    "length"    "countsFromAbundance"
# 查看abundance数据，代表tpm值
head(txi$abundance,3)
                      test1_rep1 test1_rep2 test2_rep1 test2_rep2
ENSMUSG00000000001.5    62.21182   55.25616   53.45963   60.46954
ENSMUSG00000000003.16    0.00000    0.00000    0.00000    0.00000
ENSMUSG00000000028.16   52.10330   39.86649   35.43731   37.04251
# 查看counts数据
head(txi$counts,3)
                      test1_rep1 test1_rep2 test2_rep1 test2_rep2
ENSMUSG00000000001.5    2695.758   2605.763   2311.899   2838.876
ENSMUSG00000000003.16      0.000      0.000      0.000      0.000
ENSMUSG00000000028.16   1332.969   1139.000    919.746   1047.000
# 查看有效基因长度
head(txi$length,3)
                      test1_rep1 test1_rep2 test2_rep1 test2_rep2
ENSMUSG00000000001.5    2962.254   2998.074   2957.301   2981.350
ENSMUSG00000000003.16    512.772    512.772    512.772    512.772
ENSMUSG00000000028.16   1748.917   1816.368   1774.842   1794.938

Part2差异分析

后续可接三款主流的差异分析软件，下面分别介绍：

1、DESeq2 差异分析代码

# 加载包
library(DESeq2)
# 构建样本信息表
sampleTable <- data.frame(condition = factor(rep(c("control", "treat"), each = 2)))
rownames(sampleTable) <- colnames(txi$counts)

# 使用DESeqDataSetFromTximport函数传入给DESeq2
dds <- DESeqDataSetFromTximport(txi, sampleTable, ~condition)
using counts and average transcript lengths from tximport

# 差异分析
dds <- DESeq(dds)
estimating size factors
using 'avgTxLength' from assays(dds), correcting for library size
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
res <- results(dds)
# 把res转为data.frame格式
diff_res <- as.data.frame(res)
# 查看差异结果
head(diff_res,3)
                      baseMean log2FoldChange     lfcSE       stat    pvalue      padj
ENSMUSG00000000001.5  2610.399    -0.02393157 0.2183073 -0.1096233 0.9127081 0.9673216
ENSMUSG00000000003.16    0.000             NA        NA         NA        NA        NA
ENSMUSG00000000028.16 1110.566    -0.32164574 0.2437125 -1.3197754 0.1869100 0.4593252

# 导出保存文件
diff_res$gene_id <- rownames(diff_res)
write.table(diff_res,file = 'DESeq2_results.csv',sep = ',',row.names = FALSE)

2、edgeR 差异分析代码

# 加载包
library(edgeR)
# 取出counts
cts <- txi$counts
# 取出基因有效长度
normMat <- txi$length

# 得到每个基因长度的缩放因子
normMat <- normMat/exp(rowMeans(log(normMat)))
normCts <- cts/normMat
# 从缩放的counts来计算有效文库大小，避免样本间的组成偏差
eff.lib <- calcNormFactors(normCts) * colSums(normCts)

# Combining effective library sizes with the length factors, and calculating offsets for a log-link GLM.
# 结合有效文库大小和长度因子，用广义线性模型（GLM)计算偏移量
normMat <- sweep(normMat, 2, eff.lib, "*")
normMat <- log(normMat)

# 构建 DGEList object
y <- DGEList(cts)
y <- scaleOffset(y, normMat)
# 自动过滤低表达基因
keep <- filterByExpr(y)
y <- y[keep, ]

# 设计实验
sampleTable <- data.frame(condition = factor(rep(c("control", "treat"), each = 2)))
rownames(sampleTable) <- colnames(txi$counts)

design_full <- model.matrix(~condition, data = sampleTable)
           (Intercept) conditiontreat
test1_rep1           1              0
test1_rep2           1              0
test2_rep1           1              1
test2_rep2           1              1
attr(,"assign")
[1] 0 1
attr(,"contrasts")
attr(,"contrasts")$condition
[1] "contr.treatment"
# 估计离散度
y <- estimateDisp(y, design_full)
# 拟合模型
fit <- glmFit(y, design_full)
# 统计检验
lrt <- glmLRT(fit)
# 提取差异结果
tt <- topTags(lrt, n=nrow(y), sort="none")[[1]]
# 查看差异结果
head(tt,3)
                            logFC   logCPM         LR      PValue        FDR
ENSMUSG00000000001.5  -0.03700919 7.126297 0.03409904 0.853496315 0.93655267
ENSMUSG00000000028.16 -0.33464260 5.895729 2.18261856 0.139576806 0.36365238
ENSMUSG00000000031.17  2.67117263 3.563905 9.90557402 0.001647789 0.01718821

# 输出保存差异结果
tt$gene_id <- rownames(tt)
write.table(tt,file = 'edgeR_results.csv',sep = ',',row.names = FALSE)

3、limma-voom 差异分析代码

limma 包最早开发用来分析微整列芯片的表达数据（microarray），后来又开发了 voom 功能来针对 RNA-seq 数据进行差异分析。
因为 lima-voom 不使用存储在y$offset中的偏移矩阵，所以建议使用从丰度生成的缩放计数，scaledTPM或lengthScaledTPM。

# 加载R包
library(limma)
# 重新导入数据
txi <- tximport(file_dir, type = "salmon", tx2gene = tx2gene, countsFromAbundance = "lengthScaledTPM")
reading in files with read_tsv
1 2 3 4
summarizing abundance
summarizing counts
summarizing length

# 构建 DGEList object
y <- DGEList(txi$counts)
# 自动过滤低表达基因
keep <- filterByExpr(y)
y <- y[keep, ]
# 计算标准化因子
y <- calcNormFactors(y)

# 设计实验
sampleTable <- data.frame(condition = factor(rep(c("control", "treat"), each = 2)))
rownames(sampleTable) <- colnames(txi$counts)

design <- model.matrix(~condition, data = sampleTable)
v <- voom(y, design)
# 差异分析，先拟合在统计检验
fit <- lmFit(v,design)
fit2 <- eBayes(fit)
lim_res = topTable(fit2, coef=2, n=Inf)
# 查看差异结果
head(lim_res,3)
                          logFC  AveExpr         t      P.Value  adj.P.Val        B
ENSMUSG00000031548.8  -2.488670 5.957044 -14.20956 2.321833e-06 0.01446637 5.525443
ENSMUSG00000000037.18  2.028352 6.782320  12.33597 5.956494e-06 0.01446637 4.739079
ENSMUSG00000010796.12  3.766355 4.153263  13.39382 3.446536e-06 0.01446637 4.622380

# 输出保存差异结果
lim_res$gene_id <- rownames(lim_res)
write.table(lim_res,file = 'limma_results.csv',sep = ',',row.names = FALSE)