查看原文
其他

用LeafCutter探索转录组数据的可变剪切

jimmy 生信技能树 2022-06-07

该软件早在2016年就公布了,发表在biorxiv预印本上面,但直到2017年的双11,才发表在NG上面,文章是 : Annotation-free quantification of RNA splicing using LeafCutter 最大的特点应该是不需要参考基因组的基因注释信息了吧,就是gtf/gff文件可以省略,当然,比对还是需要的。它还有另外一个非常重要的功能,splicing quantitative trait loci (sQTLs) 但是跟我目前关系不大, 就不介绍了。

本教程,首发于 生信菜鸟团博客:http://www.bio-info-trainee.com/2949.html

背景介绍

目前主流的探究转录组数据的可变剪切的算法要么是基于estimate isoform ratios 或者 exon inclusion levels ,但是挑战还是蛮多的,可变剪切本跟正常转录本重合的比例很大,技术误差也是有的,依赖于基因现有的注释信息,既不准确,也不完全。所以作者开发了LeafCutter工具。

LeafCutter workflow.

  • First, short reads are mapped to the genome. When SNP data are available, WASP should be used to filter allele-specific reads that map with a bias.

  • Next, LeafCutter extracts junction reads from.bam files, identifies alternatively excised intron clusters, and summarizes intron usage as counts or proportions.

  • Finally, LeafCutter identifies intron clusters with differentially excised introns between two user-defined groups by using a Dirichlet-multinomial model, or maps genetic variants associated with intron excision levels by using a linear model.

作者在Genotype-Tissue Expression (GTEx) Consortium数据集上面测试了,并且把结果跟 GENCODE v19, Ensembl, and UCSC 着3大主流的基因注释信息数据库比较。还在其它数据库里面验证了,数据下载地址是:dbGaP under accession phs000424.v6.p1 (GTEx), GEO under accession GSE41637 (RNA-seq data from mammalian organs), and ENA under accession PRJEB3366 (Geuvadis).

软件下载地址:

  • LeafCutter software, https://github.com/davidaknowles/leafcutter;

  • LeafViz visualizations, https://leafcutter.shinyapps.io/leafviz/;

  • rheumatoid arthritis summary statistics, http://plaza.umin.ac.jp/yokada/datasource/software.htm.

软件安装及使用

最简单的就是conda进行安装了:

  1. conda install -c davidaknowles r-leafcutter

如果安装失败,可能需要单独为它创建一个环境。

不过,它本身就是一个R包,所以在个人电脑里面的rstudio里面安装即可。

  1. if (!require("devtools")) install.packages("devtools", repos='http://cran.us.r-project.org')

  2. devtools::install_github("davidaknowles/leafcutter/leafcutter")

但是源代码里面有一些脚本和测试数据,所以还是要下载看看

  1. mkdir -p ~/biosoft

  2. cd ~/biosoft

  3. git clone https://github.com/davidaknowles/leafcutter

  4. cd leafcutter

  5. ## 需要修改里面的一个脚本 scripts/bam2junc.sh 把软件路径增添进去即可

里面又是perl又是python的,感觉他们团队开发环境不统一。

第一步:bam2junc

比对一般来说,优先选择STAR等支持跨越内含子的转录组比对工具得到bam文件,运行下面的脚本即可进行批量转换:

  1. cat bam_path.txt |while read id

  2. do

  3. file=$(basename $id )

  4. sample=${file%%.*}

  5.    echo Converting $id to $sample.junc

  6.    sh /public/biosoft/leafcutter/scripts/bam2junc.sh  $id $sample.junc

  7. done

得到的junc文件如下:

  1. chr7    134840725   134843893   .   1   -

  2. chr2    234355442   234355737   .   1   +

  3. chr4    37828435    37831585    .   13  +

  4. chr19    39101772    39101882    .   5   +

  5. chr11    109735445   109827551   .   19  +

  6. chr18    48458730    48465939    .   8   -

  7. chr12    82751048    82752457    .   12  -

  8. chr15    51018323    51018517    .   14  -

  9. chr1    247323115   247335149   .   2   +

  10. chr10    92920631    92982445    .   1   +

这个步骤有点耗时,所有的junc文件地址需要保存给下一步使用

第二步:Intron clustering

这个步骤,需要python2.7版本,这个是python的一个大坑,到现在版本仍然不统一。

  1. ls *.junc >test_juncfiles.txt

  2. python /public/biosoft/leafcutter/clustering/leafcutter_cluster.py -j test_juncfiles.txt -m 50 -o testYRIvsEU -l 500000

几分钟就运行完毕。

得到的比较重要的文件如下:

  1. 1.3M Jan  4 17:45 testYRIvsEU_perind.counts.gz

  2. 680K Jan  4 17:45 testYRIvsEU_perind_numers.counts.gz

  3. 5.0M Jan  4 17:45 testYRIvsEU_pooled

  4. 540K Jan  4 17:45 testYRIvsEU_refined

  5. 877 Jan  4 17:45 testYRIvsEU_sortedlibs

  6. 854 Jan  4 17:43 test_juncfiles.txt

值得注意的是 testYRIvsEU_perind_numers.counts.gz 文件,里面每一行都是一个内含子,每一列都是一个样本,写明了它们的表达值,这些数值就可以用来做可变剪切分析。

  1. #  zcat testYRIvsEU_perind_numers.counts.gz |tail

  2. chr8:145651155:145651305:clu_6538 21 14 19 8 9 0 13 33 0 0 4 0 5 8 12 0 12 34 15 0 0 10 11

  3. chr8:145651155:145651409:clu_6538 1021 611 186 190 294 284 681 89 222 57 257 363 694 807 523 44 469 812 926 71 80 260 214

  4. chr8:145652362:145653872:clu_6539 1265 694 132 74 302 71 178 34 44 12 63 122 230 218 472 6 146 1421 1084 16 14 83 46

  5. chr8:145652654:145653872:clu_6539 48 24 56 0 26 0 13 0 2 5 2 0 3 19 17 0 2 8 64 0 0 3 0

  6. chr8:145652674:145653872:clu_6539 18 26 0 0 0 7 2 0 5 0 0 0 1 6 11 0 3 34 37 0 0 9 6

  7. chr8:146017525:146017630:clu_6540 2 3 44 0 2 12 4 0 0 0 22 5 9 10 2 0 1 9 11 0 0 1 0

  8. chr8:146017525:146017751:clu_6540 1067 671 620 41 295 347 224 89 62 33 262 136 229 223 356 17 288 480 1842 9 35 70 23

  9. chr8:146076780:146078224:clu_6541 18 3 0 0 17 17 8 0 0 3 2 3 16 6 12 0 4 45 29 9 0 10 2

  10. chr8:146076780:146078378:clu_6541 22 17 0 0 0 3 1 0 0 0 3 2 15 7 2 0 7 62 55 0 0 4 0

  11. chr8:146076780:146078757:clu_6541 10 1 16 0 12 52 0 0 11 0 24 9 27 3 0 0 7 0 28 0 0 2 0

第三步:制作分组矩阵进行差异分析

避免暴露我真实的项目,这里就给作者的示例文件吧:

  1. RNA.NA18486_YRI.chr1.bam YRI

  2. RNA.NA18487_YRI.chr1.bam YRI

  3. RNA.NA18488_YRI.chr1.bam YRI

  4. RNA.NA18489_YRI.chr1.bam YRI

  5. RNA.NA18498_YRI.chr1.bam YRI

  6. RNA.NA06984_CEU.chr1.bam CEU

  7. RNA.NA06985_CEU.chr1.bam CEU

  8. RNA.NA06986_CEU.chr1.bam CEU

  9. RNA.NA06989_CEU.chr1.bam CEU

  10. RNA.NA06994_CEU.chr1.bam CEU

很简单的两列文件,说明每一个样本属于哪个组即可。

  1. /public/biosoft/leafcutter/scripts/leafcutter_ds.R --num_threads 4 \

  2. --exon_file=/public/biosoft/leafcutter/leafcutter/data/gencode19_exons.txt.gz \

  3. testYRIvsEU_perind_numers.counts.gz group_info.txt

这里的 group_info.txt 就是自己制作好的分组矩阵。值得提醒的是,上面的文件有且只能有2个分组,这样软件才知道怎么样去比较,如果自己的分组很多,可以考虑制作多个分组文件,运行多次。

当然,上面的脚本已经没有必要在linux服务器里面运行啦。

既然有了内含子的表达矩阵,又有了分组信息,差异分析根本就不会消耗多少计算资源,全部下载到自己的电脑里面去做吧。

自己打开文件 /public/biosoft/leafcutter/scripts/leafcutter_ds.R 就明白了整个流程。

也是几分钟就完成了全部结果。

  1. Running differential splicing analysis...

  2. Differential splicing summary:

  3.                                             statuses Freq

  4. 1 <2 introns used in >=min_samples_per_intron samples  425

  5. 2                          <=1 sample with coverage>0   62

  6. 3               <=1 sample with coverage>min_coverage  939

  7. 4                            Not enough valid samples 3047

  8. 5                                             Success 2068

  9. Saving results...

  10. Loading exons from /Users/jmzeng/biosoft/leafcutter/leafcutter/data/gencode19_exons.txt.gz

  11. All done, exiting

得到的文件里面,需要详细了解的是 leafcutter_ds_cluster_significance.txt 主要靠自己看readme啦。

第四步:可视化那些可变剪切

也是包装好的脚本。

  1. /Users/jmzeng/biosoft/leafcutter/scripts/ds_plots.R -e  /Users/jmzeng/biosoft/leafcutter/leafcutter/data/gencode19_exons.txt.gz testYRIvsEU_perind_numers.counts.gz   group_info.txt leafcutter_ds_cluster_significance.txt -f 0.05

所有的可变剪切形式都会可视化在一张图里面。

可视化函数还是有很多参数是可以自己调整的,值得学习


leafcutter得到剪切形式


上面的剪切形式略微有点复杂。


leafcutter得到剪切形式


这个稍微简单一些,可以看到虽然剪切形式没有变化,但是不同的转录本表达量差异很大。


leafcutter得到剪切形式


最后这个有可变剪切形式。


您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存