用Expedition来分析单细胞转录组数据的可变剪切

了解我的应该都知道我最近几个月都在奋战一个陌生的领域，单细胞转录组数据处理。真的很有挑战性，笔记累积了一大堆了，但是没有太值得分享的，大多是利用bulk转录组数据处理的经验而已，但是下面这个是单细胞转录组数据独有的，简单分享一下吧。

工具发表于2017，文章是 Single-cell alternative splicing analysis with Expedition reveals splicing dynamics during neuron differentiation  为了展示他们的工具，测了各种细胞，包括 62 iPSCs, 69 NPCs, and 60 MNs 数据库是可以在GSE85908下载 ：

• induced pluripotent stem cells (iPSCs)

• neural progenitor cells (NPCs)

• motor neurons(MNs)

软件地址：(http://github.com/YeoLab/Expedition)  其实是一个软件合集，包括下面的3个软件需要分别安装分别使用。

• (i) outrigger, a de novo splice graph transversal algorithm to detect AS from single cell RNA-seq;

• (ii) anchor, a Bayesian approach to assign splicing modalities

• (iii) bonvoyage, using non-negative matrix factorization to visualize modality changes.

当然，主要是第一个软件的运行及解读，后面的两个锦上添花，可以忽略，不影响分析单细胞转录组数据的可变剪切的基本需求。

值得一提的是，该软件只关注mutually exclusive exon (MXE)和skipped exon (SE)这两种 alternative splicing (AS) events ！

一、outrigger

软件github地址是：https://github.com/YeoLab/outrigger

软件说明书在：http://yeolab.github.io/outrigger/

可以用conda非常方便的安装

conda config --add channels r
conda config --add channels bioconda
conda create --name outrigger-env outrigger
source activate outrigger-env

使用方法也非常简单，说明书讲解的非常清楚，就三个步骤，针对比对好的bam文件(需要STAR软件比对，才能输出splice junction (SJ.out.tab)文件)

cd ~/projects/tasic2016/analysis/tasic2016_v1
outrigger index --sj-out-tab *SJ.out.tab \
   --gtf /projects/ps-yeolab/genomes/mm10/gencode/m10/gencode.vM10.annotation.gtf
outrigger validate --genome mm10 \
   --fasta /projects/ps-yeolab/genomes/mm10/GRCm38.primary_assembly.genome.fa
outrigger psi

看懂算法需要一点点时间和耐心，https://pypi.python.org/pypi/outrigger， 图文并茂的讲解了。但是作者说这个软件有个很严重的缺点，太耗时间

• outrigger index: This will run for 24-48 hours.

• outrigger validate: This will take 2-4 hours.

• outrigger psi: This will run for 4-8 hours.

起初为了加快软件运行速度，我给定了 --n-jobs 5 多线程，25G的内存，但是内存爆表了，重新分配50G内存，一晚上终于运行完毕了。

得到的结果如下：

outrigger_output/
├── index
│   ├── exon_direction_junction.csv
│   ├── gtf
│   │   ├── gencode.vM12.annotation.gtf
│   │   ├── gencode.vM12.annotation.gtf.db
│   │   └── novel_exons.gtf
│   ├── mxe
│   │   ├── event.bed
│   │   ├── events.csv
│   │   ├── exon1.bed
│   │   ├── exon2.bed
│   │   ├── exon3.bed
│   │   ├── exon4.bed
│   │   ├── intron.bed
│   │   ├── splice_sites.csv
│   │   └── validated
│   │       └── events.csv
│   └── se
│       ├── event.bed
│       ├── events.csv
│       ├── exon1.bed
│       ├── exon2.bed
│       ├── exon3.bed
│       ├── intron.bed
│       ├── splice_sites.csv
│       └── validated
│           └── events.csv
├── junctions
│   ├── metadata.csv
│   └── reads.csv
└── psi
   ├── mxe
   │   ├── psi.csv
   │   └── summary.csv
   ├── outrigger_psi.csv
   ├── outrigger_summary.csv
   └── se
       ├── psi.csv
       └── summary.csv

10 directories, 29 files

二、anchor

软件github地址是：https://github.com/YeoLab/anchor

软件说明书在： https://yeolab.github.io/anchor/

可以用conda非常方便的安装

conda create -n anchor-env pandas scipy numpy matplotlib seaborn
pip install anchor-bio

使用方法也有点诡异，看起来需要进入python的交互界面

import anchor

bm = anchor.BayesianModalities()
modalities = bm.fit_transform(data)

还没看懂干什么的。

三、bonvoyage

软件github地址是： https://github.com/YeoLab/bonvoyage

软件说明书在： http://yeolab.github.io/bonvoyage/

可以用conda非常方便的安装

conda create -n anchor-env pandas numpy matplotlib seaborn scikit-learn
pip install bonvoyage

这个主要是用来可视化上面步骤的anchor推断的可变剪切形式

import bonvoyage

wp = bonvoyage.Waypoints()
waypoints = wp.fit_transform(data)
import bonvoyage

bonvoyage.waypointplot(waypoints)

背景知识

mutually exclusive exon (MXE)skipped exon (SE)这些名词的解释我就不搬运了，但是下面这个原理图还是值得瞧一瞧：

总的来说，你看完这个教程应该是没办法学到什么技巧，只能收藏一下咯，毕竟大部分人根本木有单细胞数据，也没有足够强大的服务器来hold住这样的大数据。

但是，万一，你也正好在做方面的探索，欢迎来信跟我讨论具体细节。