cellranger count 计算的结果只能作为错略观测的结果，如果需要进一步分析聚类细胞，还需要进行下游分析，这里使用官方推荐 R 包（Seurat 3.0）
流程参考官方外周血分析标准流程（ https://satijalab.org/seurat/v3.0/pbmc3k_tutorial.html ）
Rstudio操作过程：
## 安装seurat
install.packages('Seurat')
## 载入seurat包
library(dplyr)
library(Seurat)
## 读入pbmc数据（文件夹路径不能包含中文，注意“/“的方向不能错误，这里读取的是10x处理的文件，也可以处理其它矩阵文件，具体怎样操作现在还不知道，文件夹中的3个文件分别是：barcodes.tsv，genes.tsv，matrix.mtx，文件的名字不能错，否则读取不到）
pbmc.data <- Read10X(data.dir = "D:/pbmc3k_filtered_gene_bc_matrices/filtered_gene_bc_matrices/hg19/")
## 查看稀疏矩阵的维度,即基因数和细胞数
dim(pbmc.data)
pbmc.data[1:10,1:6]
## 创建Seurat对象与数据过滤，除去一些质量差的细胞（这里读取的是单细胞 count 结果中的矩阵目录；在对象生成的过程中，做了初步的过滤；留下所有在>=3 个细胞中表达的基因 min.cells = 3；留下所有检测到>=200 个基因的细胞 min.genes = 200。)
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)
pbmc
##计算每个细胞的线粒体基因转录本数的百分比（%）,使用[[ ]] 操作符存放到metadata中，mit-开头的为线粒体基因
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
##展示基因及线粒体百分比（这里将其进行标记并统计其分布频率，"nFeature_RNA"为基因数，"nCount_RNA"为细胞数，"percent.mt"为线粒体占比）
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
CombinePlots(plots = list(plot1, plot2))
## 过滤细胞：根据上面小提琴图中基因数"nFeature_RNA"和线粒体数"percent.mt"，分别设置过滤参数，这里基因数 200-2500，线粒体百分比为小于 5%，保留gene数大于200小于2500的细胞；目的是去掉空GEMs和1个GEMs包含2个以上细胞的数据；而保留线粒体基因的转录本数低于5%的细胞,为了过滤掉死细胞等低质量的细胞数据。
pbmc <- subset(pbmc, subset = nFeature_RNA >200 &nFeature_RNA <2500 &percent.mt <5)
## 表达量数据标准化,LogNormalize的算法：A = log( 1 + ( UMIA ÷ UMITotal ) × 10000
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
#pbmc <- NormalizeData(pbmc) 或者用默认的
## 鉴定表达高变基因(2000个）,用于下游分析,如PCA；
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
## 提取表达量变化最高的10个基因；
top10 <- head(VariableFeatures(pbmc), 10)
top10
plot1 <- VariableFeaturePlot(pbmc)
plot2 <- LabelPoints(plot = plot1, points = top10)
CombinePlots(plots = list(plot1, plot2))
plot1<-VariableFeaturePlot(object=pbmc)
plot2<-LabelPoints(plot=plot1,points=top10,repel=TRUE)
CombinePlots(plots=list(plot1,plot2))
## PCA分析：
# PCA分析数据准备,使用ScaleData()进行数据归一化；默认只是标准化高变基因（2000个）,速度更快,不影响PCA和分群,但影响热图的绘制。
#pbmc <- ScaleData(pbmc,vars.to.regress ="percent.mt")
## 而对所有基因进行标准化的方法如下：
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
pbmc <- ScaleData(pbmc, vars.to.regress = "percent.mt")
## 线性降维（PCA）,默认用高变基因集,但也可通过features参数自己指定；
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
## 展示 pca 结果（最简单的方法）
DimPlot(object=pbmc,reduction="pca")
## 检查PCA分群结果, 这里只展示前5个PC,每个PC只显示5个基因；
print(pbmc[["pca"]], dims = 1:5, nfeatures = 5)
##PC_ 1
##Positive: RPS27, MALAT1, RPS6, RPS12, RPL13
##Negative: CSTA, FCN1, CST3, LYZ, LGALS2
##PC_ 2
##Positive: NKG7, GZMA, CST7, KLRD1, CCL5
##Negative: RPL34, RPL32, RPL13, RPL39, LTB
##PC_ 3
##Positive: MS4A1, CD79A, BANK1, IGHD, CD79B
##Negative: IL7R, RPL34, S100A12, VCAN, AIF1
##PC_ 4
##Positive: RPS18, RPL39, RPS27, MALAT1, RPS8
##Negative: PPBP, PF4, GNG11, SDPR, TUBB1
##PC_ 5
##Positive: PLD4, FCER1A, LILRA4, SERPINF1, LRRC26
##Negative: MS4A1, CD79A, LINC00926, IGHD, FCER2
## 展示主成分基因分值
VizDimLoadings(pbmc, dims = 1:2, reduction = "pca")
## 绘制pca散点图
DimPlot(pbmc, reduction = "pca")
## 画第1个或15个主成分的热图；
DimHeatmap(pbmc, dims = 1, cells = 500, balanced = TRUE)
DimHeatmap(pbmc, dims = 1:15, cells = 500, balanced = TRUE)
## 确定数据集的分群个数
# 鉴定数据集的可用维度，方法1：Jackstraw置换检验算法；重复取样（原数据的1%）,重跑PCA,鉴定p-value较小的PC；计算‘null distribution’(即零假设成立时)时的基因scores。虚线以上的为可用维度，也可以调整 dims 参数，画出所有 pca 查看。
#pbmc <- JackStraw(pbmc, num.replicate = 100)
#pbmc <- ScoreJackStraw(pbmc, dims = 1:20)
#JackStrawPlot(pbmc, dims = 1:15)
# 方法2：肘部图（碎石图）,基于每个主成分对方差解释率的排名。
ElbowPlot(pbmc)
## 细胞聚类：分群个数这里选择10,建议尝试选择多个主成分个数做下游分析,对整体影响不大；在选择此参数时,建议选择偏高的数字（为了获取更多的稀有分群,“宁滥勿缺”）；有些亚群很罕见,如果没有先验知识,很难将这种大小的数据集与背景噪声区分开来。
## 非线性降维（UMAP/tSNE)基于PCA空间中的欧氏距离计算nearest neighbor graph,优化任意两个细胞间的距离权重（输入上一步得到的PC维数） 。
pbmc <- FindNeighbors(pbmc, dims = 1:10)
## 接着优化模型,resolution参数决定下游聚类分析得到的分群数,对于3K左右的细胞,设为0.4-1.2 能得到较好的结果(官方说明)；如果数据量增大,该参数也应该适当增大。
pbmc <- FindClusters(pbmc, resolution = 0.5)
## 使用Idents（）函数可查看不同细胞的分群；
head(Idents(pbmc), 5)
## 结果：AAACCTGAGGTGCTAG AAACCTGCAGGTCCAC AAACCTGCATGGAATA AAACCTGCATGGTAGG AAACCTGCATTGGCGC
1 3 0 10 2
Levels: 0 1 2 3 4 5 6 7 8 9 10 11
## Seurat提供了几种非线性降维的方法进行数据可视化（在低维空间把相似的细胞聚在一起）,比如UMAP和t-SNE,运行UMAP需要先安装'umap-learn'包，这里不做介绍，两种方法都可以使用，但不要混用，如果混用，后面的结算结果会将先前的聚类覆盖掉，只能保留一个。
## 这里采用基于TSNE的聚类方法。
pbmc <- RunTSNE(pbmc, dims = 1:10)
## 用DimPlot()函数绘制散点图,reduction = "tsne",指定绘制类型；如果不指定,默认先从搜索 umap,然后 tsne, 再然后 pca；也可以直接使用这3个函数PCAPlot()、TSNEPlot()、UMAPPlot()； cols,pt.size分别调整分组颜色和点的大小；
DimPlot(pbmc,reduction = "tsne",label = TRUE,pt.size = 1.5)
## 这里采用基于图论的聚类方法
pbmc<-RunUMAP(object=pbmc,dims=1:10)
DimPlot(object=pbmc,reduction="umap")
## 细胞周期归类
pbmc<- CellCycleScoring(object = pbmc, g2m.features = cc.genes$g2m.genes, s.features = cc.genes$s.genes)
head(x = [email protected])
DimPlot(pbmc,reduction = "tsne",label = TRUE,group.by="Phase",pt.size = 1.5)
## 存储结果
saveRDS(pbmc, file = "D:/pbmc_tutorial.rds")
save(pbmc,file="D:/res0.5.Robj")
## 寻找cluster 1的marker
cluster1.markers <- FindMarkers(pbmc, ident.1 = 1, min.pct = 0.25)
head(cluster1.markers, n = 5)
## 结果: p_val avg_logFC pct.1 pct.2 p_val_adj
MT-CO1 0.000000e+00 -0.6977083 0.985 0.996 0.000000e+00
RPS27 2.182766e-282 0.3076454 1.000 0.999 3.480202e-278
MT-CO3 2.146399e-274 -0.4866429 0.995 0.997 3.422218e-270
DUSP1 2.080878e-247 -1.7621662 0.376 0.745 3.317752e-243
RPL34 8.647733e-244 0.3367755 1.000 0.997 1.378795e-239
##寻找每一cluster的marker
pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
pbmc.markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_logFC)
# A tibble: 24 x 7
# Groups: cluster [12]
p_val avg_logFC pct.1 pct.2 p_val_adj cluster gene
<dbl> <dbl><dbl><dbl> <dbl><fct> <chr>
1 2.29e-123 0.636 0.344 0.097 3.65e-119 0 CD8B
2 7.62e-113 0.487 0.632 0.305 1.22e-108 0 LEF1
3 2.04e- 74 0.483 0.562 0.328 3.25e- 70 1 LEF1
4 1.39e- 61 0.462 0.598 0.39 2.22e- 57 1 ITM2A
5 0. 2.69 0.972 0.483 0. 2 GNLY
6 0. 2.40 0.964 0.164 0. 2 GZMB
7 1.31e-121 0.768 0.913 0.671 2.09e-117 3 JUNB
8 2.06e- 94 0.946 0.426 0.155 3.28e- 90 3 RGS1
9 2.05e-255 1.57 0.586 0.09 3.27e-251 4 GZMK
10 2.94e-140 1.57 0.69 0.253 4.68e-136 4 KLRB1
# ... with 14 more rows
## 存储marker
write.table(pbmc.markers,file="D:/allmarker.txt")
## 各种绘图
## 绘制Marker 基因的tsne图
FeaturePlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP", "CD8A"),cols = c("gray", "red"))
## 绘制Marker 基因的小提琴图
VlnPlot(pbmc, features = c("MS4A1", "CD79A"))
VlnPlot(pbmc, features = c("NKG7", "PF4"), slot = "counts", log = TRUE)
## 绘制分cluster的热图
top10 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)
DoHeatmap(pbmc, features = top10$gene) + NoLegend()
剩下的便是寻找基因 marker 并对细胞类型进行注释(见下回分解)
数据科学入门丨选Python还是R对于想入门数据科学的新手来说，选择学Python还是R语言是一个难题，本文对两种语言进行了比较，希望能帮助你做出选择。
我是德勤的数据科学家主管，多年来我一直在使用Python和R语言，并且与Python社区密切合作了15年。本文是我对这两种语言的一些个人看法。
第三种选择
针对这个问题，Studio的首席数据科学家Htley Wickham认为，比起在二者中选其一，更好的选择是让两种语言合作。因此，这也是我提到的第三种选择，我在文本最后部分会探讨。
如何比较R和Python
对于这两种语言，有以下几点值得进行比较：
· 历史：
R和Python的发展历史明显不同，同时有交错的部分。
· 用户群体：
包含许多复杂的社会学人类学因素。
· 性能：
详细比较以及为何难以比较。
· 第三方支持：
模块、代码库、可视化、存储库、组织和开发环境。
· 用例：
根据具体任务和工作类型有不同的选择。
· 是否能同时使用：
在Python中使用R，在R中使用Python。
· 预测：
内部测试。
· 企业和个人偏好：
揭晓最终答案。
历史
简史：
ABC语言 - >Python 问世(1989年由Guido van Rossum创立) - >Python 2(2000年) - >Python 3(2008年)
Fortan语言 - >S语言(贝尔实验室) - >R语言问世(1991年由Ross Ihaka和Robert Gentleman创立) - >R 1.0.0(2000年) - >R 3.0.2(2013年)
用户群体
在比较Python与R的使用群体时，要注意：
只有50％的Python用户在同时使用R。
假设使用R语言的程序员都用R进行相关“科学和数字”研究。可以确定无论程序员的水平如何，这种统计分布都是真实。
这里回到第二个问题，有哪些用户群体。整个科学和数字社区包含几个子群体，当中存在一些重叠。
使用Python或R语言的子群体：
· 深度学习
· 机器学习
· 高级分析
· 预测分析
· 统计
· 探索和数据分析
· 学术科研
· 大量计算研究领域
虽然每个领域几乎都服务于特定群体，但在统计和探索等方面，使用R语言更为普遍。在不久之前进行数据探索时，比起Python，R语言花的时间更少，而且使用Python还需要花时间进行安装。
这一切都被称为Jupyter Notebooks和Anaconda的颠覆性技术所改变。
Jupyter Notebook：增加了在浏览器中编写Python和R代码的能力
Anaconda：能够轻松安装和管理Python和R。
现在，你可以在友好的环境中启动和运行Python或R，提供开箱即用的报告和分析，这两项技术消除了完成任务和选择喜欢语言间的障碍。Python现在能以独立于平台的方式打包，并且更快地提供快速简单的分析。
社区中影响语言选择的另一个因素是“开源”。不仅仅是开源的库，还有协作社区对开源的影响。讽刺的是，Tensorflow和GNU Scientific Library等开源软件(分别是Apache和GPL)都与Python和R绑定。虽然使用R语言的用户很多，但使用Python的用户中有很多纯粹的Python支持者。另一方面，更多的企业使用R语言，特别是那些有统计学背景的。
最后，关于社区和协作，Github对Python的支持更多。如果看到最近热门的Python包，会发现Tensorflow等项目有超过3.5万的用户收藏。但看到R的热门软件包，Shiny、Stan等的收藏量则低于2千。
性能
这方面不容易进行比较。
原因是需要测试的指标和情况太多。很难在任何一个特定硬件上测试。有些操作通过其中一种语言优化，而不是另一种。
循环
在此之前让我们想想，如何比较Python与R。你真的想在R语言写很多循环吗？毕竟这两种语言的设计意图不太相同。
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as npn",
"%load_ext rpy2.ipython"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"def do_loop(u1):n",
"n",
"# Initialize `usq`n",
"usq = {}n",
"n",
"for i in range(100):n",
" # i-th element of `u1` squared into `i`-th position of `usq`n",
" usq[i] = u1[i] * u1[i]n"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"%%Rn",
"do_loop <- function(u1) {n",
"n",
"# Initialize `usq`n",
"usq <- 0n",
"n",
"for(i in 1:100) {n",
" # i-th element of `u1` squared into `i`-th position of `usq`n",
" usq[i] <- u1[i]*u1[i]n",
"}n",
"n",
"}"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"1.58 ms ± 42.8 ?s per loop (mean ± std. dev. of 7 runs, 1000 loops each)n"
]
}
],
"source": [
"%%timeit -n 1000n",
"%%Rn",
"u1 <- rnorm(100)n",
"do_loop(u1)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"36.9 ?s ± 5.99 ?s per loop (mean ± std. dev. of 7 runs, 1000 loops each)n"
]
}
],
"source": [
"%%timeit -n 1000n",
"u1 = np.random.randn(100)n",
"do_loop(u1)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.3"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Python为0.000037秒，R为0.00158秒
包括加载时间和在命令行上运行：R需要0.238秒，Python需要0.147秒。强调，这并不是科学严谨的测试。
测试证明，Python的运行速度明显加快。通常这并没有太大影响。
除了运行速度外，对于数据科学家而言哪种性能更重要？两种语言之所以受欢迎是因为它们能被用作命令语言。例如，在使用Python时大多时候我们都很依赖Pandas。这涉及到每种语言中模块和库，以及其执行方式。
第三方支持
Python有PyPI，R语言有CRAN，两者都有Anaconda。
CRAN使用内置的install.packages命令。目前，CRAN上有大约1.2万个包。其中超过1/2的包都能用于数据科学。
PyPi中包的数量超过前者的10倍，约有14.1万个包。专门用于科学工程的有3700个。其中有些也可以用于科学，但没有被标记。
在两者中都有重复的情况。当搜索“随机森林”时，PyPi中可以得到170个项目，但这些包并不相同。
尽管Python包的数量是R的10倍，但数据科学相关的包的数量大致相同。
运行速度
比较DataFrames和Pandas更有意义。
我们进行了一项实验：比较针对复杂探索任务的执行时间，结果如下：
在大多数任务中Python运行速度更快。
http://nbviewer.jupyter.org/gist/brianray/4ce15234e6ac2975b335c8d90a4b6882
可以看到，Python + Pandas比原生的R语言DataFrames更快。注意，这并不意味着Python运行更快，Pandas 是基于Numpy用C语言编写的。
可视化
这里将ggplot2与matplotlib进行比较。
matplotlib是由John D. Hunter编写的，他是我在Python社区中最敬重的人之一，他也是教会我使用Python的人。
Matplotlib虽然不易学习但能进行定制和扩展。ggplot难以进行定制，有些人认为它更难学。
如果你喜欢漂亮的图表，而且无需自定义，那么R是不错的选择。如果你要做更多的事情，那么Matplotlib甚至交互式散景都不错。同样，R的ShinnyR能够增加交互性。
是否能同时使用
可能你会问，为什么不能同时使用Python和R语言?
以下情况你可以同时使用这两种语言：
· 公司或组织允许；
· 两种都能在你的编程环境中轻松设置和维护；
· 你的代码不需要进入另一个系统；
· 不会给合作的人带来麻烦和困扰。
一起使用两种语言的方法是：
· Python提供给R的包：如rpy2、pyRserve、Rpython等；
· R也有相对的包：rPython、PythonInR、reticulate、rJython，SnakeCharmR、XRPython
· 使用Jupyter，同时使用两者，例子如下：
之后可以传递pandas的数据框，接着通过rpy2自动转换为R的数据框，并用“-i df”转换：
http://nbviewer.jupyter.org/gist/brianray/734bd54f468d9a6db9171b2cfc98405a
预测
Kaggle上有人对开发者使用R还是Python写了一个Kernel。他根据数据发现以下有趣的结果：
· 如果你打算明年转向Linux，则更可能是Python用户；
· 如果你研究统计数据，则更可能使用R；如果研究计算机科学，则更可能使用Python；
· 如果你还年轻(18-24岁)，则更可能是Python用户；
· 如果你参加编程比赛，则更可能是Python用户；
· 如果你明年想使用Android，则更可能是Python用户；
· 如果你想在明年学习SQL，则更可能是R用户；
· 如果你使用MS office，则更可能是R用户；
· 如果你想在明年使用Rasperry Pi，则更可能是Python用户；
· 如果你是全日制学生，则更可能是Python用户；
· 如果你使用的敏捷方法(Agile methodology)，则更可能是Python用户；
· 如果对待人工智能，比起兴奋你更持担心态度，则更可能是R用户。
企业和个人偏好
当我与Googler和Stack Overflow的大神级人物Alex Martelli交流时，他向我解释了为什么Google最开始只官方支持少数几种语言。即使是在Google相对开发的环境中，也存在一些限制和偏好，其他企业也是如此。
除了企业偏好，企业中第一个使用某种语言的人也会起到决定性作用。第一个在德勤使用R的人他目前仍在公司工作，目前担任首席数据科学家。我的建议是，选择你喜欢的语言，热爱你选择的语言，起到领导作用，并热爱你的事业。
当你在研究某些重要的内容时，犯错是难以避免的。然而，每个精心设计的数据科学项目都为数据科学家留有一些空间，让他们进行实验和学习。重要的是保持开放的心态，拥抱多样性。
最后就我个人而言，我主要使用Python，之后我期待学习更多R的内容。
![r语言acf[-1]是什么意思](/aiimages/r%E8%AF%AD%E8%A8%80acf%5B-1%5D%E6%98%AF%E4%BB%80%E4%B9%88%E6%84%8F%E6%80%9D.png)
