本日写一篇关于单细胞转录因子的分析-scenic
在做分析之前你应该看一下这篇文献，看看别人做了什么scenic分析
《Single-cell RNA sequencing highlights the role of inflammatory cancer-associated fibroblasts in bladder urothelial carcinoma》
分析原理
对于这个包我想发表一下感想，这个包刚开始对于我来说比力难理解，其实跟着官方教程running无非就是debug而已，这没啥，紧张的是这个分析做了些啥，末了得到啥，我能讲啥，这是最难的，反正我是不敢做那些我本身都不相识过程的分析，以是起首我们要相识这个包的分析流程：
GENIE3或GRNBoost→RcisTarget→AUCell（以是起首你得先相识这三个包做了啥，不外不太想相识也可以以后看，我会讲讲我的理解，盼望对你们有所资助）
1.利用GENIE3或GRNBoost(Gradient Boosting)基于共表达推断转录因子与候选靶基因之间的共表达模块
https://www.jianshu.com/p/d71dcd4cff5a (可以学一下这个包，体验一下，大概我理解的不对哈）
是不是看到这个就头大 ，所有的教程都是如许，其实简朴的说第一步推断出转录因子(TF)与其作用的靶点/gene(一个转录因子大概调控很多的gene哈)
2.RcisTarget对每个共表达模块举行顺式调控基序(cis-regulatory motif）分析
https://www.jianshu.com/p/c9df97f9e6e0 (RcisTarget包）
由于GENIE3模子只是基于共表达，会存在很多假阳性和间接靶标，为了辨认直接团结靶标(direct-binding targets)，利用RcisTarget对每个共表达模块举行顺式调控基序(cis-regulatory motif）分析。举行TF-motif富集分析，辨认直接靶标。(仅生存具有精确的上游调治子且显着富集的motif modules，并对它们举行修剪以撤除缺乏motif支持的间接靶标。）这些处理惩罚后的每个TF及其潜伏的直接targets gene被称作一个调治子(regulon)，以是第二步就得到了转录因子(TF)与其对应的直接作用的靶点(targets genes),而且称为regulon(以是每一个regulon都是1个TF和这个TF调控的的targets genes)
3.利用AUCell算法对每个细胞的每个regulon活性举行打分
https://www.jianshu.com/p/6a6705f12842（AUCell包）
对于一个regulon来说，比力细胞间的AUCell得分可以判断出该regulon在哪种细胞显着更高的，这将确定regulon在哪些细胞中处于"打开"状态。简朴来说第三步得到了每一个细胞对每一个regulon的得分，通过得分我们可以知道每一个细胞的每一个regulon的激活情况，即TF的激活情况
以是综上所得，我的浅近的认知：SCENIC能让我们知道每一个细胞他们差别的TF激活情况，不要跟RcisTarget包弄混了，固然都是得到都是TF的差别，但是RcisTarget包是通过二者的差别基因来判断是什么TF造成的，而SCENIC是可以大概得出每一个细胞的TFs的激活情况。
ok 关于分析流程就将到这，固然这只是我个人浅近的理解，更尺度，正直表明的还是要看原文献和各个生物信息学大佬的表明哈。（接待各人讨论与品评）
分析流程（R+linux）
我用过单纯R跑，超等慢，而且跑不出来，以是看别人的教程写着用linux+R一起
1.在R语言中获得单细胞的表达矩阵 scRNA是我本身的Seurat对象哈
write.csv(t(as.matrix(scRNA@assays$RNA@counts)), file = "scRNA.csv")#提取表达矩阵，并定名为scRNA.csv，留意矩阵肯定要转置，否则会报错2.从这里就要转战场去linux了，莫慌，学就是了
个人的linux学习流程：
1.https://www.bilibili.com/video/BV1pE411C7ho?p=42&spm_id_from=333.880.my_history.page.click
2.https://www.bilibili.com/video/BV1ds411g7eg?p=7&spm_id_from=333.880.my_history.page.click
发起先学第一个视频再学jimmy老师的，由于jimmy老师讲的都是干货，但是呢比力散，发起先把根本视频学一遍，在听jimmy老师的，会更加有印象。
当你学完了就可以run接下来的流程了，
从如今开始都在linux运行
起首你要在Linux安装conda，（这是第一道难关）conda是什么，怎么安装，它醒目什么：
conda是什么：我的理解，conda就像手机上的应用超市，我们在手机上进入应用超市点击下载就把app下载了，假如没有应用超市的话，我们要安装一个app则必要下载安装包再解压之类的，因此有了应用超市，我们可以更便捷的安装app，那么conda就像服务器上的应用超市，让我们更方便的下载各种软件和包。
conda安装：
#下载安装包，解压，激活wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh #下载安装包bash Miniconda3-latest-Linux-x86_64.sh #利用bash下令解压安装，记得是一起yes下去source ~/.bashrc #激活安装的conda#设置国内镜像conda config --add channels r conda config --add channels conda-forge conda config --add channels biocondaconda 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醒目什么：1.conda可以大概更好的安装软件与r包 2.conda能创建一个情况
创建情况是什么意思呢，起重要搞清晰创建情况跟mkdir不一样，对于情况的理解就是，雷同于我们看书必要一个安静的情况，以是我们可以conda一个情况，这个情况让我们更好做这次分析，因此对于SCENIC分析来说 ，由于r的分析过程很慢，而python的分析比力快，以是我们在linux要创建一个python的情况
conda创建一个python情况
conda create -n pyscenic python=3.8 #创建一个名为pyscenic的python3.8情况，这是一个难点，这里提示一下，假如安装failed的话，大概是.condarc文件配置的题目哈。conda info -e #查察全部的情况conda activate pyscenic #激活我们工作的情况，进入到该情况#在这个情况中配置一些依赖的东西conda install pipconda install -y numpy conda install -y -c anaconda cytoolzconda install -y scanpy #留意scanpy这个包，假如安装不了用pip安装：pip install -i https://mirrors.bfsu.edu.cn/pypi/web/simple scanpy #假如还是安装不了，就去Googlepip install pyscenic=0.11.1 #必须安装这个版本哈 由于在2022年5月份更新了，以是要安装从前的版本 否则后面运行会报错linux上分析流程
1.在服务器上创建一个工作文件夹，并往文件夹中传输4个文件，进入该文件夹，开始工作
mkdir pyscience #创建名为pyscience的工作文件夹，如许我们就在一个名叫pyscenic 的情况下创建了一些pyscience的工作文件夹。cd pyscience #进入该文件夹2.利用文件传输软件（Xftp）往我们pyscience文件夹中传送4个文件
- hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.feather
- hs_hgnc_tfs.txt
- motifs-v9-nr.hgnc-m0.001-o0.0.tbl
- scRNA.csv
表明这些文件
123这三个文件都是该上述那三个流程必要用到的参比数据集，假如分析不是人的就不是这三个哈，而且必须确保这三个数据集是齐备无损的
4这个文件还记得吗，是第一步在r语言操纵获得的单细胞表达矩阵
假如各人必要123文件可以简信我哈（佛系复兴）
3.run Linux流程
3.1 将scRNA.csv这个文件酿成一个loom文件
pip install ipython #安装ipython，便于运行python代码 #假如慢的话加其他镜像ipython #启动ipython 为了后面用python语言 #一行一行的输入以下代码，import os, sys os.getcwd()os.listdir(os.getcwd())import loompy as lpimport numpy as npimport scanpy as scx=sc.read_csv("scRNA.csv")row_attrs = {"Gene": np.array(x.var_names),}col_attrs = {"CellID": np.array(x.obs_names)}lp.create("sample.loom",x.X.transpose(),row_attrs,col_attrs)exit #退出ipython包当我们run完这些代码后，在我们工作的文件夹下会多一个sample.loom的文件，如许就转换成功了（这里跟jimmy大佬的代码有一些不一样，大佬直接用python脚本，而我怕堕落大概有一些依赖没有安装，以是用ipython包一行一行的run）
ipython运行示例图
记得肯定要exit退出ipython包
3.2 pyscenic 的3个步调之 GRN(linux里是用这个算法，R中是用GENIE3） 复制全部粘贴后运行即可 #这一步千万不要改num_workers 设置20就行
pyscenic grn \--num_workers 20 \--output adj.sample.tsv \--method grnboost2 \sample.loom \hs_hgnc_tfs.txt #转录因子文件，1839 个基因的名字列表3.3 pyscenic 的3个步调之 cistarget 复制全部粘贴后运行即可 #这一步改num_workers 设置99也行
pyscenic ctx \adj.sample.tsv \hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.feather \--annotations_fname motifs-v9-nr.hgnc-m0.001-o0.0.tbl \--expression_mtx_fname sample.loom \--mode "dask_multiprocessing" \--output reg.csv \--num_workers 3 \--mask_dropouts3.4 pyscenic 的3个步调之 AUCell 复制全部粘贴后运行即可 #这一步改num_workers 设置99也行
pyscenic aucell \sample.loom \reg.csv \--output sample_SCENIC.loom \--num_workers 3三个步调run完后，用ls查察一下发现，多了一个sample_SCENIC.loom文件，这就是我们SCENIC分析的结果，将其下载到本身电脑并导入R中，举行plot了，by the way， GRN和cistarget这两个步调很慢，发起看文献度过哈哈哈哈。
4.R中plot
#在这里我们要理解SCENIC做了什么分析#对每个细胞举行转录因子富集分析 筛选出较为真实的转录因子 并以细胞为单元 即每一个细胞对每一个TF举行打分（AUG) #以是我们可以按照细胞范例， 族，样原泉源，将相应的TFplot出来#sample_SCENIC.loom导入R里面举行可视化library(SCENIC)packageVersion("SCENIC")library(Seurat)library(SCopeLoomR)library(AUCell)library(SCENIC)library(dplyr)library(KernSmooth)library(RColorBrewer)library(plotly)library(BiocParallel)library(grid)library(ComplexHeatmap)library(data.table)library(scRNAseq)rm(list=ls())# 提取sample_SCENIC.loom 信息loom <- open_loom('sample_SCENIC.loom') #导入sample_SCENIC.loom文件#起首我们要把导入的loom处理惩罚成R中的数据#获取regulon regulon界说：TF与作用的genes#1.提取每一个TF与每一个gene作用系数regulons_incidMat <- get_regulons(loom, column.attr.name="Regulons") regulons_incidMat[1:4,1:4] #在这里就可以看出 每一个TF与每一个gene的作用数值regulons <- regulonsToGeneLists(regulons_incidMat) #做成一个list TF与其作用gene的list TF+genes 个人感觉这里假如后面想分析这个TF，则这里可以画这个TF与其作用的gene的网络图 #2.获得regulon的AUC 即TF在每一个细胞的激活程度regulonAUC <- get_regulons_AUC(loom,column.attr.name='RegulonsAUC')regulonAUC[1:4,1:4] #regulonAUC这个文件含有每一个TF在各个细胞中的表达量 列名为细胞名 行名为TF#3.找出在这单细胞数据中 高表达的TFregulonAucThresholds <- get_regulon_thresholds(loom) tail(regulonAucThresholds[order(as.numeric(names(regulonAucThresholds)))]) #这里可以看出哪一些TF是在这个单细胞数据中高表达的#4.这两步不知道是啥 embeddings <- get_embeddings(loom) #好像是什么嵌入 必须要做 否则后面会报错close_loom(loom)#如许就处理惩罚完成了#接下来我们就可以挑选本身感爱好的TF，举行个性化分析#流程：1.查察富集到的全部的TF 2.从这些TF中挑选我们本身感爱好的大概与课题相关的#1.查察富集到的全部的TF#必要每一个细胞的每一个TF的表达量#导入单细胞数据scRNA <- readRDS("scRNA.rds") #这个rds是我的seruat对象#将每一个细胞的每一个TF的表达量与单细胞的每一个细胞匹配sub_regulonAUC <- regulonAUC[,match(colnames(scRNA),colnames(regulonAUC))]dim(sub_regulonAUC)#确认是否划一identical(colnames(sub_regulonAUC), colnames(scRNA))#[1] TRUE 肯定要为TRUE，思索一下为啥，假如你用过r跑的话，r是要求你除了表达矩阵还要细胞的meta信息，而我们用linux跑的话，只用了细胞的表达矩阵，以是我们要判断我们我们做的分析是这个单细胞数据的分析，而且必须千篇一律，否则你后面将分析的结果添加到seruat对象里就是不match的#2.从这些TF中挑选我们本身感爱好的大概与课题相关的#挑选感爱好的TF（以是要符合 1.在这个数据中有表达 2.不在所有细胞中都高表达，否则无法做差别分析）names(regulons) #我们可以通过这个函数来看在这个单细胞数据中 有哪些TF表达 在此中挑选感爱好的TF#设置感爱好的TF regulonsToPlot = c("TWIST1(+)","NKX2-1(+)","ZNF831(+)","SIX4(+)","FOSB(+)","TBX21(+)")#将感爱好的TF的表达量到场到seruat对象中scRNA@meta.data = cbind(scRNA@meta.data,t(assay(sub_regulonAUC[regulonsToPlot,])))#3.可视化#设置画图的分组 用idents 也可以用其他的分组Idents(scRNA) <- sce$celltypetable(Idents(scRNA))DotPlot(scRNA, features = unique(regulonsToPlot)) + RotatedAxis()RidgePlot(scRNA, features = regulonsToPlot , ncol = 1)VlnPlot(scRNA, features = regulonsToPlot,pt.size = 0 ) FeaturePlot(scRNA, features = regulonsToPlot )更多的plot形式各人去Google哦，终于长舒一口气，这个包学了大概2个星期.........
写在末了，各人还是多去学习jimmy大佬的教程，真的很有资助，办理题目时面红耳赤，办理完后心情舒坦，这该死的结果感！！！
References:
https://cn.bing.com/search?q=pyscenic&form=QBLHCN&sp=-1&pq=pyscenic&sc=8-8&qs=n&sk=&cvid=875AE9ECDA3D4EE5BEE138BA629C438E
https://www.jianshu.com/p/0a29ecfaf21e
https://blog.csdn.net/qq_45688354/article/details/108014189 |
