入门必读:https://www.jianshu.com/p/85e2e38df673
Q1:有的群2个或者3个不同细胞类型marker高表达,这时候如何给群定义具体哪种细胞呢?
随着单细胞技术的成熟,仅靠几个基因就判定一种细胞类型的模式可能要重新思考。在单细胞数据分析中有一个关键的步骤FindClusters(分群,以启发样本中可能有的细胞类型数量),但是这个目前用的方法是非监督聚类,也就是数据驱动的,不依赖生物学背景。而某一细胞类型的marker肯定是要依据生物学背景,您说的这种情况,就属于两者的矛盾的发生现场。有两个思路供您参考
- 一个群内部有不同细胞类型的marker表达,是不是这个群本身有异质性(即有亚群)。
- 找的marker是不是特异的,某些细胞类型是几个marker(有阳性有阴性)的。还有一种情况是分的不同的群有着相同的marker,这时可以把它们当做一个群。以上,基于FindClusters这个结果是符合生物学意义这一前提。
Q2:如果一个没有相似疾病报道的疾病做单细胞,细胞marker怎么确定?
可延伸为如何鉴定一个新的细胞类型(即,之前并没有人明确命名的细胞类型)。这个单纯靠单细胞转录组可能有些困难,可能要辅之以湿实验的验证(手动找marker)。单细胞已有的分析可以为确定marker的方法有哪些呢?第一个就是差异分析,检查marker的特异性;第二个是轨迹推断,检查该细胞类型与已知细胞类型之间的分化关系。鉴于单细胞技术之前细胞类型基本靠流式分选鉴定,建议做单细胞膜蛋白(传统是按细胞功能鉴定的)会比单纯的转录组好一些。随着单细胞技术的成熟,仅靠几个基因就判定一种细胞类型的模式可能要重新思考
Q3:怎么实现Seurat对象 与monocle bioconductor, scanpy,loom之间的转换
参考:单细胞转录组数据分析||Seurat3.1教程:Interoperability between single-cell object formats
关于Seurat2和Seurat3的转换:参考scRNA代码仓库,有官方说明书
https://www.jianshu.com/p/396345566479?from=singlemessage&tdsourcetag=s_pcqq_aiomsg
Q4:怎么做细胞与细胞互作
celltalker:单细胞转录组配受体分析: https://www.jianshu.com/p/2c2f94b4f072?tdsourcetag=s_pcqq_aiomsg
Q5:关于单细胞基因调控网络GRN
活性TFs及其靶基因的组合通常表示为GRNs
。探索GRNs是基因组研究领域的主要挑战之一。一旦确定了驱动并维持细胞状态行为的关键regulators,它们最终就可以用来做干扰这些调控过程的切入点。比如,结合一组特异的TFs组合,将成纤维细胞重编程为诱导性多能干细胞(iPS);许多其他的重编程途径也是通过特定的TFs组合驱动一个GRN来促使细胞状态的改变;最近在癌症治疗中进行了尝试,使癌细胞转成一种易于受特定药物侵害的状态。
- 使用聚类算法将细胞分为不同的细胞类型或状态
- 通过轨迹推断方法沿伪时间轴对细胞进行排序
从转录组学数据推断GRNs通常依赖于可以从表达模式中提取调控信息的假设。例如,具有相似行为的那些基因受共同机制(例如特异性的TF)的调节。依据这样的假设,调控网络(GRN)推断的目的可以是:
- 对导致细胞从一种状态转变成另一种状态的TF激活事件顺序进行建模
- 确定TFs潜在的靶标基因
- 鉴定一个细胞状态能够维持所依赖的特定主要regulators(或regulators组合)。
一类GRN推断方法着重于解密在动态过程中,细胞从一种状态转换到另一种状态所需的TF逻辑组合
这通常是通过布尔网络模型实现的。例如Single-Cell Network Synthesis (SCNS) toolkit和BoolTraineR,通过将每个细胞进行状态分类(基于TF表达)并连接有限数量差异的细胞来构建布尔网络。生成的状态图能够让我们找到参与细胞状态改变的关键TFs,并可用于预测TF的过表达或敲除后的效果。但是,这不涉及有关靶基因的信息。
另外,网络规模的增加会导致计算量的迅速增加。因此,这些工具只能模拟少数基因(<100)。所以通常选择相关TFs子集,再结合这些方法应用于动态过程的轨迹推断步骤。
调控网络推断的另一种方法是:将TF与候选靶基因连接,最终目的是确定驱动特定细胞状态的“主要regulators”
此类别中的一种主要方法是基于共表达分析并已广泛用于bulk基因表达数据,例如GENIE3和WGCNA。最近的研究已经成功地将相似的方法应用于单细胞数据。使用这些方法时需要考虑的要素包括以下假设:
regulator表达水平的变化直接影响下游靶标的表达
忽略了转录后调控
忽略了并非所有co-varying基因都一定是直接作用的靶标
这些方法对normalization和batch-effects敏感,可能会引入人工协变量
基于共表达算法的其中一部分算法专门针对动态过程中的单细胞转录组数据构建GRN模型。这些方法结合细胞沿时间轴的初始顺序(或预测的轨迹),同时对基因和调节子之间的表达动态进行建模,使用的技术包括:(non-linear) correlation、regression、covariance analysis 、multivariant information 、ordinary differential equation (ODE)模型以及其他的模型。
SCENIC workflow 包含3个主要步骤:
1.用GENIE3(随机森林) 或GRNBoost (Gradient Boosting梯度下降) 推断转录因子与候选靶基因之间的共表达模块。每个模块包含一个转录因子及其靶基因,纯粹基于共表达。2.使用RcisTarget分析每个共表达模块中的基因,以鉴定enriched motifs;仅保留TF motif富集的模块和targets,每个TF及其潜在的直接targets gene被称作一个调节子(regulon)
3.使用AUCell评估每个细胞中每个regulon的活性,AUCell分数用于生成Regulon活性矩阵,通过为每个regulon设置AUC阈值,可以将该矩阵进行二值化(0|1,on|off),这将确定Regulon在哪些细胞中处于“打开”状态。
使用RcisTarget是SCENIC不同于大多共表达算法的重要区别。由于GENIE3模块仅基于共表达,因此结果可能包含许多误报和间接target,为了鉴定推断的直接结合的靶标基因,使用RcisTarget对每个共表达模块进行顺式调控基序(motif)分析。仅保留具有正确基因上游调节子且显着富集TF motif的模块,并对它们进行修剪以除去缺乏基序支持的间接靶标,这些处理后的模块才称为regulon。
本文Ref参考很重要,作为进一步数据挖掘分析的基础
References
[1]
癌症治疗: https://www.sciencedirect.com/science/article/pii/S0959437X17300096
[2] 使用聚类算法将细胞分为不同的细胞类型或状态: https://www.sciencedirect.com/science/article/pii/S0098299717300493
[3] 通过轨迹推断方法沿伪时间轴对细胞进行排序: https://onlinelibrary.wiley.com/doi/full/10.1002/eji.201646347
[4] 1: https://www.sciencedirect.com/science/article/abs/pii/S0010482514000420
[5] 2: https://www.embopress.org/doi/full/10.1038/msb4100120
[6] 3: https://www.nature.com/nmeth/journal/v9/n8/abs/nmeth.2016.html
[7] SCNS: https://link.springer.com/chapter/10.1007/978-3-319-21690-4_38
[8] BoolTraineR: https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-016-1235-y
[9] 对iPS的重编程建模: https://www.sciencedirect.com/science/article/pii/S0092867412010215
[10] Moignard等人的工作: https://www.nature.com/articles/nbt.3154
[11] GENIE3: https://orbi.uliege.be/handle/2268/107580
[12] WGCNA: https://link.springer.com/article/10.1186/1471-2105-9-559
[13] SINCERA: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4658017/
[14] ACTION: https://www.nature.com/articles/s41467-018-03933-2
[15]SCENIC: https://www.nature.com/articles/nmeth.4463
Q6:怎么检测fusion? 我们的fusion转录本在bulk RNA-seq已经detect到,表达丰度也是不错的,那是不是我测序的深度足够就有可能捕获到fusion转录本呢?
测序深度是一个,10X的建库策略也是需要考虑的,不同转录本的fusion位置和长度不同,不好说这种可能性有多大
单细胞结构, SNV/CNV/snp以及fusion都是值得思考的。
10*毕竟不是全长测序,适不适合以及准确性需要考虑
STAR-solo:http://www.bio-info-trainee.com/5498.html?tdsourcetag=s_pcqq_aiomsg
Can we detect splice variants or isoforms using single cell 3’ RNA-seq data?:https://kb.10xgenomics.com/hc/en-us/articles/360021902671-Can-we-detect-splice-variants-or-isoforms-using-single-cell-3-RNA-seq-data-?tdsourcetag=s_pcqq_aiomsg
关于空间转录组
含源代码:https://github.com/SpatialTranscriptomicsResearch/st_pipeline大多是python写的
其他参考文章:https://www.jianshu.com/p/20df850ed6ca
Seurat新版教程:https://www.jianshu.com/p/f6da86489784?utm_campaign=haruki
公开课
如果按细胞大小10um来看,往上的是分辨率更高的亚显微层面(分辨率和基因检出数)有权衡
目前唯一的商业化空转产品是10x Visium,优势在于:可以获得无偏量的高通量基因表达,空间分辨率较高(55um,每个spot的大小(6.5x6.5mm,5000个spot,每个不同颜色的spot的barcode不同,一个spot可以有上百万条oligoUMI探针),检测面积大(42.5um,像ISH相关只能检测1-几个mm^2),检测效率(灵敏度)高,完整、简易的操作流程,可和现有实验方法整合
组织优化过程很重要,组织通透化时间决定了测序检测基因的灵敏度,组织优化过程中cDNA添加荧光标签,显微拍照观察,通透时间决定了mRNA的弥散程度,可以设置时间梯度选择优化条件
测序深度与样本有关,胚胎、大脑等复杂样本建议测深,肺脏等其他可以浅一些,每个spot 50k个reads pair
https://satijalab.org/seurat/v3.1/spatial_vignette.html官网上除了空间转录组还包含了其他各种vigenettes比如多组学数据整合的guide
组织样本的厚度一般10um-20um
Practical
参考:https://satijalab.org/seurat/v3.1/spatial_vignette.html
这个地方我在下载的时候报错了:除非运气很好不然大概率会报错的,最好是直接找上Github:https://github.com/satijalab/seurat-data 这个专门的SeuratData包可以当做是一个meta package,管理Seurat支持的数据集
devtools::install_github("satijalab/seurat", ref = "spatial")
查看可下载的数据集然后下载,如果还不行就手动把url对应的包下载到lib,下载好解压再重新载入数据就成功了(可以看到载入以后显示出已下载的数据):
> library(SeuratData)
-- Installed datasets ------------------------------------- SeuratData v0.2.1 --
√ stxBrain 0.1.1
-------------------------------------- Key -------------------------------------
√ Dataset loaded successfully
> Dataset built with a newer version of Seurat than installed
(?) Unknown version of Seurat installed
注意Visium是以spot为基本单位的,meta.data的每一行为barcode,记录了每一个barcode对应的UMI count和测得的基因数(一个spot测几千个基因)
> head(brain@meta.data)
orig.ident nCount_Spatial nFeature_Spatial slice region
AAACAAGTATCTCCCA-1 anterior1 13069 4242 1 anterior
AAACACCAATAACTGC-1 anterior1 37448 7860 1 anterior
AAACAGAGCGACTCCT-1 anterior1 28475 6332 1 anterior
AAACAGCTTTCAGAAG-1 anterior1 39718 7957 1 anterior
AAACAGGGTCTATATT-1 anterior1 33392 7791 1 anterior
AAACATGGTGAGAGGA-1 anterior1 20955 6291 1 anterior
完成SCT和cluster后meta.data中会多出几列:
brain@meta.data$seurat_clusters
brain@meta.data$SCT_snn_res.0.8
brain@meta.data$nFeature_SCT
brain@meta.data$nCount_SCT
但是我发现因为我并没能下载他们pre-released版本,所以没有需要的函数,然后在github上找解决方案:
https://github.com/search?q=seurat+v3.2&type=Issues不过没什么用,所以我唯一能想到的就是改天再做
#github网站居然总是出现安全链接失败。只能多试几次,看什么时候人品回复正常
> devtools::install_github("satijalab/seurat", ref = "spatial")
Error: Failed to install 'Seurat' from GitHub:
schannel: failed to receive handshake, SSL/TLS connection failed
这里需要注意一个问题:
install_github()安装包出错: https://www.jianshu.com/p/58e33b71b466
但是没有解决这个问题,于是我尝试了本地安装,这样可以绕开github安装
[https://codeload.github.com/satijalab/seurat/legacy.tar.gz/spatial](https://codeload.github.com/satijalab/seurat/legacy.tar.gz/spatial) # 下载地址
install.packages("H:/singlecell/Seurat/satijalab-seurat-v3.1.1-302-g1cb8a3d.tar.gz", repos = NULL, type = "source")
#需要编译,可能会报错,Rcpp系列的包该装的都装上
#成功
packageVersion("Seurat")
[1] ‘3.1.4.9902’#这个pre-release版本编号比较奇特
「R包」如何实现Seurat 2 和 Seurat3 在同一个环境下共存:https://www.baidu.com/link?url=aw0bEQi-VjJP2wC7m3jEAJwsOhqUoeHyf1M2wPr5p6iLTMQj1SeMepP7e6-wA3Y83Of_GVXvkqZ0cB-p7wXFSYOaTbpE-YjhFrGyXYSBNKe&wd=&eqid=bd40068700049ff0000000065e706341
下面这种方法是在没有预注释的情况下根据基因表达的情况,找出那些具有强烈空间定位倾向的基因:
#在FindSpatiallyVariables中
#实现的另一种方法是在没有预注释的情况下搜索显示空间模式的特性。
#默认的方法(method = 'markvariogram ')受到 Trendsceek,的启发,
#后者将空间转录组学数据建模为标记点过程,并计算一个' variogram ',
#它识别其表达水平取决于其空间位置的基因。更具体地说,
#这个过程计算伽玛(r)值,测量两个点之间一定的“r”距离的相关性。
#默认情况下,我们在这些分析中使用的r值为‘5’
#注意这种计算会耗时较长
brain <- FindSpatiallyVariableFeatures(brain, assay = "SCT", features = VariableFeatures(brain)[1:1000],
selection.method = "markvariogram")
对于brain数据集我们可以提取我们感兴趣的区域做进一步分析
> head(brain@images$anterior1@coordinates)
tissue row col imagerow imagecol
AAACAAGTATCTCCCA-1 1 50 102 7475 8501
AAACACCAATAACTGC-1 1 59 19 8553 2788
AAACAGAGCGACTCCT-1 1 14 94 3164 7950
AAACAGCTTTCAGAAG-1 1 43 9 6637 2099
AAACAGGGTCTATATT-1 1 47 13 7116 2375
AAACATGGTGAGAGGA-1 1 62 0 8913 1480
同样在空间转录组中也可以实现数据集的整合,但是需要一个reference dataset: a reference scRNA-seq dataset of ~14,000 adult mouse cortical cell taxonomy from the Allen Institute, generated with the SMART-Seq2 protocol
# note that setting ncells=3000 normalizes the full dataset but learns noise models on 3k cells
# this speeds up SCTransform dramatically with no loss in performance
library(dplyr)
allen_reference <- SCTransform(allen_reference, ncells = 3000, verbose = FALSE) %>% RunPCA(verbose = FALSE) %>%
RunUMAP(dims = 1:30)
# After subsetting, we renormalize cortex 注意,reference和query data都要进行归一化和cluster操作
cortex <- SCTransform(cortex, assay = "Spatial", verbose = FALSE) %>% RunPCA(verbose = FALSE)
# the annotation is stored in the 'subclass' column of object metadata
DimPlot(allen_reference, group.by = "subclass", label = TRUE)
anchors <- FindTransferAnchors(reference = allen_reference, query = cortex, normalization.method = "SCT")
predictions.assay <- TransferData(anchorset = anchors, refdata = allen_reference$subclass, prediction.assay = TRUE,
weight.reduction = cortex[["pca"]])
cortex[["predictions"]] <- predictions.assay
获得anchors以后我们就可以进行数据的整合工作:
>anchors
An AnchorSet object containing 1236 anchors between 1 Seurat objects
This can be used as input to IntegrateData or TransferData.
上述步骤获得了anchors并进行了注释的映射(TransferData),这样对于cortex的每一个spot我们都有了subclass注释的预测性得分。
可以取出不同layer的细胞:
DefaultAssay(cortex) <- "predictions"
SpatialFeaturePlot(cortex, features = c("L2/3 IT", "L4"), pt.size.factor = 1.6, ncol = 2, crop = TRUE)
基于prediction的得分,同样能找出那些有强烈空间定位倾向的细胞类型,不过和之前找variable features不同,这次用的是prediction score作为 marker而不是基因表达:
cortex <- FindSpatiallyVariableFeatures(cortex, assay = "predictions", selection.method = "markvariogram",
features = rownames(cortex), r.metric = 5, slot = "data")
top.clusters <- head(SpatiallyVariableFeatures(cortex), 4)
SpatialPlot(object = cortex, features = top.clusters, ncol = 2)
同一个Seurat对象还可以包含不同的slice,比如说,我们取另一区域的brain数据集并归一化:
brain2 <- LoadData("stxBrain", type = "posterior1")
brain2 <- SCTransform(brain2, assay = "Spatial", verbose = FALSE)
然后和原来的合并成一个Seurat对象:
brain.merge <- merge(brain, brain2)
接下来就可以对多个slice同时做PCA降维、聚类等分析了
DefaultAssay(brain.merge) <- "SCT"
VariableFeatures(brain.merge) <- c(VariableFeatures(brain), VariableFeatures(brain2))
brain.merge <- RunPCA(brain.merge, verbose = FALSE)
brain.merge <- FindNeighbors(brain.merge, dims = 1:30)
brain.merge <- FindClusters(brain.merge, verbose = FALSE)
brain.merge <- RunUMAP(brain.merge, dims = 1:30)
#其他的图展示:
SpatialDimPlot(brain.merge)
SpatialFeaturePlot(brain.merge, features = c("Hpca", "Plp1"))
单细胞拟时分析
https://www.jianshu.com/p/e94cff521ebc?tdsourcetag=s_pcqq_aiomsg
UMAP与t-NSE的区别
UMAP 与 t-SNE 均是非线性降维,多用于数据可视化
UMAP 结构与t-SNE一致
UMAP 计算更快
UMAP 能更好地反映高纬结构,比t-SNE有着更好的连续性
这种连续性反映到单细胞分析中就是能更好滴可视化细胞的分化状态(UMAP better represents the multi-branched continuous trajectory of hematopoietic development)
参考:https://www.jianshu.com/p/5aa1e2467339
We conclude that both techniques are similar in their visualisation capabilities, but UMAP has a clear advantage over t-SNE in runtime, making it highly plausible to employ UMAP as an alternative to t-SNE in mIF data analysis.
反卷积法推断bulk RNA数据中的细胞组成
整体思路:细胞相互作用的复杂网络支配免疫系统与肿瘤细胞之间的互动,了解实体瘤的特定免疫细胞组成,对于预测患者对免疫疗法有何反应显得至关重要。在这篇文章中,作者使用肿瘤单细胞RNA测序数据的适应症特异性和细胞类型特异性参考基因表达谱(RGEP),深入分析如何通过数学反卷积从bulk数据中得出实体瘤的细胞组成数据。证明了肿瘤衍生的RGEP对成功的去卷积至关重要,而来自外周血的RGEP则不足。我们区分了9种主要细胞类型以及3种T细胞亚型。使用源自肿瘤的RGEP,我们可以估计许多与肿瘤相关免疫细胞和基质细胞的含量,治疗相关的比例,以及完善的恶性细胞基因表达谱
https://www.jianshu.com/p/74e8e4fa4643
从理论上讲,如果可以为每个肿瘤相关细胞建立参考基因表达谱(RGEP),则可以从其整体基因表达谱推断出实体瘤的免疫,肿瘤和基质细胞含量。从数学上讲,这类反推问题称为反卷积
我们发现,来自不同患者的适应症特异性免疫细胞RNA-seq谱图彼此足够相似,可以为每种细胞类型定义一个共有谱图,并且这些共有谱图可以对肿瘤bulk谱图进行准确的反卷积。
还可以参考:对混合细胞类型的转录组数据去卷积 (综述)https://www.jianshu.com/p/87caf98f3ca0
Gene expression distribution deconvolution in single-cell RNA sequencing:https://www.pnas.org/content/115/28/E6437
单细胞转录组高级分析介绍
10xGenomics公司提供了一种将转录组与V(D)J联合分析的平台,我们可以对同一个细胞,同时获得其转录组数据和V(D)J数据,从而进一步深入研究其免疫特征。
VDJ结果可以反映细胞亚群的相似性,如果两个亚群的VDJ表型一致,那么两个亚群的相似性越高。
单细胞RNA-seq、ATAC-seq分析工具集锦
https://github.com/seandavi/awesome-single-cell
Cell丨首个哺乳动物单细胞染色质调控图谱——颉伟点评
http://www.360doc.com/content/18/0807/19/28251648_776437557.shtml