可变剪切背景知识可变剪接(Alternative Splicing,AS)是指从一个mRNA前体中通过不同的剪接方式,对外显子和内含子进行组合,产生不同的mRNA剪接异构体的过程。高等真核生物中的可变剪接极大地拓展了基因功能的多样性,是调节基因表达和产生蛋白质组多样性的重要机制。 可变剪接受到具有特殊结构域的顺势调控元件(RNA motif)和识别这些motif的RNA结合蛋白(RNA binding protein,RBP)调控 。RNA-seq通常是二代转录组,可以通过高深度的测序数据组装构建转录本序列,预测外显子与内含子的结构并识别出可变剪接模式,假阳性不小。三代全长转录组利用其读长更长的优势,可以直接读取转录本的全长序列。相比二代测序,无需打断、无需组装,直接获得全长转录本的结构信息,可以更准确的分析生物体内存在可变剪接事件。 Multivariate Analysis of Transcript Splicing (MATS)软件列出5种可变剪切,分别是skipped exon (SE)外显子跳跃,alternative 5′ splice site (A5SS)第一个外显子可变剪切,alternative 3′ splice site (A3SS)最后一个外显子可变剪切,mutually exclusive exons (MXE)外显子选择性跳跃和 retained intron (RI)内含子滞留,示意图在官网,如下: SplicSeq软件可以拿到7种可变剪切形式信息:
如果你看TCGA的ISOexpresso数据库:http://wiki./ISOexpresso/ 上面有8种: 可变剪切在癌症领域研究例子可变剪切在癌症领域研究不多,但重要性毋庸置疑,主要是因为技术手段限制,值得介绍的经典文章,比如邵志敏教授的题为“PHF5A Epigenetically Inhibits Apoptosis to Promote Breast Cancer Progression.”利用CRISPR-Cas9文库技术对RNA结合蛋白进行系统性功能筛选,发现了乳腺癌生存依赖的的剪切因子PHF5A . TCGA、METABRIC、KMplot以及本中心的研究队列显示,PHF5A在乳腺癌组织中普遍存在高表达,并且其高表达指示着较差的预后。生物学功能上,PHF5A缺失会导致癌细胞增殖、迁移和成瘤能力显著受损。该文章的测序数据在 SRP139147 ,然后走RNA-seq分析流程,比较cancer cells (CA1a or DCIS) 和noncancer cells (MCF10AT or MCF10A)差异分析,虽然作者文章主要是关心那1,542 manually curated RBPs ( RNA-binding proteins) ,不过最后落脚点仍然是可变剪切: 文章写作落脚点是PHF5A knockdown appeared to switch full-length FASTK (FASTK-L) to an intron 5-retained variant (herein termed FASTK short, FASTK-S) in CA1a cells 是通过差异分析筛选出来的。 这个时候,有两个非常经典的软件leafcutter和rMATS,我都在生信技能树写过教程,两年前过去了,现在又需要重新使用,是时候更新一下软件和用法了。 首先看软件各自的官网安装leafcutter 软件相关地址:
rMATS 软件相关地址:
全部使用conda一键式安装吧: conda create -n isoform 因为conda可以自动帮你解决全部的依赖: 比如 r-leafcutter就依赖于一大波的R包,而且leafcutter的GitHub源代码里面有一些脚本和测试数据,所以还是要下载看看 mkdir -p ~/biosoft 再次吐槽一下,这个软件包里面涵盖了,perl,python,r,sh代码,开发流程非常的不统一! 然后看看软件各自的安装和使用首先让我们先回顾一下leafcutter 软件的4个标准步骤:第一个步骤是shell脚本bam2junc.sh把bam文件转为junc文件,可以构建好bam_path.txt文件,存储全部的bam文件路径然后批量处理,第一个步骤全部的bam文件输出的junc文件路径保存在 all_juncfiles.txt 。 第二个步骤是使用python脚本 leafcutter_cluster.py做 Intron clustering,就是把第一个步骤全部的bam文件输出的junc文件合并得到 ls ~/data/public/star/*bam > bam_path.txt 第3个步骤是使用R脚本 leafcutter_ds.R :对整个项目的样本进行分组,然后对它们的 chr start end strand gene_name 自己写脚本制作,比如: zcat ~/biosoft/leafcutter/leafcutter/data/gencode19_exons.txt.gz|head 然后,分组文件 group_info.txt 也需要自己制作,是两列的格式, 样本名和分组,举例如下: SRR2016934 control 运行第三步代码如下: exonbed=/home/yb77613/biosoft/leafcutter/leafcutter/data/gencode_v32_hg38_exons.txt 第4个步骤是使用R脚本 ds_plots.R 对差异分析结果进行可视化。 exonbed=/home/yb77613/biosoft/leafcutter/leafcutter/data/gencode_v32_hg38_exons.txt 所有的统计学显著的可变剪切形式都会可视化在一张图里面,也可以自己别出心裁的使用其它可视化方法。我们随机展示其中一个结果(MFF基因): chr2:clu_4932_NA Success 4.6012973874349 7 0.238436066745941 0.634806885002018 MFF 看外显子表达量是: 软件就是根据分组,做一下统计检验看看是否有差异。 接着让我们看看rMATS v4.0.2 (turbo)的用法它只有一个步骤,直接根据分组对bam文件进行可变剪切的差异分析,通过统计模型对不同样本(有生物学重复的)进行可变剪切事件的表达定量,然后以likelihood-ratio test计算P value来表示两组样品在IncLevel(Inclusion Level)水平上的差异(从公式上来看,IncLevel跟PSI的定义也是类似的),lncLevel并利用Benjamini Hochberg算法对p value进行校正得FDR值。 gtf="/home/yb77613/reference/gtf/gencode/gencode.v32.annotation.gtf" 主要运行过程中,是一些python模块,或者系统库文件缺失会报错,比如我遇到的: Running the statistical part. 搜索了一些解决方案,最后采纳: echo $LD_LIBRARY_PATH 我们运行的rMATS 程序首先会解析gtf文件: There are 60609 distinct gene ID in the gtf file 然后读取全部的bam文件。 Done processing each gene from dictionary to compile AS events 但是,很明显,仍然是报错的, 因为图便利,采用了conda安装,代码是: 可以看到,conda虽然说能解决初学者99%的软件安装问题,但也不是万能的。 比较两个分析的结果rMATS运行失败,懒得去解决它的bugs了,以后再说。 和salmon加DRIMSeq流程比较前面我们介绍过,不需要走bam这个文件格式做中间产物,在Swimming downstream: statistical analysis of differential transcript usage following Salmon quantification 流程里面也可以进行转录本差异分析
可变剪切的结果,其实就是转录本差异分析,所以也可以互相比较的。比如前面的MFF基因,搜索: gene_id "ENSG00000168958.19"; transcript_id "ENST00000456345.1"; gene_type "protein_coding"; gene_name "MFF"; 对应的转录本ID,发现也是差异的: 这个基因的3个差异转录本,可视化如下: 可以看到,寻找差异转录本,或者差异外显子的分析思路,是基于已有的注释,在人类这个物种的研究当然是非常顺畅的。 既然不同的算法,拿到的结果大同小异,我们可以放心的走下一步啦! 生信技能树可变剪切相关教程节选因为做目录确实很浪费时间,差不多就下面这些,大家先学习吧: 可以选择参加11月23号周六下午3-6点,上海张江高科的义诊,https://mp.weixin.qq.com/s/xIOzD4XT8_G8LCZP0kZxgw 或者11月24号周日上午9-12点的拍卖行,两种形式都可以帮助你解决生物信息学问题,https://mp.weixin.qq.com/s/wNi-Ips8WKffAAZAmpo3NA 同时我们还会提供一些招聘见面会,比如 https://mp.weixin.qq.com/s/0i51Fv5pdu4f8W06RGIdGQ 号外:生信技能树全国巡讲11月在福州和上海,点击了解报名哈:(福州、上海见!)全国巡讲第19-20站(生信入门课加量不加价) 福州已经停止报名,上海还可以报名,不过位置也不多了,然后下一站是12月底的长沙,欢迎咨询! |
|