Py学习  »  Python

使用基于python的velocyto软件做RNA速率分析

生信技能树 • 2 年前 • 1280 次点击  

单细胞领域的一个高级分析是RNA速率分析,使用velocyto软件可以做,我们同样的把它区分为上下游分析。

上游分析需要在Linux操作环境里面,前面对10x的测序数据fq文件完成了 cellranger命令之后会有一个outputs文件夹。在该文件夹运行conda安装好的Python版本的velocyto软件即可,输出loom文件,供下游R里面操作。

安装自己的conda,每个用户独立操作

安装方法代码如下:

# 首先下载文件,20M/S的话需要几秒钟即可
wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh
# 接下来使用bash命令来运行我们下载的文件,记得是一路yes下去
bash Miniconda3-latest-Linux-x86_64.sh 
#  安装成功后需要更新系统环境变量文件
source ~/.bashrc

安装好conda后需要设置镜像。

conda config --add channels r 
conda config --add channels conda-forge 
conda config --add channels bioconda
conda config --add channels https://mirrors.bfsu.edu.cn/anaconda/cloud/bioconda/
conda config --add channels https://mirrors.bfsu.edu.cn/anaconda/cloud/conda-forge/
conda config --add channels https://mirrors.bfsu.edu.cn/anaconda/pkgs/free/
conda config --add channels https://mirrors.bfsu.edu.cn/anaconda/pkgs/main/
conda config --set show_channel_urls yes 

使用conda安装velocyto的一些依赖

暂不支持conda直接安装 velocyto,将会在1.0版本使用 。

#需要一些依赖
conda create -n velocyto 
conda activate velocyto 
conda install samtools
conda install numpy scipy cython numba matplotlib scikit-learn h5py click
pip install pysam
 # `PyPI`安装
pip install velocyto
# Successfully installed loompy-3.0.6 numpy-groupies-0.9.13 pandas-1.2.5 pytz-2021.1 velocyto-0.17.17

下载特定物种的特殊gtf文件

基因组注释文件主要是:GENCODE 或者Ensembl  ; 不过, 我们这个单细胞转录组使用cellranger流程的话,需要重复数据的gtf文件,rmsk

 

我们这里是 hg38_repeat_rmsk.gtf 文件,需要自己制作并且下载到指定工作目录!

从cellranger得到loom文件

这里需要 使用基于python的velocyto软件,它需要3个参数,其中两个是gtf文件,一个是前面cellranger命令的outputs目录哦,完整的命令如下:

rmsk_gtf=$HOME/pipeline/velocyto/hg38_repeat_rmsk.gtf # 从genome.ucsc.edu下载 
cellranger_outDir=HSY-fushui # 前面cellranger命令的outputs目录 
cellranger_gtf=$HOME/pipeline/refdata-gex-GRCh38-2020-A/genes/genes.gtf # 这个是cellranger官网提供的

ls -lh $rmsk_gtf  $cellranger_outDir $cellranger_gtf

nohup velocyto run10x -m $rmsk_gtf  $cellranger_outDir $cellranger_gtf & 
# 如果是其它单细胞数据,可以换参数,比如run_smartseq2   

理论上,这个一句话命令,就可以完成了python的velocyto软件并且拿到loom文件,但是,因为 samtools问题,上面的流程可能是会失败。

这个时候可以自行先运行 samtools sort 命令处理得到 cellsorted_possorted_genome_bam.bam  文件,如下所示的命令:

cd $cellranger_outDir/out 
nohup samtools sort -@ 10  -t CB -O BAM -o cellsorted_possorted_genome_bam.bam possorted_genome_bam.bam &

# 这个  samtools sort  的 速度很快

# 下次就可以把前面的 velocyto run10x 重新跑一次,因为  samtools sort 命令已经是成功了。

# The file /home/---/outs/cellsorted_possorted_genome_bam.bam already exists. 
# The sorting step w ill be skipped and the existing file will be used.

在同样的前面cellranger命令的outputs  文件夹目录下面,会输出一个 velocyto 文件夹,如下所示:


$ ls -lh  velocyto/ 
158M 7月   3 22:38 HSY-fushui.loom

有了这个loom文件,后续就是回到R语言里面的统计可视化啦!

如果是多个10x样品都需要运行velocyto

首先,把全部的bam文件循环进行sort,代码如下:

 ls  */outs/possorted_genome_bam.bam|while read id;do  new=${id/possorted_genome_bam.bam/cellsorted_possorted_genome_bam.bam}
echo $new 
nohup samtools sort -@ 4  -t CB -O BAM -o $new   $id  & 
 done

这个代码难度有点高哦,需要精通Linux才能理解它。

同样的,然后可以批量走python的velocyto软件,代码如下所示:

rmsk_gtf=$HOME/pipeline/velocyto/hg38_repeat_rmsk.gtf # 从genome.ucsc.edu下载 
#cellranger_outDir=HSY-fushui # 前面cellranger命令的outputs目录 
cellranger_gtf=$HOME/pipeline/refdata-gex-GRCh38-2020-A/genes/genes.gtf # 这个是cellranger官网提供的
#ls -lh $rmsk_gtf  $cellranger_outDir $cellranger_gtf

# 同样的一个简单的 循环即可
ls -d *-*|while read cellranger_outDir;do 
nohup velocyto run10x -m $rmsk_gtf  $cellranger_outDir $cellranger_gtf & 
done 

每个样品都会输出如下所示的loom文件:


158M 7月  19 17:39 HSY-fushui.loom
 79M 7月  19 17:39 HSY-PBMC.loom
132M 7月  19 17:39 HSY-yi.loom
163M 7月  19 17:39 HSY-yuan.loom
246M 7月  19 17:39 LS-Endo-Pro.loom
139M 7月  19 17:39 LS-PBMC-Pro.loom
211M 7月  19 17:39 RDX-PBMC.loom
205M 7月  19 17:39 RDX-yuan.loom
123M 7月  19 17:39 WY-PBMC.loom
 91M 7月  19 17:39 WY-yi.loom
136M 7月  19 17:39 WY-yuan.loom
253M 7月  19 17:39 YF-fushui.loom
188M 7月  19 17:39 YF-PBMC.loom
135M 7月  19 17:39 YF-yi.loom
148M 7月  19 17:39 YF-yuan.loom
243M 7月  19 17:39 YX-Endo-Decidu.loom
128M 7月  19 17:39 YX-PBMC-Decidu.loom
 95M 7月  19 17:39 ZZX-PBMC.loom
123M 7月  19 17:39 ZZX-yuan-2.loom

拿到了这么多 loom文件,就可以进行下游 velocyto.R 这个R包进行后续统计可视化啦!明天我们就讲解velocyto.R 这个R包用法!

关于RNA velocity (gene expression trajectory)

RNA velocity是基于真实的转录动力学,可用于细胞基因表达的动态分化的研究。

如上左图,刚转录出的mRNA包含外显子和内含子,经过splicing切除内含子后,得到用于编码 蛋白的spliced mRNA。spliced mRNA的丰度由未成熟mRNA的splicing速度和降解速率共同决 定。如上中图:每个点代表一个细胞,在拟时间轴上,未经过剪切的mRNA的出现始终早于经 过剪切的mRNA。如上右图:红色代表未经过剪切的mRNA,蓝色代表经过剪切的mRNA,可以 看出,这些细胞的应该是从左往右分化的,因此Velocity可以用于定义轨迹的起点分支和终点。也就是说,Velocity可以在不知发育过程的前提下,预测谱系的方向(如下图)。

Velocity可以用于周期的轨迹

号外:  我们提供单细胞数据分析服务哦

Python社区是高质量的Python/Django开发社区
本文地址:http://www.python88.com/topic/117303
 
1280 次点击