10a. Identify marker peaks in the integrated manually annotated clusters

Using the manually called clusters with pseudo-bulk replicates, identify marker peaks based on these cluster identifiers.

# Distribution of cells across manual cell annotation clusters
fwrite(as.data.table(table(proj5$ManualClusters)),
       "scATAC-seq_ManualCluster_CellFrequency.csv", sep=",")

markersPeaks <- getMarkerFeatures(
    ArchRProj = proj5, 
    useMatrix = "PeakMatrix", 
    groupBy = "ManualClusters",
    bias = c("TSSEnrichment", "log10(nFrags)"),
    testMethod = "wilcoxon"    
)

markersPeaks
markerList <- getMarkers(markersPeaks, cutOff = "FDR <= 0.01 & Log2FC >= 1")
markerList

# Save cell type marker peak lists
for (celltype in names(markerList)) {
    name <- paste0(celltype,
                   "_scATAC-seq_Integrated_ManualClusters_marker_peaks_",
                   Sys.Date(), ".csv")
    fwrite(as.data.table(markerList[[celltype]]), name, sep=",")
}

Plot expression UMAPs for genes of interest to the Gomez group

genes_of_interest <- c("Hopx", "Tulp1", "Bhlhe41", "Zscan26", "Gata3", "Tbx2",
                       "Litaf", "Camta1", "Mgp", "Pbx1", "Atf3", "Foxd1",
                       "Mef2c", "Creb1", "Ren1", "Akr1b7", "Myh11", "Cnn1",
                       "Acta2", "Itgb1", "Foxd1", "Mef2a", "Mef2b", "Mef2c",
                       "Mef2d", "Smarcc1", "Nfix", "Rbpj", "Tal1", "Tagln",
                       "Piezo1", "Piezo2", "Hoxb7", "Six2", "Cited1", "Smtn",
                       "Bach1", "Bach2", "Jund", "Junb", "Jun", "Fos", "Nfic")
for (gene in genes_of_interest) {
    if (!file.exists(paste0(
        sample_name, "_scATAC-seq_GeneIntegrationMatrix_",
        gene, "_", Sys.Date(), ".svg"))) {
        p1 <- plotEmbedding(
            ArchRProj = proj5, 
            colorBy = "GeneIntegrationMatrix_LSI_Harmony", 
            name = gene, 
            continuousSet = "horizonExtra",
            embedding = "UMAP_LSI_Harmony",
            imputeWeights = getImputeWeights(proj5)
        )

    p2 <- plotEmbedding(
        ArchRProj = proj5, 
        colorBy = "GeneScoreMatrix_no_ribosomal", 
        continuousSet = "horizonExtra",
        name = gene, 
        embedding = "UMAP_LSI_Harmony",
        imputeWeights = getImputeWeights(proj5)
    )

    svglite(paste0(sample_name,
                   "_scATAC-seq_GeneIntegrationMatrix_", gene,
                   "_", Sys.Date(), ".svg"))
    print(p1)
    dev.off()

    svglite(paste0(sample_name,
                   "_scATAC-seq_GeneScoreMatrix_", gene,
                   "_", Sys.Date(), ".svg"))
    print(p2)
    dev.off()
    }
}