RNA velocity分析练习(一)文件下载以及预处理

目前有多种方法可以进行RNA velocity的分析:

  1. scVelo(参考文章:单细胞转录组数据分析|| scVelo 教程:RNA速率分析工具
  2. velocyto (官网:Welcome to velocyto.py!
  3. 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文件需要哪些软件。

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
禁止转载,如需转载请通过简信或评论联系作者。
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 194,524评论 5 460
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 81,869评论 2 371
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 141,813评论 0 320
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 52,210评论 1 263
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 61,085评论 4 355
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 46,117评论 1 272
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 36,533评论 3 381
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 35,219评论 0 253
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 39,487评论 1 290
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 34,582评论 2 309
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 36,362评论 1 326
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 32,218评论 3 312
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 37,589评论 3 299
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 28,899评论 0 17
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 30,176评论 1 250
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 41,503评论 2 341
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 40,707评论 2 335