Topic 1 Single_Cell_analysis（1）

Part1 Plan

Part2 原始数据复现——LUAD

Part3 Scissor——Cox regression

scissor的输入是三个数据源：单细胞表达矩阵，批量表达矩阵和interested的表型。每个批量的样本表型注释可以是连续变量、二元组向量或者临床生存数据。

1.prework

• miniconda安装：略，blog上有详细教程
• 在miniconda 安装seurat：seurat-package用于下载LUAD的singcell数据

2.准备单细胞数据

location <- "https://xialab.s3-us-west-2.amazonaws.com/Duanchen/Scissor_data/"
load(url(paste0(location, 'scRNA-seq.RData')))
dim(sc_dataset)

说明共有33694个基因和4102个cell 。对于Scissor中使用的单细胞数据，首选是包含预处理数据和构建网络的Seurat，使用Seurat-package预处理这些数据，构建一个cell-cell之间的similarity网络
注：Seurat_preprocessing()是author自己定义的function，contained in Scissor-package

2.1 Seurat_preprocessing（）

(1)realize functions

• 对于每行数据创建一个Seurat 对象（即对于每个基因创建一个Seurat对象）
• 标准化需要处理的数据
• 特征识别
• 确定数据尺度和中新功能
• 对于sc_dataset构造一个SNN（shared nearest neighbor）图
• 通过基于SNN模块化优化的聚类算法识别细胞集群
• 降维选择特征数据（PCA，t-SNE，UMAP）

sc_dataset <- Seurat_preprocessing(sc_dataset, verbose = F)
class(sc_dataset)

sc_dataset

names(sc_dataset)

分簇——umap

DimPlot(sc_dataset, reduction = 'umap', label = T, label.size = 10)

分簇——pca

DimPlot(sc_dataset,reduction='pca',label=T,label.size=10)

分簇——t-sne

DimPlot(sc_dataset,reduction='tsne',label=T,label.size=10)

(2) mainly parameters

normalization methods :

• LogNormalize：每个cell的特征数除以该cell的所有特征数，并乘以Scale.Factor——使用log1p转换为自然数
• CLR：应用于居中的log比转换。
• RC：relative count.每个cell的特征数除以该cell的所有特征数，并乘以Scale.Factor
  注：scale.factor是cell-level的归一化比例，可以设置

selectin.method：——变量特征选择的方法

• vst：首先根据local的loss进行多项式回归。挑选符合log(方差)和log(平均值)的关系，然后使用观察到的方差和均值标准化特征值。去除最大值后在标准化值后计算特征变量
• mvp：首先函数计算每个特征的平均表达和dispersion，之后依据平均表达将特征划分进num.bin的bins中，并计算每个bin中dispersion的z-score(这样做的目的是识别不同的特征在编译和表达之间的关系强弱)
• dispersion：选择具有高dispersion值得基因

3.准备bulk-cell和表型phenotype数据

load(url(paste0(location, 'TCGA_LUAD_exp1.RData')))
load(url(paste0(location, 'TCGA_LUAD_survival.RData')))

判断bulk-cell中的每个细胞是否都存在于表型数据中的细胞

all(colnames(bulk_dataset) == bulk_survival$TCGA_patient_barcode)

4.执行scissor选择重要的带有信号的细胞

infos1 <- Scissor(bulk_dataset, sc_dataset, phenotype, alpha = 0.05, 
                 family = "cox", Save_file = 'Scissor_LUAD_survival.RData')

scissor首先输出五个single-cell和bulk-sample之间相关性的summary。这些值都是正值，并且这些值都不接近于0。

注：在使用过程中，如果使用数据的中间相关值小于0.01，将会给出warning，因为这会导致细胞和表型之间的关联度不可靠

4.1 查看scissor后的数据

4.2 使用Scissorh后筛选出的positive和negative基因分别存在Scissor_pos和Scissor_neg中

将positive和negative分别连同整个sc_dataset绘制在一张map上

Scissor_select <- rep(0, ncol(sc_dataset))
names(Scissor_select) <- colnames(sc_dataset)
Scissor_select[infos1$Scissor_pos] <- 1
Scissor_select[infos1$Scissor_neg] <- 2
sc_dataset <- AddMetaData(sc_dataset, metadata = Scissor_select, col.name = "scissor")
DimPlot(sc_dataset, reduction = 'umap', group.by = 'scissor', cols = c('grey','indianred1','royalblue'), pt.size = 1.2, order = c(2,1))

其中：1 代表positive，2代表negative

4.3 model建立所需要的parameter

5. 调整模型参数

在上述过程中得到了新的参数，因此将参数设置为alpha=0.05——参数alpha平衡了L1范数和基于网络的惩罚影响。较大的alpha倾向于强调L1范数以鼓励稀疏性，使得在Scissor选择的cell数少于使用较小alpha的结果。在实际的应用中，Scissor可以自动将回归输入保存到RData中，方便用户调整不同的alpha。可以设置Load_file='Scissor_LUAD_survival.RData’，并使用默认的alpha=NULL运行Scissor：（这时会根据cutoff的不同将回归得到的值不断返回model，重新构建新的model，直到所选择的cell占比达到一定值的时候停止）

infos2 <- Scissor(bulk_dataset, sc_dataset, phenotype, alpha = NULL, cutoff = 0.03, family = "cox", Load_file = 'Scissor_LUAD_survival.RData')

附 存储下载数据并转化

将数据存在本地，下次直接使用。存的时候直接用write.table(file,"D:/Topic/data/sicssor_data/")，看了一下每个三个矩阵的class
sc_dataset：原来的应该也是matrix吧，我已经转成Seurat了，下次再看吧
bulk_dataset：matrix
bulk_survival：data.frame
都存为table了，下次直接转：

• 将table转为matrix：bulk_dataset<-as.matrix(bulk_dataset)
• 将table转为data.frame：bulk_survival<-as.data.frame.array(bulk_survival)