ise生成msk文件的用处_RNA velocity分析练习（三）生成loom文件

Velocyto软件针对不同测序平台，有不同的方法进行loom文件的提取，你可以参考官网：here来进行操作。

我练习的数据是Smartseq2平台测序结果，所以这里只介绍这一种。

velocyto includes a shortcut to perform the read counting for UMI-less, not stranded, full-length techniques such as SmartSeq2.

根据我之前写的笔记，你需要准备好基因组的annotation文件，和repeat_msk.gtf文件。这一步所需时间取决于你数据的测序深度和电脑配置，但一般不会超过６小时。(实际上我的破电脑真的跑了6个小时。。。)

已经准备好的bam文件：

下载好的基因组注释文件(小鼠)：

我这是去年下载的，现在已经更新到v25了

下载并解压好的repeat_msk.gtf文件：

(一)生成loom文件

进入bam文件所在的文件夹：

$ velocyto run_smartseq2 -o /media/yanfang/FYWD/scRNA_seq/RNA_relocity/loom_files/ -m /media/yanfang/FYWD/RNA_seq/ref_genome/mm10_repeat_msk.gtf -e MyTissue /media/yanfang/FYWD/scRNA_seq/RNA_relocity/GSE99933_bam/*.bam /media/yanfang/FYWD/RNA_seq/ref_genome/gencode.vM22.annotation.gtf

(代码里的"MyTissue"这个名字是一会儿生成loom文件的前缀)

-o ：输出文件的文件夹

run_smartseq2：指定哪个平台的测序结果

运行过程会弹出很多信息，每读取一个bam文件都会显示几行信息，例如说：

2020-05-19 17:13:10,863 - DEBUG - Reading /media/yanfang/FYWD/scRNA_seq/RNA_relocity/GSE99933_bam/E13.5_P9.bam

2020-05-19 17:13:11,197 - DEBUG - Read first 0 million reads

2020-05-19 17:13:22,762 - DEBUG - Counting for batch 768, containing 1 cells and 419294 reads

2020-05-19 17:13:37,676 - DEBUG - 42428 reads in repeat masked regions

2020-05-19 17:13:37,676 - DEBUG - 191333 reads overlapping with features on plus strand

2020-05-19 17:13:37,677 - DEBUG - 188152 reads overlapping with features on minus strand

2020-05-19 17:13:37,677 - DEBUG - 25321 reads overlapping with features on both strands

2020-05-19 17:13:41,095 - WARNING - The barcode selection mode is off, no cell events will be identified by <80 counts

2020-05-19 17:13:41,095 - WARNING - 0 of the barcodes where without cell

bam文件都读取完成会显示：

2020-05-19 17:13:41,161 - DEBUG - Counting done!

2020-05-19 17:13:41,161 - DEBUG - Example of barcode: E13.5_D9.bam and cell_id: MyTissue:E13.5_D9.bam

2020-05-19 17:13:41,161 - DEBUG - Generating output file /media/yanfang/FYWD/scRNA_seq/RNA_relocity/loom_files/MyTissue.loom

2020-05-19 17:13:41,161 - DEBUG - Collecting row attributes

2020-05-19 17:13:41,324 - DEBUG - Generating data table

2020-05-19 17:13:44,062 - DEBUG - Writing loom file

2020-05-19 17:13:53,902 - DEBUG - Terminated Succesfully!

得到的loom文件：

(二)读取loom文件

这里用R来进行举例：

> library(devtools)

> install_github("velocyto-team/velocyto.R")

> library(velocyto.R)

#load data

> input_loom

> adata

可以来看一下adata这个对象里都有什么：

其中spliced和unspliced就是我们需要的成熟的mRNA和未成熟mRNA的数值；spanning的意思是落在intron+exon上的reads count

# Use the spliced data as the input data提取spliced数据

> spliced_adata

> dim(spliced_adata)

[1] 55487 768　＃5万多个基因，768个细胞

> path

> files

> names(files)

# parse gene annotation, annotate bam file reads

> dat

到这一步，我们所需要的count矩阵就有了。后面就是RNA velocity的分析过程了。