单细胞数据挖掘 P2.3 挑选hvg基因，可视化

「生信技能树」单细胞数据挖掘_哔哩哔哩_bilibili生信技能树又一单细胞系列课程上线啦~小伙伴们快来学习呀https://www.bilibili.com/video/BV1pa4y1s76J?p=6

一、构建hvg并查看是否有MT/ERCC基因混杂情况

hvg基因为高变化基因，即在各个样本中，表达量差异最为明显的基因。

#highly Variable gene:简单理解sd大的
  scRNA <- FindVariableFeatures(scRNA, selection.method = "vst", nfeatures = 1500) 
  #根据文献原图，挑选变化最大的1500个hvg
  top10 <- head(VariableFeatures(scRNA), 10) 
  top10
  plot1 <- VariableFeaturePlot(scRNA) 
  #标记top10 hvg
  p1_3 <- LabelPoints(plot = plot1, points = top10, repel = TRUE, size=2.5) +
    theme(legend.position = c(0.1,0.8)) +
    labs(tag = "C")
  p1_3

findvariablefeatures函数是seurat包中的一个函数，其提取出的高变基因作为相关信息也是作为一个参数存储在scRNA矩阵中的。

nfeatures决定选出几个基因。挑选前10个看一下分布。

用tag来标记挑出的10个基因。

为了查看差异基因会不会有ERCC/MT这种线粒体基因，进行检查，并进行标注。

这里需要先加载ggplot2才能进行绘图。

 #看看ERCC
  ERCC <- rownames(scRNA)[grep("^ERCC-",rownames(scRNA))]
  LabelPoints(plot = plot1, points = ERCC, repel = TRUE, 
              size=2.5,colour = "blue") +
    theme(legend.position = c(0.1,0.8)) +
    labs(tag = "C")
  #可以直观看到ERCC均不是高变基因，而且部分的ERCC基因表达量确实很高
  p1_1 | p1_2 | p1_3 #上图

---

二、标准化

后续作者使用了scater包进行标准化了。

scater包会产生一个sce对象，这是独立于seurat产生的对象。

初学SingleCellExperiment对象 - 简书https://www.jianshu.com/p/9bba0214844b

sce构建的是一个SingleCellExperiment对象

核心分为4块：

1、assay：核心部分（盖了几层楼）：表达矩阵部分
2、colData：sample/cell annotation：meta信息样本注释部分，细胞相关信息
3、rowData：feature annotation：基因信息
4、reducedDims：PCA，tSNE，uMAP等细胞降维聚类信息，降维信息

后续是先进行log转化，标准化，之后再用plot来进行基因占不同样本的可视化。

分别使用CPM方法以及lognormalcounts来进行。主要是提取了senrat构成的表达矩阵（scRNA）中的counts以及metadata进行标准化。

这里用的标准化方法是：

CPM标准化，再进行+1的log转换。

 # 这里展开介绍一下 scater 
  # https://bioconductor.org/packages/release/bioc/html/scater.html
  # Single-cell analysis toolkit工具箱 for expression 
  #scater contains tools to help with the analysis of single-cell transcriptomic data, 
  #focusing on low-level steps such as quality control, normalization and visualization.
  #based on the SingleCellExperiment class (from the SingleCellExperiment package)
  #关于sce对象，https://www.jianshu.com/p/9bba0214844b
  library(scater)
  ct=as.data.frame(scRNA@assays$RNA@counts)
  pheno_data=scRNA@meta.data
  sce <- SingleCellExperiment(
    assays = list(counts = ct), 
    colData = pheno_data
  )
  #SingleCellExperiment是SingleCellExperiment包的函数；在加载scater包时会一起加载
  sce
  ?stand_exprs
  stand_exprs(sce) <- log2(
    calculateCPM(sce) + 1)  #只对自己的文库的标准化
  assays(sce)
  sum(counts(sce)[,1])
  head(counts(sce)[,1])
  log2(1*10^6/422507+1)
  #logcounts(sce)[1:4,1:4]
  #exprs(sce)[1:4,1:4]
  stand_exprs(sce)[1:4,1:4]

CPM标准化

每个细胞中的文库的量的标准化。

原理为：细胞中的某基因的表达量占总细胞的总基因的表达量，后再+1，防止log之后0的出现。

后面几句是对calculateCPM的分解步骤说明，直接使用calculate也可以直接产出。

  sce <- logNormCounts(sce) #可以考虑不同细胞的文库差异的标准化
  assays(sce)
  logcounts(sce)[1:4,1:4]
  
  #https://osca.bioconductor.org/normalization.html#spike-norm
  #基于ERCC的标准化方式也有许多优势（不同cell的量理论上是一样的），详见链接
  #关于一些常见的FPKM等方式在番外篇会有简单的介绍与学习

lognormcounts

lognormcounts是对不同的细胞间的文库差异进行标准化

此处使用calculateCPM以及lognormcounts是直接在sce这个对象中加入了两个两个参数，

CPM产生的是stand_exprs，log产生的是logcounts。

在运行两个标准化之后，sce会产生相应标准化之后的三个list。

使用head获取头部数据，如果表达量是1的这个基因，果然在后续标准化后的表达量中，与算出的表达量一致。

---

三、可视化

 

  #观察上面确定的top10基因在四个样本的分布比较
  # top10是一个list是用来指代sce这个文件中的需要展示的基因
  plotExpression(sce, top10 ,
                 x = "Patient_ID",  colour_by = "Patient_ID", 
                 exprs_values = "logcounts") 
  # 下面的绘图非常耗时：(保存为本地文件查看比较高效，建议为pdf文件)
  p1 <- plotHighestExprs(sce, exprs_values = "logcounts")
  ggsave("out/2.3HighestExprs.pdf", plot = p1, width = 15, height = 18) 
  #如果按照ERCC 40%的过滤标准，ERCC表达量也十分大
  ?plotHighestExprs
  # Sometimens few spike-in transcripts may also be present here, 
  # though if all of the spike-ins are in the top 50, 
  # it suggests that too much spike-in RNA was added
  
  #后续步骤暂时还按照pctERCC=40的过滤标准的结果进行分析
  
  #tips:后面主要还是基于Seurat对象，此处可以删除该变量，节约内存。
  save(scRNA, file="2.3.Rdata")
  rm(list=ls())