单细胞day3

原创

用户10300557

发布于 2024-06-19 15:28:15

2270

文章被收录于专栏：生信学习111生信学习111

1 加载数据

还是放在原来的文件夹下没有动

load("obj.Rdata") #载入数据
> p1 = DimPlot(seu.obj,reduction = "umap",label = T)
> p1
> rm(list = ls()) #清空列表区域省的乱
> library(Seurat)#加载这个Seurat这个包
> load("obj.Rdata") #载入数据
> p1 = DimPlot(seu.obj,reduction = "umap",label = T) #昨天画图的口令，有一个label = true
> p1

这张图的每一个点都是一个细胞，同一个颜色的点被认为时一类细胞，那末到底是什么细胞呢，可以通过marker基因进行分析。

2marker 基因

可以把marker基因理解为特征基因，特征基因在某一类细胞中表达量比较高，在其它类细胞中表达量低。

理论上，上调和下调都算marker基因，但我们会设置 only.pos = TRUE，只关注表达量变高的。

FindAllMarkers，又是一个限速步骤，随着细胞数量的增大会更慢，它会计算出所有细胞簇的marker基因，比如计算0这一簇时，就是0 vs 剩下的1-10，1-10全部算一个组。我理解的就是找某一类和其他剩余类别的不同，工作量不小。

还有一个参数min.pct，min.pct = 0.25的意思是，以0簇为例，就是说要么在0这一簇里25%以上的细胞中表达，要么在剩下的全部簇里25%以上的细胞中表达。这样可以去掉一些表达比例太少（就是表达占比不高）的基因，比例太少的话不能代表这个簇的整体情况，其实很好理解就是你这个在这个群体里占比比较高才可以说你是marker基因。

Seurat开发了一个新的方法加快FindAllMarkers的计算速度，在初次运行时它会给出提示让你安装这个包，我的给了提示证明我并没有安装过，可以安装一下来提提速

(没有提示就说明已经安装过啦)。

2.1 R语言知识补充github的R包安装

有一些包你可能特别想用，但是不在各种管理网站上，而是放在开发者的github仓库里。这样的R包也是可以安装和使用的。

两种安装方法

1 在线安装,注意安装的时presto这个包

if(!require(devtools))install.packages("devtools")
if(!require(presto))devtools::install_github("immunogenomics/presto",upgrade = F,dependencies = T)

2 本地下载一下

在github网站下载压缩包放在工作目录下,然后运行这个口令，这样安装可能会有一些依赖包安不上，那还得在进行安装，注意看报错

if(!require(presto))devtools::install_github("presto-master.zip",upgrade = F,dependencies = T)

2.2 取每个细胞簇logFC排名前n的marker基因。logFC用于衡量两个样本或条件下基因表达水平的相对变化

> mks = markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_log2FC)
> g = unique(mks$gene)

n=2就是每个簇取前2,这个可以自己设置，看情况吧，然后为什么要g = unique(mks$gene)这个呢，因为可能有些基因在不同簇中都被记为marker基因，有重复会报错，因此需要去除一下。

3可视化

3.1 热图

> library(ggplot2)
> DoHeatmap(seu.obj, features = g) + NoLegend()+
+   scale_fill_gradientn(colors = c("#2fa1dd", "green","#f87669"))
Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.

3.2 气泡图

 DotPlot(seu.obj, features = g,cols = "RdYlBu") +
  RotatedAxis()

3.3 小提琴图，每一个基因都会有一个小提琴，所以少画点

VlnPlot(seu.obj, features = g[1:6]) #提取了前六个

3.4 Featureplot，本质上还是umap图，只是换了一种着色方式，一张图对应一个基因，展示了该基因在所有细胞里的表达情况。

3.5 峰峦图,峰峦图，或者是山脊图，其实就是密度图

RidgePlot(seu.obj, features = g[1:2])

4 注释亚群

手动注释很复杂，手动注释需要自行查找每种可能的细胞类型已知的marker基因，可以从文献里查，可以从数据库里查，真的很复杂得自己看，还是自动注释把，自动注释也有问题：没有把组织、疾病类型等背景知识考虑在内 > 现成的参考数据只有人类和小鼠的 > 切换不同的参考数据就会给出不同的注释结果

4.2自动注释

常用的自动注释R包是singleR，使用的参考数据是celldex里面的

花花老师的有7个，但是我有15个很奇怪。解密了是更新啦

library(celldex)
library(SingleR)
ls("package:celldex")
z = "ref_Human_all.RData" #花花老师用的是f = "ref_BlueprintEncode.RData"
if(!file.exists(z)){
  ref <- celldex::ref_Human_allData() #上面用啥这块改啥保留RData前面的内容，放在Data前面
  save(ref,file = z)
}
ref <- get(load(z))
library(BiocParallel)
library(SingleR) #报错显示没找到，得加载一下
scRNA = seu.obj
test = scRNA@assays$RNA$data
pred.scRNA <- SingleR(test = test, 
                      ref = ref,
                      labels = ref$label.main, 
                      clusters = scRNA@active.ident)
pred.scRNA$pruned.labels #看一下分类
new.cluster.ids <- pred.scRNA$pruned.labels
names(new.cluster.ids) <- levels(scRNA)
library(Seurat) #"RenameIdents"是Seurat里面的
scRNA <- RenameIdents(scRNA,new.cluster.ids)
p3 <- DimPlot(scRNA, reduction = "umap",label = T,pt.size = 0.5) + NoLegend()
p3 #没操作手动注释的所以没有p2