• 主页
  • 课程

    关于课程

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

    同等学历教学

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

      关于课程

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

      同等学历教学

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

      未分类

      • 首页
      • 博客
      • 未分类
      • 鉴定差异翻译效率基因之 Riborex

      鉴定差异翻译效率基因之 Riborex

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

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

      1、介绍

      今天介绍一个 Ribo-seq 中用于快速鉴定差异翻译效率基因的软件:Riborex。文章在 2017 年发表在 Bioinformatics 上:

      标题:

      Riborex: fast and flexible identification of differential translation from Ribo-seq data[1]

      Riborex 基于在已有的差异分析软件 DESeq2、edgeR、Voom 的框架下,根据负二项模型(negative binomial distribution)和广义线性模型(GLM)来寻找翻译效率差异的基因。作者对不同的翻译效率计算软件和 Riborex 里 3 种不同的方法鉴定的基因进行 overlap,都有大部分的重叠:

      此外对不同数据集的 Ribo-seq 数据使用不同软件分析,相对于之前的几款软件,所耗费的时间也是最快的:

      不同的数据集之间,Riborex 的准确率相对于其它软件也有着良好的效果,此外 Xtail 软件也表现不错:

      总之,Riborex 是作为一个开源的 R 包实现的,它的准确性与现有的最精确的方法不相上下,但是速度要快得多。

      2、安装

      包的信息及软件存放在 githup 网址上:

      smithlabcode/riborex:[2] https://github.com/smithlabcode/riborex

      作者在 3 个月前还更新了 README 文件,此外在 vignetts 文件夹里有个 manual 使用帮助文档,可以参考使用。

      该软件依赖 DESeq2、edgeR 和 fdrtool 3 个包,提前装好,作者强烈建议使用 conda 安装:

      $ conda install -c bioconda r-riborex

      R 里安装:

      # 1、安装包
      BiocManager::install("DESeq2")
      BiocManager::install("edgeR")
      install.packages('fdrtool')
      # 从githup安装
      library(devtools)
      options(unzip='internal')
      devtools::install_github('smithlabcode/riborex')

      # 2、加载R包
      library(riborex)

      3、使用

      输入文件需要 RNA-seq 和 Ribo-seq 的定量 count 表达矩阵,可以使用 HT-seq 或 featureCount 等软件得到,接下来使用软件来鉴定差异翻译效率的基因:

      # 加载测试数据
      data(riborexdata)
      # 获取rna-seq和ribo-seq count
      RNACntTable <- rna
      RiboCntTable <- ribo
      # 查看数据
      head(RNACntTable,3)
                         BN_336 BN_337 BN_338 BN_339
      ENSRNOG00000000017      7     11      4      4
      ENSRNOG00000000024   2467   2478   3258   2316
      ENSRNOG00000000033    206    282    330    244
      head(RiboCntTable,3)
                         BN_341 BN_342 BN_343 BN_344
      ENSRNOG00000000017     15      5     10      2
      ENSRNOG00000000024   5206   5921   2864   1985
      ENSRNOG00000000033     30     30     23     13
      # 分组
      rnaCond <- c("control", "control", "treated", "treated")
      riboCond <- c("control", "control", "treated", "treated")
      # 差异分析,使用DESeq2
      res.deseq2 <- riborex(RNACntTable, RiboCntTable, rnaCond, riboCond, engine = "DESeq2")
      DESeq2 mode selected
      combining design matrix
      applying DESeq2 to modified design matrix
      estimating size factors
      estimating dispersions
      gene-wise dispersion estimates
      mean-dispersion relationship
      final dispersion estimates
      fitting model and testing
      Warning message:
      In DESeqDataSet(se, design = design, ignoreRank) :
        some variables in design formula are characters, converting to factors

      # 查看结果
      head(res.deseq2,3)
      log2 fold change (MLE): EXTRA1 treated vs control
      Wald test p-value: EXTRA1 treated vs control
      DataFrame with 3 rows and 6 columns
                           baseMean log2FoldChange     lfcSE      stat    pvalue      padj
                          <numeric>      <numeric> <numeric> <numeric> <numeric> <numeric>
      ENSRNOG00000000017    7.69493       1.533227  1.566415  0.978813  0.327672  0.736753
      ENSRNOG00000000024 3544.23807      -0.102151  0.202054 -0.505564  0.613163  0.901926
      ENSRNOG00000000033  111.43945       0.265252  0.503397  0.526924  0.598246  0.896407

      检查 p 值的分布:

      • 良好的 p 值分布图应该是在 0 侧有个最高的峰,往后面则是均匀分布,具体可参考:什么,你算出的 P-value 看上去像齐天大圣变的庙?:
      # 检查p值分布
      hist(res.deseq2$pvalue, main = 'DESeq2 unadjusted p-values', xlab='Unadjusted p-values',
           col = '#F4C7AB')

      统计差异结果:

      # 统计差异结果
      summary(res.deseq2)
      out of 13916 with nonzero total read count
      adjusted p-value < 0.1
      LFC > 0 (up)       : 1107, 8%
      LFC < 0 (down)     : 1217, 8.7%
      outliers [1]       : 0, 0%
      low counts [2]     : 540, 3.9%
      (mean count < 6)
      [1] see 'cooksCutoff' argument of ?results
      [2] see 'independentFiltering' argument of ?results

      保存差异结果:

      # 保存结果
      write.table(res.deseq2, "riborex_res_deseq2.txt", quote=FALSE,sep = 't')

      使用 edgeR 来差异分析:

      # 使用edgeR来差异分析
      res.edgeR <- riborex(RNACntTable, RiboCntTable, rnaCond, riboCond, "edgeR")
      edgeR mode selected
      combining design matrix
      applying edgeR to modified design matrix

      # 查看结果
      head(res.edgeR$table,3)
                               logFC    logCPM         LR    PValue       FDR
      ENSRNOG00000000017  1.36290149 -2.456968 1.80649827 0.1789289 0.4627345
      ENSRNOG00000000024 -0.30127172  6.212404 2.13670654 0.1438103 0.4037250
      ENSRNOG00000000033  0.07178854  1.313235 0.02631769 0.8711269 0.9636097

      # 使用edgeRD来对RNA-seq和Ribo-seq分别评估离散度
      res.edgeRD <- riborex(RNACntTable, RiboCntTable, rnaCond, riboCond, "edgeRD")

      使用 limma 包的 Voom 来进行差异分析:

      # 使用limma包的Voom来进行差异分析:
      res.voom <- riborex(RNACntTable, RiboCntTable, rnaCond, riboCond, "Voom")
      Voom mode selected
      combining design matrix
      applying Voom to modified design matrix

      # 查看差异结果
      head(res.voom,3)
                               logFC    AveExpr          t   P.Value adj.P.Val         B
      ENSRNOG00000000017  1.39674189 -2.8437969  1.2847048 0.2268559 0.5243467 -4.928058
      ENSRNOG00000000024 -0.30047073  6.0075571 -1.8737990 0.0893910 0.2991861 -5.682138
      ENSRNOG00000000033  0.07661457  0.7070877  0.1687028 0.8692764 0.9540355 -6.320171

      4、矫正 p 值

      对于我们的 p 值分布如果不是那么理想,很有可能鉴定出来的差异基因数量会很少,作者还提供了矫正 p 值的解决办法,让 p 值分布尽可能 “correct” :

      # 操作基本一致
      # 载入数据
      RNACntTable.corrected <- rna.null
      RiboCntTable.corrected <- ribo.null
      # 分组
      rnaCond.corrected <- c(rep('T0', 3), rep('T24',3))
      riboCond.corrected <- rnaCond.corrected
      # 差异分析
      res.deseq2.corrected <- riborex(RNACntTable.corrected, RiboCntTable.corrected,
                                      rnaCond.corrected, riboCond.corrected)
      DESeq2 mode selected
      combining design matrix
      applying DESeq2 to modified design matrix
      estimating size factors
      estimating dispersions
      gene-wise dispersion estimates
      ...

      # 检查p值分布
      hist(res.deseq2.corrected$pvalue, main = 'DESeq2 unadjusted p-values', xlab='Unadjusted p-values',
           col = '#F4C7AB')

      可以看到这个 p 值分布是向右倾斜逐渐升高的,p < 0.05 基因就会很少,接下来用 fdrtool 进行矫正:

      # 矫正
      results.corrected <- correctNullDistribution(res.deseq2.corrected)
      Step 1... determine cutoff point
      Step 2... estimate parameters of null distribution and eta0
      Step 3... compute p-values and estimate empirical PDF/CDF
      Step 4... compute q-values and local fdr
      # 检查矫正后的p值分布
      hist(results.corrected$pvalue, main = 'DESeq2 unadjusted p-values after correction',
           xlab='Corrected unadjusted p-values',col = '#F4C7AB')

      5、多因子实验设计的差异分析

      此外,此软件还支持多个实验条件设计的差异翻译效率的分析:

      # 多因子实验设计
      rna <- RNACntTable[,c(1,2,3,4,1,2,3,4)]
      ribo <- RiboCntTable[,c(1,2,3,4,1,2,3,4)]
      # 分组
      rnaCond <- data.frame(factor1=(c("control1", "control1", "treated1", "treated1",
                                       "control1", "control1", "control1", "control1")),
                            factor2=(c("control2", "control2", "control2", "control2",
                                       "control2", "control2", "treated2", "treated2")))
      riboCond <- data.frame(factor1=(c("control1", "control1", "treated1", "treated1",
                                        "control1", "control1", "control1", "control1")),
                             factor2=(c("control2", "control2", "control2", "control2",
                                        "control2", "control2", "treated2", "treated2")))

      进行差异分析需要提供指定的比较分组:

      # 进行差异分析需要提供指定的比较分组
      contrast = c("factor2", "control2", "treated2")
      # 差异分析
      res.deseq2 <- riborex(rna, ribo, rnaCond, riboCond, "DESeq2", contrast = contrast)
      # 查看统计结果
      summary(res.deseq2)
      out of 13916 with nonzero total read count
      adjusted p-value < 0.1
      LFC > 0 (up)       : 1887, 14%
      LFC < 0 (down)     : 1987, 14%
      outliers [1]       : 0, 0%
      low counts [2]     : 270, 1.9%
      (mean count < 3)
      [1] see 'cooksCutoff' argument of ?results
      [2] see 'independentFiltering' argument of ?results

      # edgeR和edgeRD类似
      res.edgeR <- riborex(rna, ribo, rnaCond, riboCond, "edgeR", contrast = contrast)
      res.edgeRD <- riborex(rna, ribo, rnaCond, riboCond, "edgeRD", contrast = contrast)

      6、绘制火山图

      画个火山图对上下调基因可视化一下:

      # 转换为data.frame
      dat <- as.data.frame(res.deseq2)
      # 定义上下调基因
      dat$type <- ifelse(abs(dat$log2FoldChange)>1 & dat$pvalue < 0.05,
                         ifelse(dat$log2FoldChange >1 & dat$pvalue < 0.05,'UP','DOWN'),'nonSig')
      # 查看上下调基因数量
      table(dat$type)
      DOWN nonSig     UP
         983  12011    922

      # 绘图
      library(ggplot2)
      library(ggprism)

      ggplot(dat,aes(x = log2FoldChange,y = -log10(pvalue))) +
        geom_point(aes(color = type,shape = type),alpha = .5,size = 3) +
        theme_prism() + ylim(0,80) + xlim(-6,6) +
        theme(legend.position = c(0.9,0.9)) +
        scale_color_manual(values = c('UP'= '#F54748','nonSig'='#F8A488','DOWN'='#5AA897')) +
        geom_hline(yintercept = -log10(0.05),linetype = 'dashed',size = 1) +
        geom_vline(xintercept = c(-1,1),linetype = 'dashed',size = 1)

      7、实验没有 replicate 能分析吗?

      对于实验没有 replicate 是不能做统计学检验分析的,但是 edgeR 是可以对 RNA-seq 在没有生物学重复的情况下进行差异分析,对于差异翻译效率分析的话,我也不知道 Riborex 是否能对于没有重复的条件下进行差异分析。

      于是昨天在 githup 上问了作者这个问题,今天中午便很快得到了答复!

      很不幸的是,不可以哟。

      参考资料

      [1]

      Riborex: fast and flexible identification of differential translation from Ribo-seq data: https://academic.oup.com/bioinformatics/article/33/11/1735/2964727#118738289

      [2]

      smithlabcode/riborex:: https://github.com/smithlabcode/riborex

      欢迎小伙伴留言评论!

      点击我留言!

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

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

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

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

      • 分享:
      作者头像
      weinfoadmin

      上一篇文章

      鉴定差异翻译效率基因之 deltaTE 上篇
      2021年9月10日

      下一篇文章

      purrr 包之 list 处理系列函数
      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年
      在线支付 激活码

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