第10章

章节实例：第10章 单细胞测序数据分析


单细胞测序技术的数据分析实例

（一）输入数据以及数据分析工具

1）输入数据

这里我们模拟了200个单细胞的30000个基因的表达值（以FPKM为单位），数据以矩阵的形式存储为SCS_RNA.csv文件。

2）选用数据分析工具

R软件平台：https://www.r-project.org/

R软件包-DESeq：http://www.bioconductor.org/packages/release/bioc/html/DESeq.html

R软件包–Statmod：https://cran.r-project.org/web/packages/statmod/

（二）数据的读入与归一化

1）数据在R软件中的读入与行名赋值（R控制台输入命令）：

## 读入csv数据

$ SCS_Count <- read.csc (“SCS_RNA.csv”, header=T)

## 为数据的行名赋值

$ rownames(SCS_Count) <- SCS_Count[ ,1]

$ SCS_Count <- SCS_Count[ ,2:ncol(SCS_Count)]

## 去除在所有样本中表达值均为0的基因行

$ SCS_Count <- SCS_Count[rowSums(SCS_Count)>0,]

2) 采用DESeq包进行数据的归一化（R控制台输入命令）：

## R平台载入DESeq包

$ require(DESeq)

## 利用DESeq包计算数据的library size

$ lib.size <- estimateSizeFactorsForMatrix(SCS_Count)

## 根据上述计算得到的library size对数据进行归一化

$ ed <- t(t(SCS_Count)/lib.size)

（三）根据归一化后的数据鉴定样本中高度差异表达的基因

## 计算方差估计值，变异系数

$ means <- rowMeans(ed)

$ vars <- apply(ed,1,var)

$ cv2 <- vars/means^2

## R平台载入Statmod包

$ require(Statmod)

## 回归曲线的拟合

$ minMeanForFit <- unname(quantile(means[which(cv2 > .3)], .95))

$ useForFit <- means >= minMeanForFit

$ fit <- glmgam.fit(cbind( a0 = 1, a1tilde = 1/means[useForFit]),cv2[useForFit])

$ a0 <- unname(fit$coefficients["a0"])

$ a1 <- unname(fit$coefficients["a1tilde"])

## 根据基因表达情况与拟合曲线的偏离程度对基因进行排序

$ afit <- a1/means+a0

$ varFitRatio <- vars/(afit*means^2)

$ varorder <- order(varFitRatio,decreasing=T)

$ oed <- ed[varorder,]

## 获取差异表达最显著的前100个基因并输出至文件

$rownames(oed)[1:100]

$ write.table(rownames(oed)[1:100], file="Top10_Variable_Genes")