代码库3-Signac代码

library(Signac)
counts <- Read10X_h5("../vignette_data/atac_v1_pbmc_10k_filtered_peak_bc_matrix.h5")
fragment_file <- '../vignette_data/atac_v1_pbmc_10k_fragments.tsv.gz'
chrom_assay <- CreateChromatinAssay(counts,sep = c(":", "-"),fragments = fragment_file)
pbmc.ATAC <- CreateSeuratObject(chrom_assay,assay = "ATAC")
# pbmc.ATAC <- CreateSeuratObject(counts, assay = "ATAC")

pbmc.ATAC <- RunTFIDF(pbmc.ATAC)
pbmc.ATAC <- FindTopFeatures(pbmc.ATAC, min.cutoff = 'q0')
pbmc.ATAC <- RunSVD(pbmc.ATAC)
DepthCor(pbmc.ATAC)
pbmc.ATAC <- RunUMAP(pbmc.ATAC, dims = 2:50, reduction = 'lsi', 
                     reduction.name = 'atac_umap', reduction.key = 'atac_umap')

pbmc.ATAC <- FindNeighbors(object = pbmc.ATAC, reduction = 'lsi', dims = 2:30)
pbmc.ATAC <- FindClusters(object = pbmc.ATAC, verbose = FALSE, algorithm = 3)
DimPlot(object = pbmc.ATAC, label = TRUE) + NoLegend()

library(EnsDb.Hsapiens.v75)
annotations <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v75)、
seqlevelsStyle(annotations) <- 'UCSC'
genome(annotations) <- "hg19"
Annotation(pbmc.ATAC) <- annotations
gene.activities <- GeneActivity(pbmc.ATAC)
pbmc.ATAC[['ACTIVITY']] <- CreateAssayObject(counts = gene.activities)
pbmc.ATAC <- NormalizeData(pbmc.ATAC,assay = 'ACTIVITY')
DefaultAssay(pbmc.ATAC) <- 'ACTIVITY'
pbmc.ATAC <- FindVariableFeatures(pbmc.ATAC, selection.method = "vst", nfeatures = 2000)          
LabelPoints(plot = VariableFeaturePlot(pbmc.ATAC), points = head(VariableFeatures(pbmc.ATAC), 10), repel = TRUE)
pbmc.ATAC <- ScaleData(pbmc.ATAC)    
pbmc.ATAC <- RunPCA(pbmc.ATAC)   
ElbowPlot(pbmc.ATAC)
pbmc.ATAC <- RunUMAP(pbmc.ATAC, dims = 1:20, reduction.name = 'ga_umap', reduction.key = 'ga_umap')   
DimPlot(pbmc.ATAC, reduction = "ga_umap",group.by = 'true',label = T) + NoLegend()

Idents(pbmc.ATAC) <- 'true'
CoveragePlot(pbmc.ATAC,region = c('chr5-54070562-54071062','chr1-212806552-212807052'),
             extend.upstream = 40000,extend.downstream = 20000)