目前有多种方法可以进行RNA velocity的分析：
- scVelo（参考文章：单细胞转录组数据分析|| scVelo 教程：RNA速率分析工具）
- velocyto （官网:Welcome to velocyto.py!）
- Seurat （参考文章：用Seurat做RNA Velocity）
在前一篇的文献学习里（RNA velocity of single cells文献学习），作者使用的是velocyto软件，也就是上面的第二个软件进行分析的，所以我也主要学习这个软件的使用。作者的详细代码见：here，里面有R和python两种版本。
由于RNA velocity分析的前提，是要我们从单细胞RNA-seq的数据里区分出未成熟mRNA（unspliced）和成熟的mRNA（spliced），所以你需要从fastq文件开始，与基因组进行对比后得到sam文件，从sam文件转成bam，再从bam文件里提取这些信息，最后你会得到.loom为后缀的文件，这个文件才是我们需要的。
你可以在这个网站(http://velocyto.org/velocyto.py/tutorial/cli.html)里学习如何从bam文件里提取我们需要的信息，你可以使用10X、SMART-seq2、dropseq等平台测序得到的fastq文件，每一种情况都有详细的说明。
这次练习的原始数据是GSE99933，有768个文件，这个project是应用SMART-Seq2平台进行测序的。
懒得从fastq文件走一遍完整流程的童鞋可以直接下载bam文件：
http://pklab.med.harvard.edu/velocyto/chromaffin/bams.tar
或者你可以直接下载loom文件进行分析：
http://pklab.med.harvard.edu/velocyto/chromaffin/dat.rds
（一）fastq文件下载以及文件名的批量修改
关于如何批量下载fastq文件，如何比对，如何从sam文件转成bam文件，我在单细胞测序实战（第一部分）写的非常详细。这里一共768个fastq文件，大约是12G左右。
（其实批量修改文件名不是必要的，完全可以跳过这一步，但是不知道为什么自己突发奇想要改名，可能是看着SRR文件名不舒服吧。。。想节约时间的童鞋可以直接跳到fastqc步骤）
下载fastq文件后，文件名都是SRR开头的，那么我想把前缀改成细胞的编号：E13.5_E之类的，应该怎么弄？首先你要先拿到细胞编号的list，在GEO页面里：
点击“Accession list truncated, click here to browse through all related public accessions”
在新页面里，点击“Export”:
选择All search results
现在你下载了所有样品的信息，是一个csv表，用excel打开是这样的：
用查找/替换功能把[single cell RNA-seq]替换成空白（注意是空白，不是空格）：
从这个表里提取第二列的细胞编号：
$ awk -F"," '{print $2}' sample.csv>> cell_name.txt
$ head cell_name.txt
E12.5_A1
E12.5_A2
E12.5_A3
E12.5_A4
E12.5_A5
E12.5_A6
E12.5_A7
E12.5_A8
E12.5_A9
E12.5_A10
在win10系统下批量修改文件名，看视频：here或者windows如何批量修改文件名
修改后的fastq文件的名字就变成了：
（二）fastqc
随便抽取几个fastq文件看一下：
整体的测序质量看起来都不错
（三）比对
进入fastq文件夹里进行比对：
$ ls *.gz | while read i;do hisat2 -p 10 -x /media/yanfang/FYWD/RNA_seq/ref_genome/index/mm10/genome -U $i -S /media/yanfang/FYWD/scRNA_seq/RNA_relocity/GSE99933_sam/${i}.sam;done
这一步要过夜（大概10个小时多一点），因为是用自己的电脑跑的，所以比较慢。生成的sam文件：
（四）生成bam文件
可以先看一下后续分析要去的bam文件是什么样的：
需要的bam文件是要sort后的，所以一定要注意！
$ ls *.fastq.gz.sam | while read i ;do (samtools sort -O bam -@ 10 -o /media/yanfang/FYWD/scRNA_seq/RNA_relocity/GSE99933_bam/$(basename ${i} ".fastq.gz.sam").bam ${i});done
（五）下载mouse_repeat_msk.gtf
除了mm10_annotation.gtf以外(可在https://www.gencodegenes.org/下载)，我们还需要一个文件进行后续的分析，这个文件叫repeat_msk.gtf，你可以在这个网站：here下载到：
选项都填好，点击左下角“get output”。就行了。
这个gtf文件是帮助我们去除表达的重复元素(expressed repetitive elements)，因为这些reads可能在下游分析中构成一个混淆因素。
到此，生成loom文件的必需文件就都准备齐全了。下一篇会介绍生成loom文件需要哪些软件。
网友评论