scATAC联合scRNA之signac分析(四): 单样本分析数据质控及降维聚类流程

granges(HC)# GRanges object with 170827 ranges and 0 metadata columns:#   seqnames        ranges strand# <Rle>     <IRanges>  <Rle>#   [1]       chr1 629474-630381      *#   [2]       chr1 633581-634476      *#   [3]       chr1 778268-779188      *#   [4]       chr1 816903-817798      *#   [5]       chr1 819473-820427      *#   ...        ...           ...    ...# [170823] GL000219.1   50733-51595      *#   [170824] GL000219.1  99273-100155      *#   [170825] GL000219.1 165071-165956      *#   [170826] GL000219.1 168659-169301      *#   [170827] KI270713.1   21470-22359      *#   -------#   seqinfo: 28 sequences from an unspecified genome; no seqlengths
#可以看到，features的基因组范围，前面表示染色体的部分，有些不是标准的chr，也就是EnsDb.Hsapiens.v86种没有#而是，chromosome scaffolds，e.g. (KI270713.1)，所以去除一下HC <- HC[as.vector(seqnames(granges(HC)) %in% standardChromosomes(granges(HC))), ]AA <- AA[as.vector(seqnames(granges(AA)) %in% standardChromosomes(granges(AA))), ]SD <- SD[as.vector(seqnames(granges(SD)) %in% standardChromosomes(granges(SD))), ]
#添加基因组注释# extract gene annotations from EnsDb# 基因和基因组注释信息是后续计算TSS富集得分，核小体含量和基因活跃度得分所必需的。# 同样，你得保证这里选择的基因组版本得和你上游数据处理时用到的基因组版本一致。library(AnnotationHub)ah <- AnnotationHub()query(ah, "EnsDb.Hsapiens.v98")#Search for the Ensembl 98 EnsDb for Homo sapiens on AnnotationHubensdb_v98 <- ah[["AH75011"]]#这个版本的注释文件可以保存一下，如果后续其他数据分析还用得到的话# save(ensdb_v98, file = 'ensdb_v98.RData')annotations <- GetGRangesFromEnsDb(ensdb = ensdb_v98)# annotation1 <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v86)seqlevels(annotations) <- paste0('chr', seqlevels(annotations))genome(annotations) <- "hg38"
# add the gene information to the objectAnnotation(HC) <- annotationsAnnotation(AA) <- annotationsAnnotation(SD) <- annotations

#接下来计算质控指标atac = list(HC,AA,SD)names(atac) <- c("HC","AA","SD")
#计算nucleosome signal score，TSS enrichment score#add fraction of reads in peaks,# add blacklist ratiofor(i in seq_along(atac)){
  atac[[i]] <- NucleosomeSignal(object = atac[[i]])  atac[[i]] <- TSSEnrichment(object = atac[[i]])  atac[[i]]$pct_reads_in_peaks <- atac[[i]]$peak_region_fragments / atac[[i]]$passed_filters * 100  atac[[i]]$blacklist_ratio <- FractionCountsInRegion(object = atac[[i]], assay = 'peaks',regions = blacklist_hg38_unified)
}
#plot一下nCount_peaks和TSS.enrichment密度图，并且展示5%，10%，90%，95%百分位上的数值，可以作为质控参考DensityScatter(atac[["HC"]], x = 'nCount_peaks', y = 'TSS.enrichment', log_x = TRUE, quantiles = TRUE)DensityScatter(atac[["AA"]], x = 'nCount_peaks', y = 'TSS.enrichment', log_x = TRUE, quantiles = TRUE)DensityScatter(atac[["SD"]], x = 'nCount_peaks', y = 'TSS.enrichment', log_x = TRUE, quantiles = TRUE)

#按照nucleosome_signal高低对细胞进行分组，并plot不同fragment频率分布图，也可以作为质控参考for(i in seq_along(atac)){  atac[[i]]$nucleosome_group <- ifelse(atac[[i]]$nucleosome_signal > 5, 'NS > 5', 'NS <5')}
FragmentHistogram(object = atac[["HC"]], group.by = 'nucleosome_group')FragmentHistogram(object = atac[["AA"]], group.by = 'nucleosome_group')FragmentHistogram(object = atac[["SD"]], group.by = 'nucleosome_group')

#plot VlnplotVlnPlot(object = atac[["SD"]],        features = c('peak_region_fragments', 'TSS.enrichment', 'blacklist_ratio', 'nucleosome_signal', 'pct_reads_in_peaks'),        pt.size = 0.1,        ncol = 5)


#QCfor(i in seq_along(atac)){
  atac[[i]] <- subset(x = atac[[i]],                      subset = peak_region_fragments > 1000 &                        peak_region_fragments < 20000 &                        pct_reads_in_peaks > 15 &                        blacklist_ratio < 0.005 &                        nucleosome_signal < 4 &                        TSS.enrichment > 3)
}

#展示质控后结果数据，这里我们与ArchR梦幻联动，主要是我觉得ArchR质控图可视化看着还可以# devtools::install_github("GreenleafLab/ArchR", ref="master", repos = BiocManager::repositories())
#TSS & nCount_peaks 密度图library(ArchR)library(cowplot)
densityplot_list <- list()for (i in seq_along(atac)) {
  df <- data.frame("nFrags"=log10(atac[[i]]$nCount_peaks),"TSSEnrichment" = atac[[i]]$TSS.enrichment)  p <- ggPoint(    x = df[,1],    y = df[,2],    colorDensity = TRUE,    continuousSet = "sambaNight",    xlabel = "Log10(Unique Fragments)",    ylabel = "TSS Enrichment")+     geom_hline(yintercept = 3, lty = "dashed") +     geom_vline(xintercept = log10(min(atac[[i]]$nCount_peaks)), lty = "dashed")+    ggtitle(paste0(names(atac)[i],"\n","Cells Pass Filter =",ncol(atac[[i]])))
  densityplot_list[[i]] <- p}

plot_grid(densityplot_list[[1]], densityplot_list[[2]],densityplot_list[[3]],ncol = 3)

#fragment频率分布图fix(FragmentHistogram)FragmentHistogram(object = atac[["HC"]], mycols = "#2488F0")+ggtitle("HC")+FragmentHistogram(object = atac[["AA"]], mycols = "#51127CFF")+ggtitle("AA")+FragmentHistogram(object = atac[["SD"]], mycols = "#B51F29")+ggtitle("SD")

for(i in seq_along(atac)){    atac[[i]] <- RunTFIDF(atac[[i]])#归一化处理  atac[[i]] <- FindTopFeatures(atac[[i]], min.cutoff = 'q5')#特征值选择，对应scRNA就是高变基因  atac[[i]] <- RunSVD(atac[[i]])  }for(i in seq_along(atac)){    atac[[i]] <- RunUMAP(object = atac[[i]], reduction = 'lsi', dims = 2:30)  atac[[i]] <- FindNeighbors(object = atac[[i]], reduction = 'lsi', dims = 2:30)  atac[[i]] <- FindClusters(object = atac[[i]], verbose = FALSE, algorithm = 3, resolution = 0.8)#resolution也可以多设置，自行调整  }

DimPlot(atac[['HC']], label = T)+ggtitle('HC')+DimPlot(atac[['AA']], label = T)+ggtitle('AA')+DimPlot(atac[['SD']], label = T)+ggtitle('SD')
#add gene activity matrix
gene.activities <- list()for(i in seq_along(atac)){  gene.activities[[i]] <- GeneActivity(atac[[i]])}
#将ene activity matrix添加到seurat，添加一个新的assay RNA。for(i in seq_along(atac)){
  atac[[i]][['RNA']] <- CreateAssayObject(counts = gene.activities[[i]])  atac[[i]] <- NormalizeData(object = atac[[i]],assay = 'RNA',normalization.method = 'LogNormalize',scale.factor = median(atac[[i]]$nCount_RNA))    }