用LeafCutter探索轉(zhuǎn)錄組數(shù)據(jù)的可變剪切

1、用LeafCutter探索轉(zhuǎn)錄組數(shù)據(jù)的可變剪切本教程，首發(fā)于 生信菜鳥團(tuán)博客： HYPERLINK /2949.html /2949.html背景介紹目前主流的探究轉(zhuǎn)錄組數(shù)據(jù)的可變剪切的算法要么是基于estimateisoform ratios 或者 exon inclusion levels , 但是挑戰(zhàn)還是 蠻多的，可變剪切本跟正常轉(zhuǎn)錄本重合的比例很大，技術(shù)誤 差也是有的，依賴于基因現(xiàn)有的注釋信息，既不準(zhǔn)確，也不 完全。所以作者開發(fā)了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 數(shù)據(jù)集上面測試了，并且把結(jié)果跟GENCODEv19, Ensembl, and UCSC 著3大主

4、流的基因注釋信息數(shù)據(jù)庫比較。還在其它數(shù)據(jù)庫里面驗(yàn)證了，數(shù)據(jù)下載地址是：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進(jìn)行安裝了： condainstall -c davidaknowles r-leafcutter如果安裝失敗， 可能需要單獨(dú)為它創(chuàng)建一個環(huán)境。不過，它本身就是一個R包，所以在個人電腦里面的rstudio里面安裝即可

6、。if (!require(devtools) install.packages(devtools,repos=)devtools:install_github( davidaknowles/leafcutter/leafcutter)但是源代碼里面有一些腳本和測試數(shù)據(jù),所以還是要下載看看 mkdir -p /biosoft cd /biosoftgit clone HYPERLINK /davidaknowles/leafcuttercd /davidaknowles/leafcuttercd leafcutter# 需要修改里面的一個腳本scripts/bam2junc.sh把軟件路徑增添

7、進(jìn)去即可里面又是 perl又是python的，感覺他們團(tuán)隊(duì) 開發(fā)環(huán)境不統(tǒng)一。第一步：bam2junc比對一般來說，優(yōu)先選擇STAR等支持跨越內(nèi)含子的轉(zhuǎn)錄組比對工具得到bam文件，運(yùn)行下面的腳本即可進(jìn)行批量轉(zhuǎn)換：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 彳導(dǎo)至由勺junc文件如下:chr71

8、34840725134843893.1 -chr22343554422343557371 +chr4 37828435378315853910188210982755148465939827524575101851724733514913 +chr195 +chr1119 +chr188 -chr1212 -chr1514 -chr12 +chr10391017721097354454845873082751048510183232473231159292063192982445.1+這個步驟有點(diǎn)耗時，所有的junc文件地址需要保存給下一步使用第二步：Intron clustering這個步驟

9、，需要python2.7版本，這個是python的一個大坑，到現(xiàn)在版本仍然不統(tǒng)一。ls *.junc test_juncfiles.txtpython /public/biosoft/leafcutter/clustering/leafcutter_cluster.py -jtest_juncfiles.txt -m 50 -o testYRIvsEU -l 500000幾分鐘就運(yùn)行完畢。得到的比較重要的文件如下：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 文件，里面每一行都是一個內(nèi)含子，每一列都是一個樣本，寫 明了它們的表達(dá)值，這些數(shù)值就可以用來做可變剪切分析。# 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第三步：制作分組矩陣進(jìn)行差異分析避免暴露我真實(shí)的項(xiàng)目，這里就給作者的示例 文件吧：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個分組,這樣軟件才知道怎么樣去比較，如果自己的分組很多，可以考慮制作多個分組文件，運(yùn)行多次。當(dāng)然，上面的腳本已經(jīng) 沒有必要在linux服務(wù)器里面運(yùn)行啦。既然有了內(nèi)含子的表達(dá)矩陣，又有了分組信息，差異分析根本就不會消耗多少計(jì) 算資源，全部下載到自己的電腦里面去做吧。自己打開文件 /public/biosoft/leafcutter/scripts/leafcutter_ds.R 就明白 了 整個流程。也是幾分鐘就完成了全部結(jié)果。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 得到的文件里面，需要詳細(xì)了 解的是 leafcutter_ds_cluster_significance.txt 主要靠 自 己 看readme啦。第四步：可視化那些可變剪切也是包裝好的 腳本。 /Users/jmzeng/biosoft/leafcutter/scripts/ds_plots.R-e/Users/jmzeng/biosoft/leafcutter/leafcutter/data/genco