参考：生信技能树公众号、单细胞天地公众号、周运来就是我（简书）。 备注：仅作为笔记学习，无任何商业用途。

过滤不合格的细胞和基因

现在一起看看拿到表达矩阵后有一个很重要的质控过程，就是根据一些阈值来过滤不合格细胞和基因。细胞的不合格很容易理解， 可能是那个细胞检测到的基因数量太少，或者太多，也许是那个细胞的文库大小异常，都是需要谨慎考虑是否管理它。
基因的不合格，就需要大家对人类参考基因组的注释信息有所背景了。走RNA-seq定量流程，拿到的表达矩阵通常是取决于gtf文件的注释程度，人类的gtf里面五万多个基因，不可能都在你的单细胞转录组项目数据里面出现。

细胞的取舍

主要看文库大小和检测基因数量。
单细胞转录组通常使用10x数据，所以细胞数量惊人，每个样本可以是3000到10000的细胞数量都可以，取决于实验设计。每个细胞平均可以是1到10万的reads文库，都没有问题。那么，每个细胞可以检测到多少基因数量呢，就取决于每个细胞分到的reads总数，1到10万的reads文库对应着200到1000的基因数量。
但是Smart-seq2技术的单细胞数据就不一样了，每个样本细胞数量通常是96的倍数，500个左右就很厉害了，然后每个细胞的reads就可以很多，百万级别都没有问题。所以检测到基因数量也很多，如下：
image.png
所以在设置阈值来过滤不合格细胞的时候，一定要先理解好单细胞转录组数据，给全部的细胞检测到的基因数量绘制一个boxplot，看看哪些细胞所检测到的基因数量偏多或者偏少，一般来说，这样的的离群细胞就是需要去除的！比如jimmy在单细胞转录组教学就介绍过一个代码：

  1. box <- lapply(colnames(sample_ann[,1:19]),function(i) {
  2.     dat <-  sample_ann[, i, drop = F] 
  3.     dat$sample = rownames(dat)
  4.     ## 画boxplot 
  5.    ggplot(dat, aes('all cells', get(i))) +
  6.           geom_boxplot() +
  7.           xlab(NULL) + ylab(i)
  8. })
  9. plot_grid(plotlist = box, ncol = 5 )

删除一些基因

前面提到过，gtf文件的注释程度不一样，拿到的表达矩阵通常是全部的基因，比如人类的gtf里面五万多个基因，表达矩阵就是5万多行。
实际上，大量的基因在所有的细胞里面都是不表达的，这就是表达量全部为0的，肯定是需要去除啦。
还有一些基因仅仅是在5%左右的细胞表达，这个时候就很主观了，可能确实这个基因是那5%的稀有细胞的marker基因，也有可能就是单细胞转录组技术导致这个基因大量的drop-out了。
这个时候，就需要卡一个阈值，到底一个基因表达多少算是表达呢，到底一个基因在多少个细胞里面有表达，你才保留下去，做下游分析呢？

过滤线粒体和核糖体基因

我们埋下了一个伏笔，就是拿到了差不多干净的表达矩阵后还有一个预处理步骤，但是这个步骤是选修，就是说你不会过滤线粒体核糖体基因也可以不做，没有人会责怪你，因为单细胞转录组数据分析本来就有难度。分析要点无穷无尽，即使你真的准备过滤线粒体核糖体基因你也会发现标准不好把握，需要看很多文献，还得结合你自己项目数据的实际情况啦。而且过滤线粒体核糖体基因并不是质控的终点，你还有细胞周期检查，单细胞活性检查，是否有两个细胞混在一起也是需要检查，非常的累。不仅仅是标准不好把握，而且因项目而异，具体情况具体判断。

直接删除

直接在表达矩阵里面，去除掉属于的那一行表达量即可，有点类似于甲基化芯片数据分析，直接去除性染色体上面的全部探针，如下所示：
image.png

选择阈值

这个阈值很大程度上取决于你对自己项目的了解程度，不同器官组织提取的单细胞，本来就线粒体基因平均水平不一样， 不能一刀切！这个只能是靠你看文献来获取认知了，多看多学习！我的经验是，多次反复查看线粒体核糖体基因的影响，分群前后看，不同batch看，多次质控图表里面显示它，判断它是否是一个主要因素。
image.png

线粒体和核糖体基因

我这里以gencode数据库的gtf文件为标准，在人类的gencode.v32.annotation.gtf 文件里面，可以查找到37个:

  1. gene_id "ENSG00000210049.1"; gene_name "MT-TF"; hgnc_id "HGNC:7481"
  2. gene_id "ENSG00000211459.2"; gene_name "MT-RNR1"; hgnc_id "HGNC:7470"
  3. gene_id "ENSG00000210077.1"; gene_name "MT-TV"; hgnc_id "HGNC:7500"
  4. gene_id "ENSG00000210082.2"; gene_name "MT-RNR2"; hgnc_id "HGNC:7471"
  5. gene_id "ENSG00000209082.1"; gene_name "MT-TL1"; hgnc_id "HGNC:7490"
  6. gene_id "ENSG00000198888.2"; gene_name "MT-ND1"; hgnc_id "HGNC:7455"
  7. gene_id "ENSG00000210100.1"; gene_name "MT-TI"; hgnc_id "HGNC:7488"
  8. gene_id "ENSG00000210107.1"; gene_name "MT-TQ"; hgnc_id "HGNC:7495"
  9. gene_id "ENSG00000210112.1"; gene_name "MT-TM"; hgnc_id "HGNC:7492"
  10. gene_id "ENSG00000198763.3"; gene_name "MT-ND2"; hgnc_id "HGNC:7456"
  11. gene_id "ENSG00000210117.1"; gene_name "MT-TW"; hgnc_id "HGNC:7501"
  12. gene_id "ENSG00000210127.1"; gene_name "MT-TA"; hgnc_id "HGNC:7475"
  13. gene_id "ENSG00000210135.1"; gene_name "MT-TN"; hgnc_id "HGNC:7493"
  14. gene_id "ENSG00000210140.1"; gene_name "MT-TC"; hgnc_id "HGNC:7477"
  15. gene_id "ENSG00000210144.1"; gene_name "MT-TY"; hgnc_id "HGNC:7502"
  16. gene_id "ENSG00000198804.2"; gene_name "MT-CO1"; hgnc_id "HGNC:7419"
  17. gene_id "ENSG00000210151.2"; gene_name "MT-TS1"; hgnc_id "HGNC:7497"
  18. gene_id "ENSG00000210154.1"; gene_name "MT-TD"; hgnc_id "HGNC:7478"
  19. gene_id "ENSG00000198712.1"; gene_name "MT-CO2"; hgnc_id "HGNC:7421"
  20. gene_id "ENSG00000210156.1"; gene_name "MT-TK"; hgnc_id "HGNC:7489"
  21. gene_id "ENSG00000228253.1"; gene_name "MT-ATP8"; hgnc_id "HGNC:7415"
  22. gene_id "ENSG00000198899.2"; gene_name "MT-ATP6"; hgnc_id "HGNC:7414"
  23. gene_id "ENSG00000198938.2"; gene_name "MT-CO3"; hgnc_id "HGNC:7422"
  24. gene_id "ENSG00000210164.1"; gene_name "MT-TG"; hgnc_id "HGNC:7486"
  25. gene_id "ENSG00000198840.2"; gene_name "MT-ND3"; hgnc_id "HGNC:7458"
  26. gene_id "ENSG00000210174.1"; gene_name "MT-TR"; hgnc_id "HGNC:7496"
  27. gene_id "ENSG00000212907.2"; gene_name "MT-ND4L"; hgnc_id "HGNC:7460"
  28. gene_id "ENSG00000198886.2"; gene_name "MT-ND4"; hgnc_id "HGNC:7459"
  29. gene_id "ENSG00000210176.1"; gene_name "MT-TH"; hgnc_id "HGNC:7487"
  30. gene_id "ENSG00000210184.1"; gene_name "MT-TS2"; hgnc_id "HGNC:7498"
  31. gene_id "ENSG00000210191.1"; gene_name "MT-TL2"; hgnc_id "HGNC:7491"
  32. gene_id "ENSG00000198786.2"; gene_name "MT-ND5"; hgnc_id "HGNC:7461"
  33. gene_id "ENSG00000198695.2"; gene_name "MT-ND6"; hgnc_id "HGNC:7462"
  34. gene_id "ENSG00000210194.1"; gene_name "MT-TE"; hgnc_id "HGNC:7479"
  35. gene_id "ENSG00000198727.2"; gene_name "MT-CYB"; hgnc_id "HGNC:7427"
  36. gene_id "ENSG00000210195.2"; gene_name "MT-TT"; hgnc_id "HGNC:7499"
  37. gene_id "ENSG00000210196.2"; gene_name "MT-TP"; hgnc_id "HGNC:7494"

所以你拿到自己的表达矩阵后，其实简单的看看基因名字是否以 MT- 开头即可哦。
所以你会看到seurat包的函数：PercentageFeatureSet 可以用来计算线粒体基因含量。

  1. sce[["percent.mt"]] <- PercentageFeatureSet(sce, pattern = "^MT-")
  2. VlnPlot(sce, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

上面的方法是修改 sce[[“percent.mt”]] ，下面我们演示 AddMetaData 函数，同样是可以增加线粒体基因含量信息到我们的seurat对象。

  1. mt.genes <- rownames(sce)[grep("^MT-", rownames(sce))]
  2. C <- GetAssayData(object = sce, slot = "counts")
  3. percent.mito <- Matrix::colSums(C[mt.genes,])/Matrix::colSums(C)*100
  4. sce <- AddMetaData(sce, percent.mito, col.name = "percent.mito")
  5. sce[["percent.mito"]]

也可以是添加核糖体基因含量，同样的你需要知道核糖体基因的名字规则：

  1. rb.genes <- rownames(sce)[grep("^RP[SL]", rownames(sce))]
  2. percent.ribo <- Matrix::colSums(C[rb.genes,])/Matrix::colSums(C)*100
  3. sce <- AddMetaData(sce, percent.ribo, col.name = "percent.ribo")

注：小鼠和人的基因名称不一样，小鼠是小写。