• 主页
  • 课程

    关于课程

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

    同等学历教学

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

      关于课程

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

      同等学历教学

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

      未分类

      • 首页
      • 博客
      • 未分类
      • rMATS可视化-rmats2sashimiplot

      rMATS可视化-rmats2sashimiplot

      • 发布者 weinfoeditor
      • 分类 未分类
      • 日期 2019年8月21日
      • 评论 0评论

      之前的对rMATS软件做笔记的时候,提到过rmats2sashimiplot这款专门用来对rMATS分析结果做可视化的软件,最近尝试用了下,操作比较简洁,但是还是留了点问题没解决

      rmats2sashimiplot,网址:rmats2sashimiplot

      下载

      git clone https://github.com/Xinglab/rmats2sashimiplot
      
      or
      
      https://pypi.org/project/rmats2sashimiplot/2.0.2/#files
      

      安装

      python setup.py install
      

      ubuntu 16.04 报错:error while loading shared libraries: libgsl.so.0: cannot open shared object file: No such file or directory

      这是因为ubuntu 16.04不支持libgsl0ldbl,而是替换为libgsl2,所以先安装libgsl2,然后在将其软链接/复制到/usr/lib目录下

      sudo apt-get install libgsl-dev
      dpkg -L libgsl2
      sudo cp /usr/lib/x86_64-linux-gnu/libgsl.so.19 /usr/lib/libgsl.so.0
      

      解决方法参照自:https://jingyan.baidu.com/article/63acb44a13cc7d61fdc17e6d.html


      如果直接使用rMATS从fastq输入文件到结果输出文件,那么rMATS是不保留中间的bam文件的

      后来看了下rmats.py中的代码,发现只要将shutil.rmtree(args.tmp)注释掉就会保留STAR比对的bam结果文件了,但是其是放在服务器上/tmp目录中,也不太方便

      而rmats2sashimiplot可视化则需要bam文件作为输入,所以需要我们先用STAR比对得到bam文件再用rMATS做差异可变剪切分析,如下所示:

      1. 用STAR对各个样本做比对生成bam文件,比对参数参考rMATS软件调用STAR时所用参数,其实就是比默认比对参数多了:

        --chimSegmentMin 2 --outFilterMismatchNmax 3 --outSAMstrandField intronMotif --alignEndsType EndToEnd --alignSJDBoverhangMin 6 --alignIntronMax 299999

        STAR --readFilesIn A1_1.clean.fq A1_2.clean.fq --chimSegmentMin 2 --outFilterMismatchNmax 3 --alignEndsType EndToEnd --runThreadN 12 --outSAMstrandField intronMotif --outSAMtype BAM SortedByCoordinate --alignSJDBoverhangMin 6 --alignIntronMax 299999 --genomeDir ~/index/STAR/Ensembl/human/GRCh38.92/ --sjdbGTFfile ~/annotation/Ensembl/human/Homo_sapiens.GRCh38.92.chr.gtf --outFileNamePrefix A1 
        
      2. rMATS做差异可变剪切分析(例如A.txt包含了A组bam的文件名)

        python rmats.py --b1 A.txt --b2 B.txt --gtf Homo_sapiens.GRCh38.92.chr.gtf --od ./test -t paired --readLength 150 --cstat 0.0001 --nthread 10
        

      rmats2sashimiplot作图

      rmats2sashimiplot --b1 A1.bam,A2.bam,A3.bam --b2 B1.bam,B2.bam,B3.bam -t SE -e ./SE.MATS.JC_top20.txt --l1 A --l2 B --exon_s 1 --intron_s 1 -o SE_plot
      

      rmats2sashimiplot的注意事项:

      • 可以预先用samtools对bam文件建索引,不然rmats2sashimiplot也会先建索引(比较慢)
      • rMATS的结果文件中,如SE.MATS.JC.txt中染色体都是默认加上chr的,而bam文件中的染色体根据基因组注释文件来源不同,有些的没chr,而是直接以数字表示,这时需要预先将两者保持一致
      • rmats2sashimiplot默认是出PDF图片的,如果需要出png,则需要自己修改下源码(由于我也不懂python作图,所以用了最无脑的方式就是将所有跟出图相关的脚本中pdf都替换为png。。。)

      rMATS可视化的图片如下:

      rmats2sashimiplot

      但是有个问题我还是没想清楚,按照rmats2sashimiplot作者的意思,这图中的线上数字代表着 junction-spanning reads,那么这跟rMATS的分析结果中的IC_SAMPLE_1和SC_SAMPLE_1一致吗?

      比如这张图中第一个样本A-1的Exon Inclusion Isoform的reads数分别为176和159,exon skipping isoform的reads为2

      而rMATS分析结果中IC_SAMPLE_1的A-1样本数值为259,SC_SAMPLE_1为2;

      我的理解图中Exon Inclusion Isoform上的两个数值加和应该等于259(有些图片确实相等),但事实好多图片都与我所想不一致,有时相差的还很奇怪,想不明白。。。


      之前在用rMATS时,对其结果的文件的表头并不是理解的很清楚,最近查了下资料(主要是rMATS的原理文献,反正也没完全看懂),再整理下:

      IncLevel1可被认作是exon inclusion level(ψ),是Exon Inclusion Isoform在总(Exon Inclusion Isoform + Exon Skipping Isoform)所占比例

      For a skipped exon, the exon inclusion level (denoted as percent spliced in or ψ) is calculated by the number of reads uniquely mapped to the exon inclusion isoform or the exon skipping isoform

      因此计算公式如下:

      ψ = (I/lI)/(S/lS)

      • I 是指mapping到exon inclusion isoform上的reads数,对应结果表格中的IC_SAMPLE_1
      • S 是指mapping到exon skipping isoform上的reads数,对应结果表格中的SC_SAMPLE_1
      • lI 是指effective length of the exon inclusion isoform,对应结果表格中的IncFormLen,主要用于对I标准化,用于计算ψ
      • lS 是指effective length of the exon skipping isoform,对应结果表格中的SkipFormLen,主要用于对S标准化,用于计算ψ

      对于effective length的计算公式则是:

      For a segment whose length is l, and the length of the read is r, the effective length of the segment is defined by the number of unique read intervals in this region, which is l − r + 1

      这里的计算可以参照https://www.biostars.org/p/256949/,我还是没搞清楚。。。

      IncLevelDifference则是指两组样本IncLevel的差异,如果一组内多个样本,那么则是各自组的均值之间差值

      rMATS采用的是Likelihood-ratio test来比较两组样本可变剪切是否有差异

      If the maximum-likelihood estimations (MLEs) of ψi1,ψi2 have a difference smaller than or equal to the user-defined cutoff (i.e., |ψi1 −ψi2|≤ c), we set the P value to be 1

      所以跟我们参数--cstat设置有关(默认是0.0001),基于marginal distribution(在零假设前提下,近似χ2 distribution),备择假设是|ψi1 −ψi2|> c,因此粗略的理解就是ψi1和ψi2差异比c越大P值就越小

      虽然上面作图可视化的问题还是没搞清楚,但是也记录下了,之后如果搞清楚了再回来修改。。

      本文出自于http://www.bioinfo-scrounger.com转载请注明出处

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

      • 分享:
      作者头像
      weinfoeditor

      上一篇文章

      R语言-Survival analysis(生存分析)
      2019年8月21日

      下一篇文章

      R语言 PCA分析
      2019年8月21日

      你可能也喜欢

      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年
      在线支付 激活码

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