sup-patient1-PBMC.Rmd

---
title: "Sup-patient1-PBMC"
author: "jmzeng1314@163.com" 
date: "`r format(Sys.time(), '%d %B, %Y')`"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
knitr::opts_chunk$set(warning = F)
knitr::opts_chunk$set(message  = F) 
```

## 引言

这里直接读取作者给定的第一个病人的Gene expression analysis: discovery patient PBMC，用的是 10x genomics 3' Chromium expression assay

Following sequence alignment and filtering, a total of 12,874 cells were analyzed. 

最后是 17,712 genes and 12,874 cells, 所以对计算机的考验很大，而且文章使用的是seurat2，所以这里我们选取并且安装seurat2

## 载入必要的R包

需要自行下载安装一些必要的R包！ 而且需要注意版本 Seurat 

因为大量学员在中国大陆，通常不建议大家使用下面的R包安装方法，建议是切换镜像后再下载R包。

参考：http://www.bio-info-trainee.com/3727.html


```{r,eval=FALSE}
# 下面代码不运行。
# Enter commands in R (or R studio, if installed)
# Install the devtools package from Hadley Wickham
install.packages('devtools')
# Replace '2.3.0' with your desired version
devtools::install_version(package = 'Seurat', version = package_version('2.3.0'))

library(Seurat)
```

加载R包

```{r}
rm(list = ls()) # clear the environment
#load all the necessary libraries
options(warn=-1) # turn off warning message globally
suppressMessages(library(Seurat))

```

## 读入文章关于第一个病人的PBMC表达矩阵

```{r}
start_time <- Sys.time()
# 如果觉得这里较慢，可以使用 data.table 包的 fread函数。
raw_dataPBMC <- read.csv('../Output_2018-03-12/GSE117988_raw.expMatrix_PBMC.csv.gz', header = TRUE, row.names = 1)
end_time <- Sys.time()
end_time - start_time
# 通常电脑一分钟可以搞定。

dim(raw_dataPBMC) 

start_time <- Sys.time()
# 按照列，对每一个细胞进行内部归一化，主要是统一文库大小。
dataPBMC <- log2(1 + sweep(raw_dataPBMC, 2, median(colSums(raw_dataPBMC))/colSums(raw_dataPBMC), '*')) # Normalization
end_time <- Sys.time()
end_time - start_time

# 下面是简单的字符串切割,  ExtractField {Seurat}
head(colnames(dataPBMC))
timePoints <- sapply(colnames(dataPBMC), function(x) ExtractField(x, 2, '[.]'))
timePoints <-ifelse(timePoints == '1', 'PBMC_Pre', 
                    ifelse(timePoints == '2', 'PBMC_EarlyD27',
                           ifelse(timePoints == '3', 'PBMC_RespD376', 'PBMC_ARD614')))
table(timePoints)
# 可以看到是治疗前，加上治疗中的3个时间点，这4个分组信息。
```

## 表达矩阵的质量控制

简单看看表达矩阵的性质，主要是基因数量，细胞数量;

以及每个细胞表达基因的数量，和每个基因在多少个细胞里面表达。

```{r}
dim(dataPBMC)

# 可以看到，近2万的基因里面，
# 绝大部分基因只在一万多细胞的200个不到的表达
fivenum(apply(dataPBMC,1,function(x) sum(x>0) ))
boxplot(apply(dataPBMC,1,function(x) sum(x>0) ))

# 可以看到，一万多细胞里面
# 绝大部分细胞只能检测不到500个基因。
fivenum(apply(dataPBMC,2,function(x) sum(x>0) ))
hist(apply(dataPBMC,2,function(x) sum(x>0) ))
```


## 然后创建Seurat的对象

```{r}
start_time <- Sys.time()
# Create Seurat object
PBMC <- CreateSeuratObject(raw.data = dataPBMC, 
                           min.cells = 1, min.genes = 0, project = '10x_PBMC') # already normalized
dim(dataPBMC)
PBMC # 17,712 genes and 12,874 cells
# 可以看到上面创建Seurat对象的那些参数并没有过滤基因或者细胞。
end_time <- Sys.time()
end_time - start_time

# Add meta.data (nUMI and timePoints)
PBMC <- AddMetaData(object = PBMC, 
                    metadata = apply(raw_dataPBMC, 2, sum),
                    col.name = 'nUMI_raw')
PBMC <- AddMetaData(object = PBMC, metadata = timePoints, col.name = 'TimePoints')

```

## 一些质控

这里绘图，可以指定分组，前提是这个分组变量存在于meta信息里面，我们创建对象后使用函数添加了 TimePoints 属性，所以可以用来进行可视化。

这里是：'TimePoints'

```{r}
sce=PBMC
VlnPlot(object = sce, 
        features.plot = c("nGene", "nUMI"), 
        group.by = 'TimePoints', nCol = 2)
GenePlot(object = sce, gene1 = "nUMI", gene2 = "nGene")

```

可以看看高表达量基因是哪些

```{r}
tail(sort(Matrix::rowSums(sce@raw.data)))
## 散点图可视化任意两个基因的一些属性（通常是细胞的度量）
# 这里选取两个基因。
tmp=names(sort(Matrix::rowSums(sce@raw.data),decreasing = T))
GenePlot(object = sce, gene1 = tmp[1], gene2 = tmp[2])

# 散点图可视化任意两个细胞的一些属性（通常是基因的度量）
# 这里选取两个细胞
CellPlot(sce,sce@cell.names[3],sce@cell.names[4],do.ident = FALSE)

```

## 最后标准聚类可视化

很简单的流程，先ScaleData，再FindVariableGenes，然后根据找到的高变异基因进行RunPCA，再根据PCA结果进行FindClusters即可，最后再RunTSNE后进行可视化。

```{r}
start_time <- Sys.time()
# 最耗费时间的步骤在这里。
PBMC <- ScaleData(object = PBMC, vars.to.regress = c('nUMI_raw'), model.use = 'linear', use.umi = FALSE)
end_time <- Sys.time()
end_time - start_time

start_time <- Sys.time()
PBMC <- FindVariableGenes(object = PBMC, 
                          mean.function = ExpMean,
                          dispersion.function = LogVMR,
                          x.low.cutoff = 0.0125,
                          x.high.cutoff = 3, 
                          y.cutoff = 0.5)
head(PBMC@var.genes)
length(PBMC@var.genes)

PBMC <- RunPCA(object = PBMC, pc.genes = PBMC@var.genes)
## 避免太多log日志被打印出来。
PBMC <- FindClusters(object = PBMC, 
                     reduction.type = "pca", 
                     dims.use = 1:10, 
                     resolution = 1, 
                     print.output = 0,
                     k.param = 35, save.SNN = TRUE) # 13 clusters
PBMC <- RunTSNE(object = PBMC, dims.use = 1:10)
# 配色这里直接使用文章配色方案。 
TSNEPlot(PBMC, colors.use = c('green4', 'pink', '#FF7F00', 'orchid', '#99c9fb', 'dodgerblue2', 'grey30', 'yellow', 'grey60', 'grey', 'red', '#FB9A99', 'black'))
end_time <- Sys.time()
end_time - start_time

```

## 输出seurat结果后面使用

```{r}

start_time <- Sys.time()
save(PBMC,file = 'patient1.PBMC.output.Rdata')
end_time <- Sys.time()
end_time - start_time
# 这个步骤会输出文件 1.75G 
```


最后，这 13 clusters要进行注释，才能发表，如下所示：
![](../figures/patient1-PMBC-tSNE-13-cluster.png)


作者文章里面是Representative marker genes shown in Supplementary Fig. 7. 如下所示

![](../figures/sup-figure-7.png)

可以看到作者对PBMC里面的细胞都挑选了一个基因就命名了。

## 显示运行环境

```{r}
sessionInfo()
```