• 主页
  • 课程

    关于课程

    • 课程归档
    • 成为一名讲师
    • 讲师信息
    同等学历教学

    同等学历教学

    免费
    阅读更多
  • 特色
    • 展示
    • 关于我们
    • 问答
  • 事件
  • 个性化
  • 博客
  • 联系
  • 站点资源
    有任何问题吗?
    (00) 123 456 789
    weinfoadmin@weinformatics.cn
    注册登录
    恒诺新知
    • 主页
    • 课程

      关于课程

      • 课程归档
      • 成为一名讲师
      • 讲师信息
      同等学历教学

      同等学历教学

      免费
      阅读更多
    • 特色
      • 展示
      • 关于我们
      • 问答
    • 事件
    • 个性化
    • 博客
    • 联系
    • 站点资源

      未分类

      • 首页
      • 博客
      • 未分类
      • RNA-seq:Salmon 定量结果差异分析

      RNA-seq:Salmon 定量结果差异分析

      • 发布者 weinfoadmin
      • 分类 未分类, 老俊俊的生信笔记
      • 日期 2021年9月10日
      • 评论 0评论

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

      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)

      Part3差异结果比较

      用 DESeq2、edgeR 和 limma 三种差异软件鉴定到的显著性结果:

      软件 DESeq2 edgeR limma
      All significant gene numbers 1918 1627 1503
      Significant upregulated genes 842 641 625
      Significant downregulated genes 1076 986 878

      筛选标准:log2FC > 1 或 log2FC < -1 且 p < 0.05

      • 看起来鉴定到的差异基因的数量 DESeq2 是最多的,其次是 edgeR,limma 在这里面鉴定到的是最少的。
      • 差异基因总数和显著上调的基因数量中,edgeR 和 limma 鉴定到的比较相近,显著性下调基因的数量中,DESeq2 和 edgeR 比较相近。

      1、差异基因韦恩图
      在线用 vennDiagram 绘制韦恩图点击:shiny VennDiagram

      • 从 overlap 的基因数量来看,edgeR 和 limma overlap 的更多一些,其次是 DESeq2 和 edgeR,limma 与 DESeq2 overlap 到的基因是这里面最少的。

      2、差异基因火山图
      在线绘制火山图点击:画个 CNS 级别火山图!

      • 三个火山图里面都标注了一样的两个基因,可以看到在三个图里分布的趋势是一致的,log2FoldChange 的话,感觉都差不多,p 值的话,DESeq2 和 edgeR 比较相近,而整体 limma 结果的 p 值都比前两者要大很多,可能因为统计检验的方法不一样吧。

      参考资料

      [1]

      Swimming downstream: statistical analysis of differential transcript usage following Salmon quantification: http://www.bioconductor.org/packages/devel/workflows/vignettes/rnaseqDTU/inst/doc/rnaseqDTU.html#dge-analysis-with-deseq2

      [2]

      https://f1000research.com/articles/7-952/v3: https://f1000research.com/articles/7-952/v3

      [3]

      tximport: https://bioconductor.org/packages/release/bioc/vignettes/tximport/inst/doc/tximport.htm

      欢迎小伙伴留言评论!

      点击我留言!

      今天的分享就到这里了,敬请期待下一篇!

      最后欢迎大家分享转发,您的点赞是对我的鼓励和肯定!

      如果觉得对您帮助很大,打赏一下吧!


      请关注“恒诺新知”微信公众号,感谢“R语言“,”数据那些事儿“,”老俊俊的生信笔记“,”冷🈚️思“,“珞珈R”,“生信星球”的支持!

      • 分享:
      作者头像
      weinfoadmin

      上一篇文章

      RNA-seq:Salmon 快速定量
      2021年9月10日

      下一篇文章

      用R提取代表转录本
      2021年9月10日

      你可能也喜欢

      2-1675088548
      lncRNA和miRNA生信分析系列讲座免费视频课和课件资源包,干货满满
      30 1月, 2023
      9-1675131201
      如何快速批量修改 Git 提交记录中的用户信息
      26 1月, 2023
      8-1678501786
      肿瘤细胞通过改变CD8+ T细胞中的丙酮酸利用和琥珀酸信号来调控抗肿瘤免疫应答。
      7 12月, 2022

      留言 取消回复

      要发表评论,您必须先登录。

      搜索

      分类

      • R语言
      • TCGA数据挖掘
      • 单细胞RNA-seq测序
      • 在线会议直播预告与回放
      • 数据分析那些事儿分类
      • 未分类
      • 生信星球
      • 老俊俊的生信笔记

      投稿培训

      免费

      alphafold2培训

      免费

      群晖配置培训

      免费

      最新博文

      Nature | 单细胞技术揭示衰老细胞与肌肉再生
      301月2023
      lncRNA和miRNA生信分析系列讲座免费视频课和课件资源包,干货满满
      301月2023
      如何快速批量修改 Git 提交记录中的用户信息
      261月2023
      logo-eduma-the-best-lms-wordpress-theme

      (00) 123 456 789

      weinfoadmin@weinformatics.cn

      恒诺新知

      • 关于我们
      • 博客
      • 联系
      • 成为一名讲师

      链接

      • 课程
      • 事件
      • 展示
      • 问答

      支持

      • 文档
      • 论坛
      • 语言包
      • 发行状态

      推荐

      • iHub汉语代码托管
      • iLAB耗材管理
      • WooCommerce
      • 丁香园论坛

      weinformatics 即 恒诺新知。ICP备案号:粤ICP备19129767号

      • 关于我们
      • 博客
      • 联系
      • 成为一名讲师

      要成为一名讲师吗?

      加入数以千计的演讲者获得100%课时费!

      现在开始

      用你的站点账户登录

      忘记密码?

      还不是会员? 现在注册

      注册新帐户

      已经拥有注册账户? 现在登录

      close
      会员购买 你还没有登录,请先登录
      • ¥99 VIP-1个月
      • ¥199 VIP-半年
      • ¥299 VIP-1年
      在线支付 激活码

      立即支付
      支付宝
      微信支付
      请使用 支付宝 或 微信 扫码支付
      登录
      注册|忘记密码?