gravatar for i.am.filippov

2 hours ago by

I was following Seurat tutorial on clustering single-cell ATAC data. Here's the code for loading the data:

input_data <- Read10X_h5(filename = "./data/pbmc_granulocyte_sorted_3k_filtered_feature_bc_matrix.h5")
atac_counts <- input_data$Peaks

metadata <- read.csv(
  file = "./data/pbmc_granulocyte_sorted_3k_per_barcode_metrics.csv",
  header = TRUE,
  row.names = 1
)

chrom_assay <- CreateChromatinAssay(
  counts = atac_counts,
  sep = c(":", "-"),
  genome = 'hg38',
  fragments = './data/pbmc_granulocyte_sorted_3k_atac_fragments.tsv.gz',
  min.cells = 10,
  min.features = 200
)

pbmc <- CreateSeuratObject(
  counts = chrom_assay,
  assay = "peaks",
  meta.data = metadata
)

Some filtering code...and then:

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

I can see the clusters on UMAP plot, so they're there. However, I can't see the clusters in metadata:

colnames([email protected])
 [1] "orig.ident"                        "nCount_peaks"                      "nFeature_peaks"                   
 [4] "gex_barcode"                       "atac_barcode"                      "is_cell"                          
 [7] "excluded_reason"                   "gex_raw_reads"                     "gex_mapped_reads"                 
[10] "gex_conf_intergenic_reads"         "gex_conf_exonic_reads"             "gex_conf_intronic_reads"          
[13] "gex_conf_exonic_unique_reads"      "gex_conf_exonic_antisense_reads"   "gex_conf_exonic_dup_reads"        
[16] "gex_exonic_umis"                   "gex_conf_intronic_unique_reads"    "gex_conf_intronic_antisense_reads"
[19] "gex_conf_intronic_dup_reads"       "gex_intronic_umis"                 "gex_conf_txomic_unique_reads"     
[22] "gex_umis_count"                    "gex_genes_count"                   "atac_raw_reads"                   
[25] "atac_unmapped_reads"               "atac_lowmapq"                      "atac_dup_reads"                   
[28] "atac_chimeric_reads"               "atac_mitochondrial_reads"          "atac_fragments"                   
[31] "atac_TSS_fragments"                "atac_peak_region_fragments"        "atac_peak_region_cutsites"        
[34] "nucleosome_signal"                 "nucleosome_percentile"             "TSS.enrichment"                   
[37] "TSS.percentile"                    "pct_reads_in_peaks"                "high.tss"

Thanks in advance!



Source link