關(guān)于單細(xì)胞轉(zhuǎn)錄組轉(zhuǎn)錄因子的分析我們之前在單細(xì)胞系列講過R語言版本的妓盲,參考:跟著Cell學(xué)單細(xì)胞轉(zhuǎn)錄組分析(十二):轉(zhuǎn)錄組因子分析矩肩,但是R語言分析起來速度非常慢椅邓,如果你動(dòng)輒上萬的單細(xì)胞可能要運(yùn)行好幾周,這顯然不現(xiàn)實(shí)访娶。pySCENIC則很好的解決了這個(gè)問題,分析速度很快。官方教程參考:
https://pyscenic.readthedocs.io/en/latest/
一蛇券、軟件安裝
老樣子,還是先說一下安裝和分析文件的準(zhǔn)備樊拓,前面環(huán)境的配置和之前cellphonedb一樣纠亚,如果已經(jīng)操作過的,可以跳過:
#安裝下載及環(huán)境設(shè)置
# 安裝一個(gè)conda筋夏,為什么安裝他可以理解為Rstuido之于R,后期在環(huán)境設(shè)置蒂胞、軟件安裝上很方便。
wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh
# bash安裝叁丧,按照指引啤誊,都選yes,這樣一些依賴的python包都安裝了拥娄。
bash Miniconda3-latest-Linux-x86_64.sh
# 激活conda環(huán)境
source ~/.bashrc
#設(shè)置鏡像
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
然后就是安裝pySCENIC及環(huán)境設(shè)置:
#安裝pyscenic蚊锹,并創(chuàng)建分析環(huán)境
conda create -n pyscenic python=3.9#創(chuàng)建一個(gè)pyscenic 的python環(huán)境,pyscenic要求python版本3.6及以上稚瘾,目前python出到3.9了牡昆,我用3.9
conda activate pyscenic #激活pyscenic 環(huán)境
#安裝依賴包
conda install -y numpy
conda install -y -c anaconda cytoolz
conda install -y scanpy
#安裝pyscenic
pip install pyscenic
pip安裝即可,如果安裝失敗或者覺得下載速度太慢摊欠,可以使用下面代碼丢烘,使用與其他軟件的pip安裝。
pip install cellphonedb -i http://pypi.douban.com/simple/ --trusted-host pypi.douban.com
二些椒、數(shù)據(jù)文件準(zhǔn)備
一些分析用的依賴數(shù)據(jù)庫及文件:
#TF注釋
#鼠的
wget https://resources.aertslab.org/cistarget/motif2tf/motifs-v9-nr.mgi-m0.001-o0.0.tbl
#人的
wget https://resources.aertslab.org/cistarget/motif2tf/motifs-v9-nr.hgnc-m0.001-o0.0.tbl
#轉(zhuǎn)錄組因子列表
下載地址
https://github.com/aertslab/pySCENIC/tree/master/resources
#人的文件名:hs_hgnc_tfs.txt播瞳,復(fù)制為txt文件即可
#鼠的文件名:mm_mgi_tfs.txt,復(fù)制為txt文件即可
#reference數(shù)據(jù)庫免糕,之前一些網(wǎng)上教程的鏈接文件已經(jīng)不行了赢乓,因?yàn)樽隽烁掠遣啵艿臅r(shí)候會出錯(cuò),我是根據(jù)報(bào)錯(cuò)選擇了下面的文件
#鼠的
wget https://resources.aertslab.org/cistarget/databases/mus_musculus/mm10/refseq_r80/mc_v10_clust/gene_based/mm10_10kbp_up_10kbp_down_full_tx_clustered.genes_vs_motifs.rankings.feather
#人的
wget https://resources.aertslab.org/cistarget/databases/homo_sapiens/hg38/refseq_r80/mc_v10_clust/gene_based/hg38_10kbp_up_10kbp_down_full_tx_v10_clust.genes_vs_motifs.rankings.feather
三牌芋、分析文件準(zhǔn)備
這一步是從我們的R種seurat對象提取表達(dá)矩陣蚓炬,并轉(zhuǎn)化為loom文件。這里使用一個(gè)python腳本躺屁,轉(zhuǎn)化一下即可肯夏。
#文件準(zhǔn)備
#pyscenic的輸入文件是行為基因名,列為細(xì)胞ID的矩陣犀暑,所以在seurat對象中導(dǎo)出矩陣的時(shí)候需要轉(zhuǎn)置一下驯击,可以用標(biāo)準(zhǔn)化矩陣,也可以用counts矩陣母怜,影響不大余耽!
#表達(dá)矩陣、meta----R中進(jìn)行
write.csv(t(as.matrix(sce@assays$RNA@counts)),file = "sce_exp.csv")
#cellInfo <- sce@meta.data[,c("celltype","nCount_RNA","nFeature_RNA")]
#colnames(cellInfo) <- c('CellType', 'nGene' ,'nUMI')
#head(cellInfo)
#write.csv(cellInfo, file = "cellInfo.csv")
#轉(zhuǎn)化為loom文件苹熏,Linux下的python腳本
#編輯腳本
vim trans.py
#輸入以下內(nèi)容
import os, sys
os.getcwd()
os.listdir(os.getcwd())
import loompy as lp;
import numpy as np;
import scanpy as sc;
x=sc.read_csv("sce_exp.csv");#R中導(dǎo)出的表達(dá)矩陣
row_attrs = {"Gene": np.array(x.var_names),};
col_attrs = {"CellID": np.array(x.obs_names)};
lp.create("sce.loom",x.X.transpose(),row_attrs,col_attrs)
#保存并退出
#運(yùn)行trans.py
python trans.py
ls
#這樣在文件夾中會出現(xiàn)sce.loom文件碟贾,就是接下來輸入pyscenic的文件。
這樣整個(gè)準(zhǔn)備工作就完成了轨域,后續(xù)的分析很簡潔袱耽,只需要三步。然后就是可視化了干发,這個(gè)下回分解朱巨,覺得分享有用對你有幫助的,點(diǎn)個(gè)贊枉长、分享下再走唄冀续!記得關(guān)注下,之后不迷路必峰!更多精彩請?jiān)L問我的公眾號《KS科研分享與服務(wù)》