source <- c(10930,10318,10595,10972,7706,6756,9092,10551,9722,10913,11151,8186,6422,
6337,11649,11652,10310,12043,7937,6476,9662,9570,9981,9331,9449,6773,6304,9355,10477,
10148,10395,11261,8713,7299,10424,10795,11069,11602,11427,9095,7707,10767,12136,12812,
12006,12528,10329,7818,11719,11683,12603,11495,13670,11337,10232,13261,13230,15535,
16837,19598,14823,11622,19391,18177,19994,14723,15694,13248,9543,12872,13101,15053,
12619,13749,10228,9725,14729,12518,14564,15085,14722,11999,9390,13481,14795,15845,
15271,14686,11054,10395,14775,14618,16029,15231,14246,12095,10473,15323,15381,14947)
srcLen<-length(source)
for(i in 1:10){ #预测最后十个数；
real <- source[srcLen-i+1] #实际值
xNum=(srcLen-i+1)%/%7 #组数
yNum=7 #每组7个数
data<-array(1:(xNum*yNum),c(xNum,yNum))
pre=srcLen-i+1
for(x in 1:xNum){ #数组赋值
for(y in 1:yNum){
data[x,y]=source[pre]
pre=pre-1
}
if(pre<7){
break
}
}
ascData<-array(1:(xNum*yNum),c(xNum,yNum)) #数组逆序
for(x in 1:xNum){
for(y in 1:yNum){
ascData[x,y]=data[xNum-x+1,y]
}
}
colnames(ascData) <- c("a","b","c","d","e","f","g") #每列列名
trainData<-data.frame(scale(ascData[,c(1:7)]))
nn<-nnet(a~b+c+d+e+f+g,trainData[1:(xNum-1),],size=10,decay=0.01,maxit=1000,linout=F,trace=F)
predict<-predict(nn,trainData[xNum,])
predict=predict*sd(ascData[,1])+mean(ascData[,1])
percent <- (predict-real)*100/real
res <- paste("预测值：",predict,"实际值：",real,"误差：",percent)
print(res)
}
第一步：获取要绘图的整洁数据（涉及到数据整洁和操作的知识）第二步：整洁数据做映射操作，确定x,y,color,size,shape,alpha等
第三步：选择合适的几何对象（根据画图的目的、变量的类型和个数）
第四步：坐标系和刻度配置
第五步：标签信息和图例信息
第六步：选择合适的主题
ggplot2的语法包括10个部件。
数据（data）
映射（mapping）
几何对象（geom）
标度（scale）
统计变换（stats）
坐标系（coord）
位置调整（Position adjustments）
分面（facet）
主题（theme）
输出（output）
前3个是必须的，其它部件ggplot2会自动配置，也可以手动配置
ggplot2基本绘图模板：
注意：
1）添加图层的加号（+）只能放在行末尾
2）红色方框里面mapping是全局域，绿色方框里面mapping是局部域，执行先后顺序，先局部域，后全局域
ggplot2画图必要部件-数据，映射和几何对象
2.1 数据
数据（Data）用于画图的整洁数据
library(tidyverse
ggplot()先只提供数据，创建一个空图形。
# ggplot()先提供整洁数据，生成一个空图形
2映射
映射，把数据变量集与图形属性库建立关联。
最常用的映射有：
x：x轴
y：y轴
color：颜色
size：大小
shape：形状
fill：填充
alpha：透明度
以mpg数据集为例，把变量displ和hwy分别映射到x和y，变量drv映射到color，此时图形就有了坐标轴和网格线，color需要在有了几何对象后才能体现出来。
# 映射操作
ggplot(data = mpg, mapping = aes(x = displ,
y = hwy, color = drv))
2.3 几何对象
几何对象是表达数据的视觉对象
不同类型的几何对象是从不同的角度表达数据。
pgglot2提供了50多种“几何对象”，均以geom_xxxx()的方式命名，常用的有：
几何对象很简单，只需要添加图层即可。
例如，以mpg数据集为例，画散点图。
ggplot(data = mpg, mapping = aes(x = displ,
y = hwy,
color = drv)) +
geom_point()层依次叠加，在上图的基础上，再添加一个几何对象：光滑曲线。
#继续增加一个几何对象：光滑曲线
# 写法1
ggplot(data = mpg, mapping = aes(x = displ,
y = hwy,
color = drv)) +
geom_point() +
geom_smooth(se=FALSE)
# 写法2
ggplot(data = mpg, mapping = aes(x = displ, y = hwy)) +
geom_point(aes(color = drv)) +
geom_smooth(se=FALSE)
思考题：
1）写法1和写法2的差异？（全局域和局部域的使用差异）
2）写法2若是要实现写法1的功能，怎么编写代码？
03
标度
ggplot2会自动根据输入变量选择最优的坐标刻度方法，若要手动设置或调整，就需要使用标度函数。
标度函数用来控制几何对象中的标度映射（x轴，y轴或者由color，fill，shape，size产生的图例）。
ggplot2提供丰富的标度函数，常用的有：
拓展功能：scales包提供很多设置刻度标签风格的函数，比如百分数、科学计数法法、美元格式等。
3.1 修改坐标轴刻度及标签
连续变量使用scale_*_continuous()函数，参数breaks设置各个刻度的位置，参数labels设置各个刻度对应的标签。
离散变量使用scale_*_discrete()函数，修改离散变量坐标轴的标签。
时间变量使用scale_x_date()函数设置日期刻度，参数date_breaks设置刻度间隔，date_labels设置标签的日期格式
以mpg数据集为例，修改连续变量坐标轴刻度及标签。
# scale_y_continuous函数
# 对比分析和观察
# 图1
ggplot(mpg, aes(displ, hwy)) +
geom_point()
# 图2
ggplot(mpg, aes(displ, hwy)) +
geom_point() +
scale_y_continuous(breaks = seq(15, 40, by = 10))
# 图3
ggplot(mpg, aes(displ, hwy)) +
geom_point() +
scale_y_continuous(breaks = seq(15, 40, by = 10),
labels = c(" 一五 "," 二五 "," 三五 "))
以mpg数据集为例，修改离散变量的标签
# scale_x_discrete函数
# 对比分析和观察
# 图1
ggplot(mpg, aes(x = drv)) +
geom_bar()
# 图2
ggplot(mpg, aes(x = drv)) +
geom_bar() +
scale_x_discrete(labels = c("4" = " 四驱 ", "f" = " 前驱 ",
"r" = " 后驱 "))
以ggplot2自带的economics数据集为例，修改日期变量。
# scale_x_date函数
# 以ggplot2自带的economics为例
economics %>% glimpse()
# 图1
ggplot(tail(economics, 45), aes(date, uempmed / 100)) +
geom_line()
# 图2
ggplot(tail(economics, 45), aes(date, uempmed / 100)) +
geom_line() +
scale_x_date(date_breaks = "6 months", date_labels = "%Y-%b") +
scale_y_continuous(labels = scales::percent)
3.2 修改坐标轴标签、图例名及图例位置
用labs()函数参数x，y或者xlab()，ylab()，设置x轴，y轴标签。
若用参数color生成了图例，可以在labs()函数用参数color修改图例名。
用theme图层的参数legend.position设置图例的位置。
以mpg数据为例。
# 修改坐标轴标签，图例名和图例位置
mpg
# 图1
ggplot(mpg, aes(displ, hwy)) +
geom_point(aes(color = drv)) +
labs(x = " 引擎大小 (L)", y = " 高速燃油率 (mpg)",
color = " 驱动类型 ") +
theme(legend.position = "top")
# 图2
ggplot(mpg, aes(displ, hwy)) +
geom_point(aes(color = drv)) +
xlab(" 引擎大小 (L)") +
ylab(" 高速燃油率 (mpg)") +
labs(color = " 驱动类型 ") +
theme(legend.position = "top")
# 图3 不需要图例
ggplot(mpg, aes(displ, hwy)) +
geom_point(aes(color = drv)) +
xlab(" 引擎大小 (L)") +
ylab(" 高速燃油率 (mpg)") +
theme(legend.position = "none")
3.3 设置坐标轴的范围
用coord_cartesian()函数参数xlim和ylim，或者用xlim()，ylim()设置x轴和y轴的范围。
以mpg数据集为例。
# 修改坐标轴的范围
# 图1 coord_cartesian()的参数xlim和ylim
ggplot(mpg, aes(displ, hwy)) +
geom_point(aes(color = drv)) +
coord_cartesian(xlim = c(5, 7), ylim = c(10, 30))
# 图2 xlim()和ylim()函数
ggplot(mpg, aes(displ, hwy)) +
geom_point(aes(color = drv)) +
xlim(5, 7) +
ylim(10, 30)
3.4 变换坐标轴
用scale_x_log10()函数变换坐标系，可以保持原始数据的坐标刻度。
# 修改坐标轴的范围
# 图1 coord_cartesian()的参数xlim和ylim
ggplot(mpg, aes(displ, hwy)) +
geom_point(aes(color = drv)) +
coord_cartesian(xlim = c(5, 7), ylim = c(10, 30))
# 图2 xlim()和ylim()函数
ggplot(mpg, aes(displ, hwy)) +
geom_point(aes(color = drv)) +
xlim(5, 7) +
ylim(10, 30)
3.5 设置图形标题
用labs()函数设置图形标题。
参数title 设置正标题
参数subtitle 设置副标题
参数caption 设置脚注标题（默认右下角）
# 设置标题
# mpg数据集为例
p <- ggplot(mpg, aes(displ, hwy)) +
geom_point(aes(color = drv)) +
geom_smooth(se = FALSE) +
labs(title = " 燃油效率与引擎大小的关系图 ",
subtitle = " 两座车 ( 跑车 ) 因重量小而符合预期 ",
caption = " 数据来自 fueleconomy.gov")
p
标题若要居中，采用theme图层设置。
p + theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5))
3.6 设置color、fill颜色
数据的某个维度信息可以通过颜色来表示。
可以直接使用颜色值，建议使用RColorBrewer（调色板）或者colorspace包。
1）连续变量
- 用scale_color_gradient()设置二色渐变色。
# 连续变量
# 图1 scale_color_gradient()函数
ggplot(mpg, aes(displ, hwy, color = hwy)) +
geom_point() +
scale_color_gradient(low = "green", high = "red")
- 用scale_color_distiller()设置调色板中的颜色
# 图2 scale_color_distiller()函数
ggplot(mpg, aes(displ, hwy, color = hwy)) +
geom_point() +
scale_color_distiller(palette = "Set1")
2）离散变量
- 用scale_color_manual()手动设置颜色，还可以修改图例及其标签信息
# 离散变量
# 图1 scale_color_manual()函数
ggplot(mpg, aes(displ, hwy, color = drv)) +
geom_point() +
scale_color_manual(" 驱动方式 ",
values = c("red", "blue", "green"),
breaks = c("4", "f", "r"))
ggplot(mpg, aes(displ, hwy, color = drv)) +
geom_point() +
scale_color_manual(" 驱动方式 ",
values = c("red", "blue", "green"),
labels = c(" 四驱 ", " 前驱 ", " 后驱 "))
-用scale_fill_brewer()调用调色板中的颜色
# 图2 scale_fill_brewer()函数
ggplot(mpg, aes(x = class, fill = class)) +
geom_bar() +
scale_fill_brewer(palette = "Dark2")
.7 添加文字标注
ggrepel包提供了geom_label_repel()函数或者geom_text_repel()函数，为图形添加文字标注。
操作步骤：
第一步：先准备好标记点的数据
第二步：增加文字标注图层，包括标记点的数据和标注的文字给label参数
# 设置文字标注信息
library(ggrepel)
# 选取每种车型 hwy 值最大的样本
best_in_class <- mpg %>%
group_by(class) %>%
slice_max(hwy, n = 1)
best_in_class %>% select(class, model, hwy)
ggplot(mpg, aes(displ, hwy)) +
geom_point(aes(color = class)) +
geom_label_repel(data = best_in_class,
aes(label = model))
04
计变换、坐标系和位置调整
.1 统计变换
统计变换是构建新的统计量而画图。
例如，条形图或直方图，是对数据分组的频数做画图；平滑曲线是对数据拟合模型的预测值画图。
gplot2可以把统计变换直接融入画图中，不必先在对数据做统计变换后再画图。
gplot2提供30多种统计，均以stats_xxx()的方式命名。
1）可在几何对象中直接使用的统计变换，直接使用几何对象就可以了。
能在几何对象创建的，而需要单独使用。
mpg数据集为例。
stat_summary()做统计绘图并汇总。
# 图1 stat_summary()做统计绘图并汇总
p <- ggplot(mpg, aes(x = class, y = hwy)) +
geom_violin(trim = FALSE, alpha = 0.5, color = "green")
p
p + stat_summary(fun = mean,
fun.min = function (x) {mean(x) - sd(x)},
fun.max = function (x) {mean(x) + sd(x)},
geom = "pointrange",
color = "red")
tat_smooth()添加光滑曲线，与geom_smooth()相同。
参数method设置平滑曲线的拟合方法，如lm线性回归、glm广义线性回归、loess多项式回归、gam广义加法模型(mgcv包)、rlm稳健回归(MASS包)等。
参数formula指定平滑曲线方程，如y ~ x, y ~ poly(x, 2)， y ~ log(x)等。
参数se设置是否绘制置信区间。
# 图2 stat_smooth()添加平滑曲线
ggplot(mpg, aes(displ, hwy)) +
geom_point() +
stat_smooth(method = "lm",
formula = y ~ splines::bs(x, 3),
se = FALSE)
ggplot(mpg, aes(displ, hwy)) +
geom_point() +
geom_smooth(method = "lm",
formula = y ~ splines::bs(x, 3),
se = FALSE)
4.2 坐标系
ggplot2默认是直角坐标系。
- coord_cartesian()
常用的其它坐标系：
以mpg数据集为例，坐标轴翻转。
# 图1 坐标轴翻转coord_flip()
p <- ggplot(mpg, aes(class, hwy)) +
geom_boxplot()
p
p + coord_flip()
直角坐标下条形图转换为极坐标下玫瑰图。
# 图2 直角坐标条形图-->极坐标玫瑰图
p <- ggplot(mpg, aes(class, fill = drv)) +
geom_bar()
p
p + coord_polar()
4.3 位置调整
条形图的位置调整
# 图1：条形图条形位置调整
ggplot(mpg, aes(class, fill = drv)) +
geom_bar()
ggplot(mpg, aes(class, fill = drv)) +
geom_bar(position = "dodge")
ggplot(mpg, aes(class, fill = drv)) +
geom_bar(position = position_dodge(preserve = "single"))
散点图的散点位置调整
# 图1：散点图的散点位置调整
ggplot(mpg, aes(displ, hwy)) +
geom_point()
ggplot(mpg, aes(displ, hwy)) +
geom_point(position = "jitter")
用patchwork包排布多个图形
library(patchwork)
p1 <- ggplot(mpg, aes(displ, hwy)) +
geom_point()
p2 <- ggplot(mpg, aes(drv, displ)) +
geom_boxplot()
p3 <- ggplot(mpg, aes(drv)) +
geom_bar()
p1 | (p2 / p3)
p1 | p2 | p3
p1 / p2 / p3
p1 / (p2 | p3)
05
分面
利用分类变量把图形分成若干“子图”（面），实际上就是对数据分组后再画图，属于数据分析里面细分和下钻的思想。
5.1 用facet_wrap()函数
封装分面，先生成一维的面板系列，再封装到二维中。
语法形式：~ 分类变量 或者 ~ 分类变量1 + 分类变量2
参数scales设置是否共用坐标刻度，fixed 默认 共用， free 不共用，还可以额通过free_x，free_y单独设置。
下游分析
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 并对细胞类型进行注释(见下回分解)
