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

1、用LeafCutter探索转录组数据的可变剪切本教程，首发于 生信菜鸟团博客： HYPERLINK /2949.html /2949.html背景介绍目前主流的探究转录组数据的可变剪切的算法要么是基于estimateisoform ratios 或者 exon inclusion levels , 但是挑战还是 蛮多的，可变剪切本跟正常转录本重合的比例很大，技术误 差也是有的，依赖于基因现有的注释信息，既不准确，也不 完全。所以作者开发了LeafCutter工具。LeafCutterworkflow.First, short reads are mapped to the genome. Wh

2、en 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 clu

3、sters 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 数据集上面测试了，并且把结果跟GENCODEv19, Ensembl, and UCSC 着3大主

4、流的基因注释信息数据库比较。还在其它数据库里面验证了，数据下载地址是：dbGaP under accession phs000424.v6.p1 (GTEx), GEO under accession GSE41637 (RNA-seq data from mammalian organs), and ENA under accessionPRJEB3366 (Geuvadis).软件下载地址：LeafCuttersoftware,https:/davidaknowles/leafcutter;LeafViz visualizations,https:/leafcutter.shinyapps

5、.io/leafviz/;rheumatoid arthritis summary statistics, HYPERLINK http:/plaza.umin.ac.jp/yokada/datasource/software.htm http:/plaza.umin.ac.jp/yokada/datasource/software.htm. 软 件安装及使用最简单的就是conda进行安装了： condainstall -c davidaknowles r-leafcutter如果安装失败， 可能需要单独为它创建一个环境。不过，它本身就是一个R包，所以在个人电脑里面的rstudio里面安装即可

6、。if (!require(devtools) install.packages(devtools,repos=)devtools:install_github( davidaknowles/leafcutter/leafcutter)但是源代码里面有一些脚本和测试数据,所以还是要下载看看 mkdir -p /biosoft cd /biosoftgit clone HYPERLINK /davidaknowles/leafcuttercd /davidaknowles/leafcuttercd leafcutter# 需要修改里面的一个脚本scripts/bam2junc.sh把软件路径增添

7、进去即可里面又是 perl又是python的，感觉他们团队 开发环境不统一。第一步：bam2junc比对一般来说，优先选择STAR等支持跨越内含子的转录组比对工具得到bam文件，运行下面的脚本即可进行批量转换：cat bam_path.txt |while read iddofile=$(basename $id )sample=$file%.*echo Converting $id to $sample.junc sh /public/biosoft/leafcutter/scripts/bam2junc.sh $id $sample.juncdone 彳导至由勺junc文件如下:chr71

8、34840725134843893.1 -chr22343554422343557371 +chr4 37828435378315853910188210982755148465939827524575101851724733514913 +chr195 +chr1119 +chr188 -chr1212 -chr1514 -chr12 +chr10391017721097354454845873082751048510183232473231159292063192982445.1+这个步骤有点耗时，所有的junc文件地址需要保存给下一步使用第二步：Intron clustering这个步骤

9、，需要python2.7版本，这个是python的一个大坑，到现在版本仍然不统一。ls *.junc test_juncfiles.txtpython /public/biosoft/leafcutter/clustering/leafcutter_cluster.py -jtest_juncfiles.txt -m 50 -o testYRIvsEU -l 500000几分钟就运行完毕。得到的比较重要的文件如下：1.3M Jan 4 17:45testYRIvsEU_perind.counts.gz680K Jan4 17:45testYRIvsEU_perind_numers.counts

10、.gz5.0M Jan 4 17:45testYRIvsEU_pooled540K Jan 4 17:45testYRIvsEU_refined 877 Jan 4 17:45testYRIvsEU_sortedlibs 854 Jan 4 17:43 test_juncfiles.txt值得注意的是 testYRIvsEU_perind_numers.counts.gz 文件，里面每一行都是一个内含子，每一列都是一个样本，写 明了它们的表达值，这些数值就可以用来做可变剪切分析。# zcattestYRIvsEU_perind_numers.counts.gz |tailchr8:145651

11、155:145651305:clu_6538 21 14 19 8 9 013 33 0 0 4 0 5 8 12 0 12 34 15 0 0 1011chr8:145651155:145651409:clu_6538 1021 611 186190 294 284 681 89 222 57 257 363 694 807 523 44 469812 926 71 80 260214chr8:145652362:145653872:clu_6539 1265 694 13274 302 71 178 34 44 12 63 122 230 218 472 6 146 14211084 16

12、 14 83 46chr8:145652654:145653872:clu_6539 4824 56 0 26 0 13 0 2 5 2 0 3 19 17 0 2 8 64 0 0 30chr8:145652674:145653872:clu_6539 18 26 0 0 0 7 2 0 50 0 0 1 6 11 0 3 34 37 0 0 96chr8:146017525:146017630:clu_6540 2 3 44 0 2 12 4 0 00 22 5 9 10 2 0 1 9 11 0 0 1 0chr8:146017525:146017751:clu_6540 1067 67

13、1 620 41295 347 224 89 62 33 262 136 229 223 356 17 288 4801842 9 35 70 23chr8:146076780:146078224:clu_6541 183 0 0 17 17 8 0 0 3 2 3 16 6 12 0 4 45 29 9 0 102chr8:146076780:146078378:clu_6541 22 17 0 0 0 3 1 0 00 3 2 15 7 2 0 7 62 55 0 0 40chr8:146076780:146078757:clu_6541 10 1 16 0 12 52 00 11 0 2

14、4 9 27 3 0 0 7 0 28 0 0 2 0第三步：制作分组矩阵进行差异分析避免暴露我真实的项目，这里就给作者的示例 文件吧：RNA.NA18486_YRI.chr1.bamYRIRNA.NA18487_YRI.chr1.bamYRIRNA.NA18488_YRI.chr1.bamYRIRNA.NA18489_YRI.chr1.bamYRIRNA.NA18498_YRI.chr1.bamYRIRNA.NA06984_CEU.chr1.bamCEURNA.NA06985_CEU.chr1.bamCEURNA.NA06986_CEU.chr1.bamCEURNA.NA06989_CEU.

15、chr1.bamCEURNA.NA06994_CEU.chr1.bam CEU很简单的两歹U文件，说明每一个样本属于哪个组即可。/public/biosoft/leafcutter/scripts/leafcutter_ds.R-num_threads 4 -exon_file=/public/biosoft/leafcutter/leafcutter/data/gencod e19_exons.txt.gz testYRIvsEU_perind_numers.counts.gz group_info.txt 这里的 group_info.txt 就是自己制作好的分组矩阵。值得提醒的是，上面的

16、文件有且只能有2个分组,这样软件才知道怎么样去比较，如果自己的分组很多，可以考虑制作多个分组文件，运行多次。当然，上面的脚本已经 没有必要在linux服务器里面运行啦。既然有了内含子的表达矩阵，又有了分组信息，差异分析根本就不会消耗多少计 算资源，全部下载到自己的电脑里面去做吧。自己打开文件 /public/biosoft/leafcutter/scripts/leafcutter_ds.R 就明白 了 整个流程。也是几分钟就完成了全部结果。Runningdifferential splicing analysisDifferential splicing summary: statuses

17、Freq1 2 introns usedin =min_samples_per_intron samples 42521 sample with coverage,。 6231sample with coveragemin_coverage 9394Not enough valid samples 30475Success 2068Saving results.Loading exons from/Users/jmzeng/biosoft/leafcutter/leafcutter/data/gencode19 _exons.txt.gzAll done, exiting 得到的文件里面，需要详细了 解的是 leafcutter_ds_cluster_significance.txt 主要靠 自 己 看readme啦。第四步：可视化那些可变剪切也是包装好的 脚本。 /Users/jmzeng/biosoft/leafcutter/scripts/ds_plots.R-e/Users/jmzeng/biosoft/leafcutter/leafcutter/data/genco