目前有多種方法可以進(jìn)行RNA velocity的分析：
- scVelo（參考文章：單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)分析|| scVelo 教程：RNA速率分析工具）
- velocyto （官網(wǎng):Welcome to velocyto.py!）
- Seurat （參考文章：用Seurat做RNA Velocity）
在前一篇的文獻(xiàn)學(xué)習(xí)里（RNA velocity of single cells文獻(xiàn)學(xué)習(xí)）刁赖，作者使用的是velocyto軟件蝶棋，也就是上面的第二個(gè)軟件進(jìn)行分析的理朋，所以我也主要學(xué)習(xí)這個(gè)軟件的使用。作者的詳細(xì)代碼見：here，里面有R和python兩種版本。
由于RNA velocity分析的前提，是要我們從單細(xì)胞RNA-seq的數(shù)據(jù)里區(qū)分出未成熟mRNA（unspliced）和成熟的mRNA（spliced）孤个，所以你需要從fastq文件開始，與基因組進(jìn)行對(duì)比后得到sam文件沛简，從sam文件轉(zhuǎn)成bam齐鲤，再從bam文件里提取這些信息，最后你會(huì)得到.loom為后綴的文件椒楣，這個(gè)文件才是我們需要的给郊。
你可以在這個(gè)網(wǎng)站(http://velocyto.org/velocyto.py/tutorial/cli.html)里學(xué)習(xí)如何從bam文件里提取我們需要的信息，你可以使用10X捧灰、SMART-seq2淆九、dropseq等平臺(tái)測(cè)序得到的fastq文件，每一種情況都有詳細(xì)的說明。
這次練習(xí)的原始數(shù)據(jù)是GSE99933炭庙，有768個(gè)文件饲窿，這個(gè)project是應(yīng)用SMART-Seq2平臺(tái)進(jìn)行測(cè)序的。
懶得從fastq文件走一遍完整流程的童鞋可以直接下載bam文件：
http://pklab.med.harvard.edu/velocyto/chromaffin/bams.tar
或者你可以直接下載loom文件進(jìn)行分析：
http://pklab.med.harvard.edu/velocyto/chromaffin/dat.rds
（一）fastq文件下載以及文件名的批量修改
關(guān)于如何批量下載fastq文件焕蹄，如何比對(duì)逾雄，如何從sam文件轉(zhuǎn)成bam文件，我在單細(xì)胞測(cè)序?qū)崙?zhàn)（第一部分）寫的非常詳細(xì)擦盾。這里一共768個(gè)fastq文件嘲驾，大約是12G左右淌哟。
（其實(shí)批量修改文件名不是必要的迹卢，完全可以跳過這一步，但是不知道為什么自己突發(fā)奇想要改名徒仓，可能是看著SRR文件名不舒服吧腐碱。。掉弛。想節(jié)約時(shí)間的童鞋可以直接跳到fastqc步驟）
下載fastq文件后症见，文件名都是SRR開頭的，那么我想把前綴改成細(xì)胞的編號(hào)：E13.5_E之類的殃饿，應(yīng)該怎么弄谋作？首先你要先拿到細(xì)胞編號(hào)的list，在GEO頁面里：
在新頁面里乎芳，點(diǎn)擊“Export”:
現(xiàn)在你下載了所有樣品的信息遵蚜，是一個(gè)csv表，用excel打開是這樣的：
用查找/替換功能把[single cell RNA-seq]替換成空白（注意是空白奈惑，不是空格）：
從這個(gè)表里提取第二列的細(xì)胞編號(hào)：
$ 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系統(tǒng)下批量修改文件名吭净，看視頻：here或者windows如何批量修改文件名
修改后的fastq文件的名字就變成了：
（二）fastqc
隨便抽取幾個(gè)fastq文件看一下：
（三）比對(duì)
進(jìn)入fastq文件夾里進(jìn)行比對(duì)：
$ 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個(gè)小時(shí)多一點(diǎn)），因?yàn)槭怯米约旱碾娔X跑的肴甸，所以比較慢寂殉。生成的sam文件：
（四）生成bam文件
可以先看一下后續(xù)分析要去的bam文件是什么樣的：
$ 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/下載)，我們還需要一個(gè)文件進(jìn)行后續(xù)的分析赚爵，這個(gè)文件叫repeat_msk.gtf庶柿，你可以在這個(gè)網(wǎng)站：here下載到：
選項(xiàng)都填好村怪，點(diǎn)擊左下角“get output”。就行了澳泵。
這個(gè)gtf文件是幫助我們?nèi)コ磉_(dá)的重復(fù)元素(expressed repetitive elements)实愚，因?yàn)檫@些reads可能在下游分析中構(gòu)成一個(gè)混淆因素。
到此，生成loom文件的必需文件就都準(zhǔn)備齊全了腊敲。下一篇會(huì)介紹生成loom文件需要哪些軟件击喂。
