LeafCutter的介绍及用法
原理和介绍引用:生信菜鸟团
介绍:利用RNA-seq数据来定量特异剪切的一个软件。其核心思想是利用剪切reads(跨内含子的reads)来定量样本间内含子的使用差异。(对这个我不理解)。这个文件标榜的是不使用注释文件,因为现有的注释文件并不准确,同时也不利用外显子来检测(为什么?),可能是利用外显子来检测的需要注释文件吧
优点:容易检测新的内含子;可以建模比exonic PSI更复杂的剪切事件;避免了估计isoform丰度的难题;简单、计算效率高的算法可以扩展到100甚至1000秒的样本
Our approach obviates the need for transcript annotations and circumvents the challenges in estimating relative isoform or exon usage in complex splicing events
这个方法不需要转录组注释,并且避免在复杂剪切事件中用相关亚型和外显子来估计特异剪切这个挑战。
用法:自己到github上下载脚本
1.将bams转化为juncs
for bamfile in `ls /datasets/bam/*.bam`
do
echo Converting $bamfile to $bamfile.junc
sh sh/public/biosoft/leafcutter/scripts/bam2junc.sh $bamfile $bamfile.junc
echo $bamfile.junc >> test_juncfiles.txt
运行这个脚本即可
2.内含子聚类
python ../clustering/leafcutter_cluster.py -j test_juncfiles.txt -m 50 -o testYRIvsEU -l 500000
值得注意的是 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
很简单的两列文件,说明每一个样本属于哪个组即可
4.差异分析
/public/biosoft/leafcutter/scripts/leafcutter_ds.R--num_threads4\--exon_file=/public/biosoft/leafcutter/leafcutter/data/gencode19_exons.txt.gz \testYRIvsEU_perind_numers.counts.gz group_info.txt
5.可视化那些可变剪切
/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-f0.05
所有的可变剪切形式都会可视化在一张PDF图里面。