作者,Evil Genius
今天我们需要实现的目标
看看方法部分
其实就是CNV分数矩阵做均一化,PCA降维,然后采用miloR进行聚类分析,分析CNV事件的相互联系。
我们来根据文献的内容实现一下,inferCNV大家自己做一下,我们从CNV矩阵开始。
libarary(Seurat)
library(miloR)
library(SingleCellExperiment)
library(scater)
library(scran)
library(dplyr)
library(patchwork)
expression_matrix <- read.csv("cnv_matrix.csv", row.names = 1)
seurat_object <- CreateSeuratObject(counts = expression_matrix)
seurat_object <- NormalizeData(seurat_object)
####这个地方看文献是均一化之后直接进行PCA分析
seurat_object <- RunPCA(seurat_object, features = rownames(seurat_object))
####获取PCA矩阵
pca_matrix = seurat_object@reductions$pca@cell.embeddings
pca_matrix
这个时候就要进行miloR分析
sce <- SingleCellExperiment(assays=list(counts=seurat_object@assays$RNA@counts, logcounts=seurat_object@assays$RNA@data),reducedDims=SimpleList(PCA=seurat_object@reductions$pca@cell.embeddings))
milo <- Milo(sce)
###Construct KNN graph
milo <- buildGraph(sce, k = 30, d = 30, reduced.dim = "PCA")
milo <- makeNhoods(milo)
####Defining representative neighbourhoods on the KNN graph
milo <- makeNhoods(milo, prop = 0.1, k = 30, d=30, refined = TRUE, reduced_dims = "PCA")
####Computing neighbourhood connectivity
milo <- calcNhoodDistance(milo, d=30, reduced.dim = "PCA")
####testing
#####注释信息
design <- data.frame(colData(embryo_milo))[,c("sample", "disease", "celltype")]
milo <- buildNhoodGraph(milo)
da_results <- testNhoods(milo, design = ~ sample + disease, design.df = design)
plotNhoodGroups(embryo_milo, da_results, layout="umap")
生活很好,有你更好