质控过滤以及标准化

生信喵

背景
       标准处理流程：读取数据后对矩阵进行标准的处理流程，包括 QC 过滤，数据标准化以及检测差异表达的基因组。
一、数据质控

       数据质控就是分别对行（基因）以及列（细胞）进行过滤，目的是过滤掉一些不符合要求的数据，例如油滴中的空细胞，多细胞，以及死细胞。
质控的参数主要有两个：
       1、每个细胞测到的 UMI 总数，过低为空细胞，过高为多细胞；同时检测表达的基因数目，过高或过低都有问题；
       2、检测每个细胞的线粒体基因的比例，理论上线粒体基因组与核基因组相比，只占很小一部分，当细胞裂解后线粒体比率会增高，cellranger 过滤中不包含该步骤，seurat 可以根据进行线粒体基因比率进行过滤。

1. #计算线粒体基因比例
   
2. pbmc[['percent.mt']] <- PercentageFeatureSet(pbmc,pattern = "^MT-")
   
3. pbmc[['percent.mt']]
   
4. #绘制小提琴图
   
5. VlnPlot(pbmc,features = c("nFeature_RNA","nCount_RNA","percent.mt"),ncol=3)

复制代码

      
       小提琴图
       nFeature_RNA 代表每个细胞测到的基因数目，nCount 代表每个细胞测到所有基因的表达量之和，percent.mt 代表测到的线粒体基因的比例。

1. #绘制散点图
   
2. plot1 <- FeatureScatter(pbmc,feature1 = "nCount_RNA",feature2 = "percent.mt")
   
3. plot2 <- FeatureScatter(pbmc,feature1 = "nCount_RNA",feature2 = "nFeature_RNA")
   
4. plot1
   
5. plot2
   
6. plot1+plot2

复制代码

      
       散点图

       去除线粒体基因表达比例过高的细胞，和一些极值细胞。
       本案例中选择表达基因在 200-3500 之间，线粒体基因比例小于 5%。其他数据集请使用不同的标准。例如当测序量增加，细胞数增多时，选择相应的参数。

1. #过滤条件根据不同测序项目进行调整
   
2. dim(pbmc) # 13714  2700
   
3. pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 3500 & percent.mt < 5)
   
4. dim(pbmc) # 13714  2643 基因数没变，细胞减少了

复制代码

      对过滤后的数据集进行可视化。

1. plot1 <- FeatureScatter(pbmc, feature1 ="nCount_RNA", feature2 = "percent.mt")+ NoLegend()
   
2. plot2 <- FeatureScatter(pbmc, feature1 ="nCount_RNA", feature2 = "nFeature_RNA")+ NoLegend()
   
3. # plot1 + plot2
   
4. CombinePlots(plots = list(plot1, plot2))

复制代码

      
       QC 之后散点图


二、标准化
       由于不同细胞中测序到不同的数据，需要对表达量进行数据标准化，标准化有多种方式，默认使用 LogNormalize 的标准化算法，还有 CLR，RC 等。
       A = log( 1 + ( UMIA ÷ UMITotal ) × 10000 )

1. pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
   
2. pbmc@commands #已经存储了标准化的数据
   
3. pbmc@commands$NormalizeData.RNA

复制代码

三、识别高变（表达差异）基因
       这一步的目的是鉴定出细胞与细胞之间表达量相差很大的基因，用于后续鉴定细胞类型，一些基因在部分细胞中高表达，部分细胞中低表达。研究表明，利用这些基因在下游分析中更容易找到关联。
       单细胞测序中的表达差异基因与 bulk RNAseq 处理组对照组比较不同，而是一个 cluster 与其他 cluster 之间的差别。

1. #默认2000个，通过nfeatures调整
   
2. pbmc1 <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
   
3. 
   
4. #前十高变基因
   
5. top10 <- head(VariableFeatures(pbmc1), 10)
   
6. top10
   
7. #可视化
   
8. plot1 <- VariableFeaturePlot(pbmc1)+guides(color="none")
   
9. plot1
   
10. plot2 <- LabelPoints(plot = plot1,points = top10,repel = TRUE,xnudge = 0,ynudge = 0)
    
11. plot2
    
12. plot1+plot2
    
13. 
    
14. #默认2000个，通过nfeatures调整成22000超过所有基因13000多。
    
15. pbmc1 <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 22000)
    
16. 
    
17. #前十高变基因
    
18. top10 <- head(VariableFeatures(pbmc1), 10)
    
19. top10
    
20. #可视化
    
21. plot1 <- VariableFeaturePlot(pbmc1)+guides(color="none")
    
22. plot1
    
23. plot2 <- LabelPoints(plot = plot1,points = top10,repel = TRUE,xnudge = 0,ynudge = 0)
    
24. plot2
    
25. plot1+plot2

复制代码