Seurat软件学习11-细胞周期内容的分析

原创

小胡子刺猬的生信学习123

修改于 2023-03-28 20:13:58

8860

修改于 2023-03-28 20:13:58

文章被收录于专栏：文献分享及代码学习

Seurat软件学习1-多个模型得数据进行整合：https://cloud.tencent.com/developer/article/2130078

Seurat软件学习2-scrna数据整合分析：https://cloud.tencent.com/developer/article/2131431

Seurat软件学习3-scrna数据整合分析注释数据集：https://cloud.tencent.com/developer/article/2133583

Seurat软件学习4-使用RPCA进行快速整合数据集：https://cloud.tencent.com/developer/article/2134684

Seurat软件学习5-scRNA-Seq和scATAC-Seq数据整合：https://cloud.tencent.com/developer/article/2136814

Seurat软件学习6-多模型参考映射的方法：https://cloud.tencent.com/developer/article/2144475

Seurat软件学习7-同胞多组学结合方法-WNN：https://cloud.tencent.com/developer/article/2152008

Seurat软件学习8-不同细胞类型样本的分析流程：https://cloud.tencent.com/developer/article/2191271

Seurat软件学习9-sctransform的使用：https://cloud.tencent.com/developer/article/2208107?areaSource=&traceId=

Seurat软件学习10-SCTransformV2版本的运行方式：https://cloud.tencent.com/developer/article/2212131?areaSource=&traceId=

大家好，又回来安心写推文了，前段时间出差，加上课题需要补一下试验，由于自己是一直做数据的，没有试验的基础，忙着跟师兄师姐学习，导致自己的账号是一片荒芜。

今天细胞周期的计算，在很多的物种中，细胞周期的计算是一个很有意思的现象，比如在大豆根瘤细胞中，根瘤细胞的在不同的侵染期发生的内复制的现象也是人们很关注的事情；在人类的癌症细胞中，癌症细胞的快速增殖引起了癌变的发生。因此研究细胞的增殖周期和细胞周期的变化是十分重要的。

代码来源：https://satijalab.org/seurat/articles/cell_cycle_vignette.html

Cell-Cycle Scoring and Regression

我们展示了如何在scRNA-seq数据中减轻细胞周期异质性的影响，方法是根据典型的标记物计算细胞周期阶段的分数，并在预处理期间将这些分数从数据中回归。我们在小鼠造血祖细胞的数据集上证明了这一点（Nestorowa等，Blood 2016）。你可以在这里下载运行这个小流程所需的文件https://www.dropbox.com/s/3dby3bjsaf5arrw/cell_cycle_vignette_files.zip?dl=1。

library(Seurat)

# Read in the expression matrix The first row is a header row, the first column is rownames
exp.mat <- read.table(file = "../data/nestorawa_forcellcycle_expressionMatrix.txt", header = TRUE,
    as.is = TRUE, row.names = 1)

# A list of cell cycle markers, from Tirosh et al, 2015, is loaded with Seurat.  We can
# segregate this list into markers of G2/M phase and markers of S phase
s.genes <- cc.genes$s.genes
g2m.genes <- cc.genes$g2m.genes

# Create our Seurat object and complete the initalization steps
marrow <- CreateSeuratObject(counts = exp.mat)
marrow <- NormalizeData(marrow)
marrow <- FindVariableFeatures(marrow, selection.method = "vst")
marrow <- ScaleData(marrow, features = rownames(marrow))

我们对数据集进行PCA分析，选用FindVariableFeatures()找到高变的基因，可以发现大部分的pcs可以由一定的系谱进行解释，但是PC8和PC10分出来一类的细胞周期基因，包括TOP2A和MKI67。我们通过对后期的数据中这一类的数据进行vars处理，这样细胞周期出现的异质性不会对后面的pca分析产生很大的影响。

marrow <- RunPCA(marrow, features = VariableFeatures(marrow), ndims.print = 6:10, nfeatures.print = 10)
## PC_ 6 
## Positive:  LY6C2, AA467197, CYBB, MGST2, ITGB2, PF4, CD74, ATP1B1, GP1BB, TREM3 
## Negative:  SELL, ARL6IP1, CCL9, CD34, ADGRL4, BPIFC, NUSAP1, FAM64A, CD244, C030034L19RIK 
## PC_ 7 
## Positive:  HDC, CPA3, PGLYRP1, MS4A3, NKG7, UBE2C, CCNB1, NUSAP1, PLK1, FUT8 
## Negative:  F13A1, LY86, CFP, IRF8, CSF1R, TIFAB, IFI209, CCR2, TNS4, MS4A6C 
## PC_ 8 
## Positive:  WFDC17, SLC35D3, ADGRL4, VLDLR, CD33, H2AFY, P2RY14, IFI206, CCL9, CD34 
## Negative:  NUSAP1, UBE2C, KIF23, PLK1, CENPF, FAM64A, CCNB1, H2AFX, ID2, CDC20 
## PC_ 9 
## Positive:  IGKC, JCHAIN, LY6D, MZB1, CD74, IGLC2, FCRLA, IGKV4-50, IGHM, IGHV9-1 
## Negative:  SLC2A6, HBA-A1, HBA-A2, IGHV8-7, FCER1G, F13A1, HBB-BS, PLD4, HBB-BT, IGFBP4 
## PC_ 10 
## Positive:  H2AFX, FAM64A, ZFP383, NUSAP1, CDC25B, CENPF, GBP10, TOP2A, GBP6, GFRA1 
## Negative:  CTSW, XKRX, PRR5L, RORA, MBOAT4, A630014C17RIK, ZFP105, COL9A3, CLEC2I, TRAT1
DimHeatmap(marrow, dims = c(8, 10))

Assign Cell-Cycle Scores

首先，我们给每个细胞基于 G2/M 和 S phase markers的表达量进行评分。这些marker的表达水平应该是anticorrelated的，表达这marker的细胞很可能不是在循环和处于G1期。

我们用CellCycleScoring()函数进行分数计算，该函数在对象元数据中存储S和G2/M的分数，以及每个细胞在G2M、S或G1阶段的预测分类。CellCycleScoring()还可以通过传递set.ident = TRUE将Seurat对象的身份设置为细胞周期阶段（原始身份被存储为old.ident）。请注意，Seurat在下游的细胞周期回归中不使用离散的分类（G2M/G1/S）。相反，它使用G2M和S期的定量分数。然而，我们提供了我们预测的结果来进行后续的分析。

marrow <- CellCycleScoring(marrow, s.features = s.genes, g2m.features = g2m.genes, set.ident = TRUE)

# view cell cycle scores and phase assignments
head(marrow[[]])
##          orig.ident nCount_RNA nFeature_RNA     S.Score  G2M.Score Phase
## Prog_013       Prog    2563089        10211 -0.14248691 -0.4680395    G1
## Prog_019       Prog    3030620         9991 -0.16915786  0.5851766   G2M
## Prog_031       Prog    1293487        10192 -0.34627038 -0.3971879    G1
## Prog_037       Prog    1357987         9599 -0.44270212  0.6820229   G2M
## Prog_008       Prog    4079891        10540  0.55854051  0.1284359     S
## Prog_014       Prog    2569783        10788  0.07116218  0.3166073   G2M
##          old.ident
## Prog_013      Prog
## Prog_019      Prog
## Prog_031      Prog
## Prog_037      Prog
## Prog_008      Prog
## Prog_014      Prog
# Visualize the distribution of cell cycle markers across
RidgePlot(marrow, features = c("PCNA", "TOP2A", "MCM6", "MKI67"), ncol = 2)

# Running a PCA on cell cycle genes reveals, unsurprisingly, that cells separate entirely by
# phase
marrow <- RunPCA(marrow, features = c(s.genes, g2m.genes))
DimPlot(marrow)

我们根据Tirosh等人2016年描述的评分策略对单细胞进行评分。更多信息见Seurat中的AddModuleScore()，该函数可用于计算任何基因列表的监督模块分数。

Regress out cell cycle scores during data scaling

我们尝试在这个数据当中去除细胞周期长生的异质性。在Seurat v1.4中，主要通过RegressOut进行去除。然而，这个结果存放在缩放的结果中，覆盖了ScaleData() 的结果，因此我们现在把结果与ScaleData()进行合并。

对于每个基因，Seurat建立了基因表达与S和G2M细胞周期分数之间的关系。这个模型的标度残差代表了一个 "修正 "的表达矩阵，可以在下游用于降维分析。

marrow <- ScaleData(marrow, vars.to.regress = c("S.Score", "G2M.Score"), features = rownames(marrow))
# Now, a PCA on the variable genes no longer returns components associated with cell cycle
marrow <- RunPCA(marrow, features = VariableFeatures(marrow), nfeatures.print = 10)
## PC_ 1 
## Positive:  BLVRB, CAR2, KLF1, AQP1, CES2G, ERMAP, CAR1, FAM132A, RHD, SPHK1 
## Negative:  TMSB4X, H2AFY, CORO1A, PLAC8, EMB, MPO, PRTN3, CD34, LCP1, BC035044 
## PC_ 2 
## Positive:  ANGPT1, ADGRG1, MEIS1, ITGA2B, MPL, DAPP1, APOE, RAB37, GATA2, F2R 
## Negative:  LY6C2, ELANE, HP, IGSF6, ANXA3, CTSG, CLEC12A, TIFAB, SLPI, ALAS1 
## PC_ 3 
## Positive:  APOE, GATA2, NKG7, MUC13, MS4A3, RAB44, HDC, CPA3, FCGR3, TUBA8 
## Negative:  FLT3, DNTT, LSP1, WFDC17, MYL10, GIMAP6, LAX1, GPR171, TBXA2R, SATB1 
## PC_ 4 
## Positive:  CSRP3, ST8SIA6, DNTT, MPEG1, SCIN, LGALS1, CMAH, RGL1, APOE, MFSD2B 
## Negative:  PROCR, MPL, HLF, MMRN1, SERPINA3G, ESAM, GSTM1, D630039A03RIK, MYL10, LY6A 
## PC_ 5 
## Positive:  CPA3, LMO4, IKZF2, IFITM1, FUT8, MS4A2, SIGLECF, CSRP3, HDC, RAB44 
## Negative:  PF4, GP1BB, SDPR, F2RL2, RAB27B, SLC14A1, TREML1, PBX1, F2R, TUBA8
# When running a PCA on only cell cycle genes, cells no longer separate by cell-cycle phase
marrow <- RunPCA(marrow, features = c(s.genes, g2m.genes))
DimPlot(marrow)

可以发现在去除细胞周期的影响后，跟上面的pca降维图比，不能较好的区分细胞周期的现象。

由于最佳的细胞周期标志物在不同组织和物种中都有很好的保守性，我们发现这个程序在不同的数据集上都能稳定可靠地工作。

Alternate Workflow

上述程序删除了所有与细胞周期相关的内容。在某些情况下，我们发现这可能会对下游分析产生负面影响，特别是在分化过程中（如小鼠造血），干细胞处于静止状态，而分化细胞正在增殖（或相反）。在这种情况下，回归所有的细胞周期效应也会模糊干细胞和祖细胞之间的区别。

作为一个替代方案，我们建议回归出G2M和S期得分之间的差异。这意味着分离非周期性细胞和周期性细胞的信号将被保留，但增殖细胞之间的细胞周期阶段差异（通常是无趣的），将从数据中回归出来。

marrow$CC.Difference <- marrow$S.Score - marrow$G2M.Score
marrow <- ScaleData(marrow, vars.to.regress = "CC.Difference", features = rownames(marrow))
# cell cycle effects strongly mitigated in PCA
marrow <- RunPCA(marrow, features = VariableFeatures(marrow), nfeatures.print = 10)
## PC_ 1 
## Positive:  BLVRB, KLF1, ERMAP, FAM132A, CAR2, RHD, CES2G, SPHK1, AQP1, SLC38A5 
## Negative:  TMSB4X, CORO1A, PLAC8, H2AFY, LAPTM5, CD34, LCP1, TMEM176B, IGFBP4, EMB 
## PC_ 2 
## Positive:  APOE, GATA2, RAB37, ANGPT1, ADGRG1, MEIS1, MPL, F2R, PDZK1IP1, DAPP1 
## Negative:  CTSG, ELANE, LY6C2, HP, CLEC12A, ANXA3, IGSF6, TIFAB, SLPI, MPO 
## PC_ 3 
## Positive:  APOE, GATA2, NKG7, MUC13, ITGA2B, TUBA8, CPA3, RAB44, SLC18A2, CD9 
## Negative:  DNTT, FLT3, WFDC17, LSP1, MYL10, LAX1, GIMAP6, IGHM, CD24A, MN1 
## PC_ 4 
## Positive:  CSRP3, ST8SIA6, SCIN, LGALS1, APOE, ITGB7, MFSD2B, RGL1, DNTT, IGHV1-23 
## Negative:  MPL, MMRN1, PROCR, HLF, SERPINA3G, ESAM, PTGS1, D630039A03RIK, NDN, PPIC 
## PC_ 5 
## Positive:  HDC, LMO4, CSRP3, IFITM1, FCGR3, HLF, CPA3, PROCR, PGLYRP1, IKZF2 
## Negative:  GP1BB, PF4, SDPR, F2RL2, TREML1, RAB27B, SLC14A1, PBX1, PLEK, TUBA8
# when running a PCA on cell cycle genes, actively proliferating cells remain distinct from G1
# cells however, within actively proliferating cells, G2M and S phase cells group together
marrow <- RunPCA(marrow, features = c(s.genes, g2m.genes))
DimPlot(marrow)