Monocle操作笔记-2：数据读取和预处理

monocle操作数据格式CellDataSet (CDs)，相较于seurat的seurat object，该格式数据更为精简

1. 数据输入输出

首先需要利用3个文件生成CellDataSet格式单细胞表达数据

(1) exprs，基因表达数据，行-基因，列-细胞

monocle支持三种类型的表达数据，UMI counts（推荐）、 TPM/FPKM、log-transformed FPKM/TPMs，需要注意的是在下游分析中monocle会自动进行normalization，所以不要数据normalized data，或提前normalize data。需要针对不同类型的表达数据选择合适的分布negbinomial.size()-UMI，negbinomial()，tobit()，gaussianff()。

(2) phenoData，描述细胞的各种元属性，相当于seurat_obj@meta.data，行-细胞，列-属性

(3) featureData，基因名，行-基因名

(4) 支持直接从seurat或scater读取（输出）数据，但是非常垃圾，试了，根本用不了，还是老老实实手动生成数据吧。

#从seurat object生成CDS
example  ##示例数据为seurat object，已经进行过注释

pd <- new("AnnotatedDataFrame", data = example@meta.data) ##pd, phenoData

gene_annotation=data.frame(gene_short_name = rownames(example[[RNA]]),
                           stringsAsFactors=F)
rownames(gene_annotation)<-gene_annotation$gene_short_name
fd <- new("AnnotatedDataFrame", data = gene_annotation)   ##fd, featureData

HSMM <- newCellDataSet(GetAssayData(example,slot="counts",assay=assay),
                       phenoData = pd,
                       featureData = fd,
                       expressionFamily=negbinomial.size())

如果所分析的数据没有UMI count，而是FPKM/TPM值，这时就需要通过一种叫Census 的方法将这些基因表达值转换为RPC（mRNAs per cell），可极大的提高分析的可行性和结果的准确性。

我没有这种数据，下面是官网代码

pd <- new("AnnotatedDataFrame", data = HSMM_sample_sheet)
fd <- new("AnnotatedDataFrame", data = HSMM_gene_annotation)

# First create a CellDataSet from the relative expression levels
HSMM <- newCellDataSet(as.matrix(HSMM_expr_matrix),
                phenoData = pd,
                featureData = fd,
                lowerDetectionLimit = 0.1,
                expressionFamily = tobit(Lower = 0.1))   ###因为是FPKM/TPM值，所以分布模型是tobit()

# Next, use it to estimate RNA counts
rpc_matrix <- relative2abs(HSMM, method = "num_genes")

# Now, make a new CellDataSet using the RNA counts
HSMM <- newCellDataSet(as(as.matrix(rpc_matrix), "sparseMatrix"),
                phenoData = pd,
                featureData = fd,
                lowerDetectionLimit = 0.5,
                expressionFamily = negbinomial.size())   ###转为了RPS，分布模型改为negbinomial.size()

2. Estimate size factors and dispersions

必须要进行的两步操作，应该是为了评估数据大小和离散程度。命令非常简单，没有太多参数。

size factor与下游分析的数据标准化有关

dispersion与下有分析的差异分析有关

HSMM <- estimateSizeFactors(HSMM)
HSMM <- estimateDispersions(HSMM)