背景：單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)的質(zhì)量控制萧锉，質(zhì)量控制去除（1）低質(zhì)量的細(xì)胞（2）雙胞（3）RNA污染的細(xì)胞

• 1.去除低質(zhì)量的細(xì)胞
  什么樣的細(xì)胞是低質(zhì)量的細(xì)胞呢？通常是檢測到的基因數(shù)目少膀篮、計數(shù)深度低或者線粒體等基因的比例比較高的細(xì)胞。
  去除低質(zhì)量細(xì)胞的閾值對象：
  （1）每個barcode（細(xì)胞）的計數(shù)數(shù)量（計數(shù)深度、測序深度）
  （2）每個barcode的基因數(shù)量
  （3）每個barcode的線粒體基因（核糖體有勾、血紅蛋白基因）比例

library(Seurat)
library(ggplot2)
library(ggpubr)
sc<-readRDS("./sc.merge.rds")
sc<-CreateSeuratObject(sc$RNA@counts,meta.data=sc@meta.data,min.cells=3) #'篩選基因，保留在大于3個細(xì)胞中表達(dá)的基因
sc[["percent.mt"]] <- PercentageFeatureSet(sc, pattern="^MT-") #'人類中線粒體基因以MT開頭古程，因此計算MT開頭的基因比例
sc[["percent.Ribo"]] <- PercentageFeatureSet(sc, pattern="^RP[SL]") #'計算以RPS或者RPL開頭的基因比例（核糖體基因）
HB.genes <- c("HBA1","HBA2","HBB","HBD","HBE1","HBG1","HBG2","HBM","HBQ1","HBZ")
HB.genes <- CaseMatch(HB.genes, rownames(sc))
sc[["percent.HB"]] <- PercentageFeatureSet(sc, features=HB.genes) #'血紅蛋白基因比例

qc_before<- VlnPlot(sc, group.by="orig.ident",
                              features=c("nFeature_RNA","nCount_RNA","percent.mt","percent.Ribo","percent.HB"),
                              pt.size=0,
                              ncol=2)
ggsave(qc_before,file="qc_before.png", width=12, height=6)
#'質(zhì)控前的相關(guān)閾值情況
Feature_ber1 <- FeatureScatter(sc, feature1="nCount_RNA",
                                   feature2="nFeature_RNA",
                                   group.by="orig.ident")
Feature_ber2 <- FeatureScatter(sc, feature1 = 'nCount_RNA',
                                   feature2="percent.mt",
                                   group.by="orig.ident")
Feature_ber3 <- FeatureScatter(sc, feature1 = 'nCount_RNA',
                                   feature2="percent.Ribo",
                                   group.by="orig.ident")
Feature_ber4 <- FeatureScatter(sc, feature1 = 'nCount_RNA',
                                   feature2="percent.HB",
                                   group.by="orig.ident")
Feature_ber1 <- Feature_ber1 +theme(legend.position="none")
Feature_ber2 <- Feature_ber2 +theme(legend.position="none")
Feature_ber3 <- Feature_ber3 +theme(legend.position="none")
Feature_ber4 <- Feature_ber4 +theme(legend.position="none")
Feature_ber <- ggarrange(Feature_ber1, Feature_ber2, Feature_ber3,Feature_ber4, ncol=4, nrow=1, widths=c(1,1,1,1))
ggsave(Feature_ber, file="Feature_ber_before.png", width=12, height=6)

qc_before.png

Feature_ber_before.png

開始質(zhì)控：

sc.sub=subset(sc,nFeature_RNA > 200 & nFeature_RNA < 2500 & nCount_RNA > 500 & nCount_RNA < 10000 & percent.mt < 5 &  percent.HB<5 & percent.Ribo<10) #'篩選細(xì)胞蔼卡，規(guī)則是細(xì)胞的基因數(shù)目在200-2500之間，測序深度（基因表達(dá)總量）在500-10000挣磨，線粒體比例低于5%雇逞，血紅蛋白基因低>于5%荤懂，核糖體基因低于10%
qc_after<- VlnPlot(sc.sub, group.by="orig.ident",
                              features=c("nFeature_RNA","nCount_RNA","percent.mt","percent.Ribo","percent.HB"),
                              pt.size=0,
                              ncol=2)
ggsave(qc_after,file="qc_after.png", width=12, height=6)
Feature_ber1 <- FeatureScatter(sc.sub, feature1="nCount_RNA",
                                   feature2="nFeature_RNA",
                                   group.by="orig.ident")
Feature_ber2 <- FeatureScatter(sc.sub, feature1 = 'nCount_RNA',
                                   feature2="percent.mt",
                                   group.by="orig.ident")
Feature_ber3 <- FeatureScatter(sc.sub, feature1 = 'nCount_RNA',
                                   feature2="percent.Ribo",
                                   group.by="orig.ident")
Feature_ber4 <- FeatureScatter(sc.sub, feature1 = 'nCount_RNA',
                                   feature2="percent.HB",
                                   group.by="orig.ident")
Feature_ber1 <- Feature_ber1 +theme(legend.position="none")
Feature_ber2 <- Feature_ber2 +theme(legend.position="none")
Feature_ber3 <- Feature_ber3 +theme(legend.position="none")
Feature_ber4 <- Feature_ber4 +theme(legend.position="none")
Feature_ber <- ggarrange(Feature_ber1, Feature_ber2, Feature_ber3,Feature_ber4, ncol=4, nrow=1, widths=c(1,1,1,1))
ggsave(Feature_ber, file="Feature_ber_after.png", width=12, height=6)

qc_after.png

可以看到留下來的細(xì)胞的基因數(shù)目都在200-2500之間，基因表達(dá)總量在300-10000之間喝峦，線粒體比例低于5%势誊，核糖體基因比例低于10%。

Feature_ber_after.png

• 2.去除雙胞
  “雙細(xì)胞被定義為在相同的細(xì)胞條形碼（barcodes）下進(jìn)行測序的兩個細(xì)胞谣蠢，例如粟耻，兩個細(xì)胞被捕獲在同一個液滴（droplet）中。這也是為什么我們一直使用barcodes而不是cells的原因眉踱〖访Γ”除了使用算法進(jìn)行雙細(xì)胞的識別外，一些別的指標(biāo)也能進(jìn)行提示：例如谈喳，同一個聚類類群中的marker基因來自不同的細(xì)胞類型册烈、某個類群中的基因數(shù)目或者基因表達(dá)量異常高......
  這里測試的是DoubletFinder：

#'obj the seurat obj
#'本過程不會對原始o(jì)bj作任何數(shù)據(jù)預(yù)處理
myDoublet<-function(obj,assay="RNA",outpath="./",group="seurat_clusters",dims=30,resolution=0.5,replace=F){
library(Seurat)
library(DoubletFinder)
library(dplyr)
ls("package:DoubletFinder")
sc<-obj
#sc$group=sc@meta.data[,group]
if(assay=="RNA"){
DefaultAssay(sc)<-"RNA"
data <- sc %>%NormalizeData() %>% FindVariableFeatures(nfeatures=2000) %>% ScaleData() %>% RunPCA(assay="RNA",verbose=FALSE)%>%FindNeighbors(reduction = "pca", dims = 1:dims)%>%FindClusters(verbose = FALSE,resolution=resolution)%>%RunTSNE(dim.use=1:dims)
}else{
DefaultAssay(sc)<-"Spatial"
data <- sc %>%SCTransform(sc, assay = "Spatial",do.scale=TRUE, do.center=TRUE,verbose = FALSE) %>% RunPCA()%>%RunTSNE(dim.use=1:dims)
}
sweep.res.list <- paramSweep(data, PCs = 1:30, sct = F,num.cores=6)
sweep.stats <- summarizeSweep(sweep.res.list, GT = FALSE)
bcmvn <- find.pK(sweep.stats) #'可以看到最佳參數(shù)的點
pK_bcmvn <- bcmvn$pK[which.max(bcmvn$BCmetric)] %>% as.character() %>% as.numeric() #'提取最佳pk值

DoubletRate = ncol(data)*8*1e-6 #'每增加1000個細(xì)胞，雙胞概率增加0.008
homotypic.prop <- modelHomotypic(data@meta.data[,group]) #'默認(rèn)seurat_clusters
nExp_poi <- round(DoubletRate*ncol(data))
nExp_poi.adj <- round(nExp_poi*(1-homotypic.prop))
data <- doubletFinder(data, PCs = 1:30, pN = 0.25, pK = pK_bcmvn,
                         nExp = nExp_poi, reuse.pANN = F, sct = F)
data <- doubletFinder(data, PCs = 1:30, pN = 0.25, pK = pK_bcmvn,
                         nExp = nExp_poi.adj, reuse.pANN = F, sct = F)

data@meta.data[,"DF_hi.lo"] <- data@meta.data[,length(colnames(data@meta.data))]
data@meta.data$DF_hi.lo[which(data@meta.data$DF_hi.lo == "Doublet" & data@meta.data[,length(colnames(data@meta.data))-3] == "Singlet")] <- "Doublet-Low Confidience"
data@meta.data$DF_hi.lo[which(data@meta.data$DF_hi.lo == "Doublet")] <- "Doublet-High Confidience"
table(data@meta.data$DF_hi.lo) #'查看細(xì)胞類型（雙胞或者單胞）分布
#'結(jié)果可視化
png(paste0(outpath,"01.",assay,"_",group,"_tsne.png"),2500,1800,res=300)
pp<-DimPlot(data, reduction = "tsne", group.by =group)
print(pp)
dev.off()
png(paste0(outpath,"01.",assay,"_",group,"_doubletFinder.png"),2500,1800,res=300)
pp<-DimPlot(data, reduction = "tsne", group.by ="DF_hi.lo",cols =c("black","red","gold"))
print(pp)
dev.off()
saveRDS(data,paste0(outpath,'01.sc_rmdoublet.rds'))#'保存數(shù)據(jù)的rds文件
if(replace==T){
sc@meta.data$DF_hi.lo<-data@meta.data[match(rownames(sc@meta.data),rownames(data@meta.data)),"DF_hi.lo"]
return(sc)
}

01.doubletFinder.png

01.RNA_seurat_clusters_tsne.png

注意：理論上去除雙胞應(yīng)該在上一步細(xì)胞的基因數(shù)目等閾值質(zhì)控之前進(jìn)行婿禽。

• 3.去除周期相關(guān)的基因

#'去除周期基因的影響
#'對于assay="RNA"
library(Seurat)
library(dplyr)
sc<-readRDS("./sc.merge.rds")
sc<-sc%>%NormalizeData()%>%ScaleData(features=rownames(sc))%>%FindVariableFeatures() %>% RunPCA() #'數(shù)據(jù)預(yù)處理
#'首先不管是sct還是log標(biāo)準(zhǔn)化赏僧，隨便選一個方法跑到RunPCA
#'然后使用seurat中內(nèi)置的資源 周期基因cc.genes
cc.genes<-cc.genes.updated.2019
g2m_genes <- cc.genes$g2m.genes
g2m_genes <- CaseMatch(search=g2m_genes, match=rownames(sc))
s_genes <- cc.genes$s.genes
s_genes <- CaseMatch(search=s_genes, match=rownames(sc))
#'這里cc.genes默認(rèn)是人類的大寫基因
#'但是小鼠的即使沒有轉(zhuǎn)換大小寫也能自動匹配
sc <- CellCycleScoring(object=sc, g2m.features=g2m_genes, s.features=s_genes,set.ident=TRUE) #'計算細(xì)胞周期
#'meta.data中多了幾列S.Score G2M.Score Phase old.ident
#'可視化一下周期基因pca
cycleplot <- DimPlot(sc, reduction="pca", group.by="Phase")
ggsave(cycleplot, file="cycleplot_before_1.png", width=12, height=6)
sc <- RunPCA(sc, features=c(s_genes,g2m_genes))
cycleplot <- DimPlot(sc, reduction="pca", group.by="Phase")
ggsave(cycleplot, file="cycleplot_before_2.png", width=12, height=6)
#'對周期基因進(jìn)行回歸處理 關(guān)鍵函數(shù)是scale
sc<-sc%>%NormalizeData()%>%ScaleData(features=rownames(sc),vars.to.regress = c("S.Score", "G2M.Score"))%>%FindVariableFeatures()
sc<-RunPCA(sc,assay="RNA",verbose=FALSE)
#'這里需要注意一下，如果是空轉(zhuǎn)的數(shù)據(jù)扭倾，應(yīng)該走sct
#'sc<-SCTransform(sc, assay = "Spatial", vars.to.regress = c("S.Score", "G2M.Score"),do.scale=TRUE, do.center=TRUE,verbose = FALSE)
if(T){ #'周期基因可視化
cyclesce<-sc
cyclesce<- RunPCA(cyclesce, assay="RNA",features = VariableFeatures(cyclesce), nfeatures.print = 10)
cycleplot<-DimPlot(cyclesce,cols = pal_npg("nrc", alpha = 0.7)(3))
ggsave(cycleplot, file="cycleplot_after_1.png", width=130, height=100,units="mm",limitsize=FALSE)
cyclesce<- RunPCA(cyclesce, assay=assay,features = c(s.genes, g2m.genes))
cycleplot<-DimPlot(cyclesce,cols = pal_npg("nrc", alpha = 0.7)(3))
ggsave(cycleplot, file="cycleplot_after_2.png", width=130, height=100,units="mm",limitsize=FALSE)
}

01.RNA_cycleplot_before_1.png

01.RNA_cycleplot_after_1.png

注意：謹(jǐn)慎使用淀零。

• 去除RNA污染
  方法一：DecontX

myDecontX<-function(obj,assay="RNA",outpath="./",my.off=0.2,dims=30,resolution=0.5,group="seurat_clusters"){
library(Seurat)
library(decontX)
library(ggplot2)
library(dplyr)
sc<-obj
if(assay=="RNA"){
my.counts<-sc$RNA@counts
}else{
my.counts<-sc$Spatial@counts
}
sc<-sc%>%SCTransform(verbose=F,do.scale=T)%>% RunPCA(assay="SCT",verbose=FALSE)%>%FindNeighbors(reduction = "pca", dims = 1:dims)%>%FindClusters(verbose = FALSE,resolution=resolution)%>%RunUMAP(dims=1:dims)
decontX_results <- decontX(my.counts)
str(decontX_results)
sc$contamination=decontX_results$contamination
head(sc@meta.data)
png(paste0(outpath,"01.",assay,"_decontx_contamination_before.png"),2500,1800,res=300)
pp<-FeaturePlot(sc,
            features = 'contamination',
            raster=FALSE       # 細(xì)胞過多時候需要加這個參數(shù)
           ) +
  scale_color_viridis_c()+
  theme_bw()+
  theme(panel.grid = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank())
print(pp)
dev.off()
png(paste0(outpath,"01.",assay,"_decontx_umap_before.png"),2500,1800,res=300)
pp<-DimPlot(sc, reduction = "umap", group.by =group)
print(pp)
dev.off()
sc = sc[,sc$contamination < my.off]
png(paste0(outpath,"01.",assay,"_decontx_contamination_afetr.png"),2500,1800,res=300)
pp<-FeaturePlot(sc,
            features = 'contamination',
            raster=FALSE       # 細(xì)胞過多時候需要加這個參數(shù)
           ) +
  scale_color_viridis_c()+
  theme_bw()+
  theme(panel.grid = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank())
print(pp)
dev.off()
png(paste0(outpath,"01.",assay,"_decontx_umap_afetr.png"),2500,1800,res=300)
pp<-DimPlot(sc, reduction = "umap", group.by =group)
print(pp)
dev.off()
if(assay=="Spatial"){
print(paste0("the filtered-spot-gene counts's dim is gene:",nrow(sc$Spatial$counts),";spot:",ncol(sc$Spatial$counts)))
nf.before=nrow(obj$Spatial$counts)
nf.after=nrow(sc$Spatial$counts)
ns.before=ncol(obj$Spatial$counts)
ns.after=ncol(sc$Spatial$counts)
spot.s=paste0(as.character(round(ns.after/ns.before, 3)*100),"%")
feature.s=paste0(as.character(round(nf.after/nf.before, 3)*100),"%")
qc.table=data.frame(type_number=c("nSpot","nFeature"),filter_before=c(as.character(ns.before),as.character(nf.before)),filter_after=c(as.character(ns.after),as.character(nf.after)),percent=c(spot.s,feature.s))
}else{
print(paste0("the filtered-cell-gene counts's dim is gene:",nrow(sc$RNA$counts),";cell:",ncol(sc$RNA$counts)))
nf.before=nrow(obj$RNA$counts)
nf.after=nrow(sc$RNA$counts)
ns.before=ncol(obj$RNA$counts)
ns.after=ncol(sc$RNA$counts)
spot.s=paste0(as.character(round(ns.after/ns.before, 3)*100),"%")
feature.s=paste0(as.character(round(nf.after/nf.before, 3)*100),"%")
qc.table=data.frame(type_number=c("nCell","nFeature"),filter_before=c(as.character(ns.before),as.character(nf.before)),filter_after=c(as.character(ns.after),as.character(nf.after)),percent=c(spot.s,feature.s))
}
write.csv(qc.table,paste0(outpath,"01.",assay,"_decontx_qc.csv"),row.names=F)
print(head(qc.table))
return(sc)
}

和soupx不同之處在于，不會修改原始矩陣表達(dá)膛壹，只是根據(jù)閾值對細(xì)胞進(jìn)行過濾驾中，因此會減少細(xì)胞的數(shù)目。
方法二：soupx

#'raw.path the 10x raw matrix
#'fit.path the 10x filtered matrix
mySoupX<-function(raw.path,fit.path,outpath="./",my.off=0.2,assay="RNA",group="seurat_clusters",dims=30,resolution=0.5){
library(SoupX)
library(Seurat)
library(DropletUtils)
library(ggplot2)
library(DoubletFinder)
library(knitr)
library(dplyr)
filt.matrix <- Read10X(fit.path)
raw.matrix  <- Read10X(raw.path)
raw.matrix<-raw.matrix[rownames(filt.matrix),]
srat  <- CreateSeuratObject(counts = filt.matrix)
soup.channel  <- SoupChannel(raw.matrix, filt.matrix)
srat <- SCTransform(srat, verbose = F,do.scale=T)
srat <- RunPCA(srat,assay="SCT", verbose = F)
srat <- RunUMAP(srat, dims = 1:dims, verbose = F)
srat <- FindNeighbors(srat,reduction="pca", dims = 1:dims, verbose = F)
srat <- FindClusters(srat, verbose = F,resolution=resolution)
meta <- srat@meta.data
umap <- srat@reductions$umap@cell.embeddings
soup.channel <- setClusters(soup.channel, setNames(meta$seurat_clusters, rownames(meta)))
soup.channel <- setDR(soup.channel, umap)
head(meta)
soup.channel <- autoEstCont(soup.channel)
soup.channel<-setContaminationFraction(soup.channel, my.off) #'閾值my.off
adj.matrix  <- adjustCounts(soup.channel, roundToInt = T) #'校正矩陣
png(paste0(outpath,"01.",assay,"_soupx_before.png"),2500,1800,res=300)
pp<-DimPlot(srat, reduction = "umap", group.by =group)
print(pp)
dev.off()
srat_1<-CreateSeuratObject(counts = adj.matrix)
srat_1<-srat_1%>%SCTransform(verbose=F,do.scale=T)%>% RunPCA(assay="SCT",verbose=FALSE)%>%FindNeighbors(reduction = "pca", dims = 1:dims)%>%FindClusters(verbose = FALSE,resolution=resolution)%>%RunUMAP(dims=1:dims)
png(paste0(outpath,"01.",assay,"_soupx_after.png"),2500,1800,res=300)
pp<-DimPlot(srat_1, reduction = "umap", group.by =group)
print(pp)
dev.off()
DropletUtils:::write10xCounts(paste0(outpath,"01.soupX_pbmc10k_filt"), adj.matrix) #'保存矩陣
}

相同的參數(shù)下模聋，可視化矩陣校正前后：

01.RNA_soupx_before.png

01.RNA_soupx_after.png

因為group.by用的是seurat_clusters,且沒有相關(guān)的背景知識肩民，所以不太能看出來優(yōu)化之處。
注意：修改了矩陣链方，所以應(yīng)該謹(jǐn)慎使用持痰。