社区所有版块导航
Python
python开源   Django   Python   DjangoApp   pycharm  
DATA
docker   Elasticsearch  
aigc
aigc   chatgpt  
WEB开发
linux   MongoDB   Redis   DATABASE   NGINX   其他Web框架   web工具   zookeeper   tornado   NoSql   Bootstrap   js   peewee   Git   bottle   IE   MQ   Jquery  
机器学习
机器学习算法  
Python88.com
反馈   公告   社区推广  
产品
短视频  
印度
印度  
Py学习  »  Python

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

生信技能树 • 5 年前 • 1958 次点击  

单细胞领域的一个高级分析是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