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

健明 2021-07-14

展开全文

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

本教程，首发于生信菜鸟团博客：http://www./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./leafviz/;
rheumatoid arthritis summary statistics, http://plaza./yokada/datasource/software.htm.

软件安装及使用

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

conda install -c davidaknowles r-leafcutter

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

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

if (!require("devtools")) install.packages("devtools", repos='http://cran.us.')
devtools::install_github("davidaknowles/leafcutter/leafcutter")

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

mkdir -p ~/biosoft
cd ~/biosoft
git clone https://github.com/davidaknowles/leafcutter
cd leafcutter
## 需要修改里面的一个脚本 scripts/bam2junc.sh 把软件路径增添进去即可

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

第一步:bam2junc

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

cat bam_path.txt |while read id
do
file=$(basename $id )
sample=${file%%.*}
echo Converting $id to $sample.junc
sh /public/biosoft/leafcutter/scripts/bam2junc.sh $id $sample.junc
done

得到的junc文件如下:

chr7 134840725 134843893 . 1 -
chr2 234355442 234355737 . 1 +
chr4 37828435 37831585 . 13 +
chr19 39101772 39101882 . 5 +
chr11 109735445 109827551 . 19 +
chr18 48458730 48465939 . 8 -
chr12 82751048 82752457 . 12 -
chr15 51018323 51018517 . 14 -
chr1 247323115 247335149 . 2 +
chr10 92920631 92982445 . 1 +

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

第二步：Intron clustering

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

ls *.junc >test_juncfiles.txt
python /public/biosoft/leafcutter/clustering/leafcutter_cluster.py -j test_juncfiles.txt -m 50 -o testYRIvsEU -l 500000

几分钟就运行完毕。

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

1.3M Jan 4 17:45 testYRIvsEU_perind.counts.gz
680K Jan 4 17:45 testYRIvsEU_perind_numers.counts.gz
5.0M Jan 4 17:45 testYRIvsEU_pooled
540K Jan 4 17:45 testYRIvsEU_refined
877 Jan 4 17:45 testYRIvsEU_sortedlibs
854 Jan 4 17:43 test_juncfiles.txt

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

# zcat testYRIvsEU_perind_numers.counts.gz |tail
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
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
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
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
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
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
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
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
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
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

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

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

RNA.NA18486_YRI.chr1.bam YRI
RNA.NA18487_YRI.chr1.bam YRI
RNA.NA18488_YRI.chr1.bam YRI
RNA.NA18489_YRI.chr1.bam YRI
RNA.NA18498_YRI.chr1.bam YRI
RNA.NA06984_CEU.chr1.bam CEU
RNA.NA06985_CEU.chr1.bam CEU
RNA.NA06986_CEU.chr1.bam CEU
RNA.NA06989_CEU.chr1.bam CEU
RNA.NA06994_CEU.chr1.bam CEU

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

/public/biosoft/leafcutter/scripts/leafcutter_ds.R --num_threads 4 \
--exon_file=/public/biosoft/leafcutter/leafcutter/data/gencode19_exons.txt.gz \
testYRIvsEU_perind_numers.counts.gz group_info.txt

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

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

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

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

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

Running differential splicing analysis...
Differential splicing summary:
statuses Freq
1 <2 introns used in >=min_samples_per_intron samples 425
2 <=1 sample with coverage>0 62
3 <=1 sample with coverage>min_coverage 939
4 Not enough valid samples 3047
5 Success 2068
Saving results...
Loading exons from /Users/jmzeng/biosoft/leafcutter/leafcutter/data/gencode19_exons.txt.gz
All done, exiting

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

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

也是包装好的脚本。

/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