15a. Calculate chromatin region co-accessibility

proj6 <- addCoAccessibility(
    ArchRProj = proj6,
    reducedDims = "LSI_Harmony"
)

proj6 <- addPeak2GeneLinks(
    ArchRProj = proj6,
    reducedDims = "LSI_Harmony",
    useMatrix = "GeneIntegrationMatrix_LSI_Harmony"
)

cA <- getCoAccessibility(
    ArchRProj = proj6,
    corCutOff = 0.5,
    resolution = 10000,
    returnLoops = TRUE
)

15b. Calculate marker features

# Exclude populations with <10 cells

markersPeaks <- getMarkerFeatures(
    ArchRProj = proj6,
    useMatrix = "PeakMatrix", 
    groupBy = "ManualClustersJG",
    bias = c("TSSEnrichment", "log10(nFrags)"),
    testMethod = "wilcoxon",
    threads=4
)

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

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

15c. Export marker features

• Isolate chromatin regions unique to JG cells.

lateJG_markerPeaks <- markerList[["late_JGs"]]
lateJG_markerPeaksGR <- makeGRangesFromDataFrame(lateJG_markerPeaks,
                                                 keep.extra.columns=T)
lateJG_markerPeaksGR <- sort(lateJG_markerPeaksGR)

earlyJG_markerPeaks <- markerList[["early_JGs"]]
earlyJG_markerPeaksGR <- makeGRangesFromDataFrame(earlyJG_markerPeaks,
                                                  keep.extra.columns=T)
earlyJG_markerPeaksGR <- sort(earlyJG_markerPeaksGR)

JG_markerPeaksGR <- subsetByOverlaps(earlyJG_markerPeaksGR,
                                     lateJG_markerPeaksGR)
export.bed(JG_markerPeaksGR, con='JG_peakset.bed')

• Export peaks at the renin locus for JG cells populations.

peakSet <- getPeakSet(proj6)
export.bed(peakSet, con='consensus_peakset.bed')

ren1_wide <- data.table(chr="chr1",
                        start="133300000",
                        end="133400000",
                        name="REN1 wide",
                        strand="*",
                        score="-")
ren1_wide <- makeGRangesFromDataFrame(ren1_wide, keep.extra.columns=T)

ren1_region_peaks_in_lateJGs <- subsetByOverlaps(lateJG_markerPeaksGR, ren1_wide)
data.table::fwrite(data.table::as.data.table(ren1_region_peaks_in_lateJGs),
                   "lateJG_marker_peaks_in_ren1-wide.csv", sep=",")

ren1_region_peaks_in_earlyJGs <- subsetByOverlaps(earlyJG_markerPeaksGR,
                                                  ren1_wide)
data.table::fwrite(data.table::as.data.table(ren1_region_peaks_in_earlyJGs),
                   "earlyJG_marker_peaks_in_ren1-wide.csv", sep=",")

ren1_overlapping_peaks <- sort(subsetByOverlaps(peakSet, ren1_wide))
export.bed(ren1_overlapping_peaks, con='ren1_peakset.bed')

JG_ren1_peakSet <- GenomicRanges::union(
    ren1_overlapping_peaks[names(ren1_overlapping_peaks)=="early_JGs"],
    ren1_overlapping_peaks[names(ren1_overlapping_peaks)=="late_JGs"])

export.bed(JG_ren1_peakSet, con='JG_ren1_peakset.bed')

• Export peaks at the Akr1b7 locus for JG cells populations.

akr1b7_wide <- data.table(chr="chr6",
                          start="34339269",
                          end="34450077",
                          name="Akr1b7 wide",
                          strand="*",
                          score="-")
akr1b7_wide <- makeGRangesFromDataFrame(akr1b7_wide, keep.extra.columns=T)

akr1b7_region_peaks_in_lateJGs <- subsetByOverlaps(lateJG_markerPeaksGR, akr1b7_wide)
data.table::fwrite(data.table::as.data.table(akr1b7_region_peaks_in_lateJGs),
                   "lateJG_marker_peaks_in_akr1b7-wide.csv", sep=",")

akr1b7_region_peaks_in_earlyJGs <- subsetByOverlaps(earlyJG_markerPeaksGR,
                                                  akr1b7_wide)
data.table::fwrite(data.table::as.data.table(akr1b7_region_peaks_in_earlyJGs),
                   "earlyJG_marker_peaks_in_akr1b7-wide.csv", sep=",")

akr1b7_overlapping_peaks <- sort(subsetByOverlaps(peakSet, akr1b7_wide))
export.bed(akr1b7_overlapping_peaks, con='akr1b7_peakset.bed')

JG_Akr1b7_peakSet <- GenomicRanges::union(
    akr1b7_overlapping_peaks[names(akr1b7_overlapping_peaks)=="early_JGs"],
    akr1b7_overlapping_peaks[names(akr1b7_overlapping_peaks)=="late_JGs"])

export.bed(JG_Akr1b7_peakSet, con='JG_akr1b7_peakset.bed')

15d. Grab FASTA sequence of peak sets

Take marker peaks and use them as the input to FIMO.

# FASTA of all JG marker peaks
bedtools getfasta -fi /project/shefflab/genomes_v04_210301/alias/mm10/fasta/default/mm10.fa -bed JG_peakset.bed > JG_peakset.fa

# FASTA of JG marker peaks that overlap the Akr1b7 locus
bedtools getfasta -fi /project/shefflab/genomes_v04_210301/alias/mm10/fasta/default/mm10.fa -bed JG_akr1b7_peakset.bed > JG_akr1b7_peakset.fa

# FASTA of all marker peaks that overlap the Akr1b7 locus
bedtools getfasta -fi /project/shefflab/genomes_v04_210301/alias/mm10/fasta/default/mm10.fa -bed akr1b7_peakset.bed > akr1b7_peakset.fa

# FASTA of JG marker peaks that overlap the REN1 locus
bedtools getfasta -fi /project/shefflab/genomes_v04_210301/alias/mm10/fasta/default/mm10.fa -bed JG_ren1_peakset.bed > JG_ren1_peakset.fa

# FASTA of all marker peaks that overlap the REN1 locus
bedtools getfasta -fi /project/shefflab/genomes_v04_210301/alias/mm10/fasta/default/mm10.fa -bed ren1_peakset.bed > ren1_peakset.fa

# FASTA of the consensus peak set
bedtools getfasta -fi /project/shefflab/genomes_v04_210301/alias/mm10/fasta/default/mm10.fa -bed consensus_peakset.bed > consensus_peakset.fa

15e. Prepare background files

module load gcc/9.2.0  openmpi/3.1.6 meme/5.3.3
export PATH="$PATH:/apps/software/standard/mpi/gcc/9.2.0/openmpi/3.1.6/meme/5.3.3/libexec/meme-5.3.3/"

# Use consensus peak set as background
fasta-get-markov -dna consensus_peakset.fa consensus_peakset.bg

15f. Extract motif set used in ArchR

library(chromVARmotifs)
library(dplyr)
library(tidyr)
library(data.table)

# Generate the PWM used from chromVAR
data("mouse_pwms_v1")
sout <- sapply(strsplit(names(mouse_pwms_v1), split = "_"), function(s) c(s[3]))
mouse_pwms_v2 <- mouse_pwms_v1[match(unique(sout), sout)]

# Convert the PWMS to MEME format for FIMO
dir.create(file.path(getwd(), "cisbp_mm_meme"))
meme_out <- file.path(getwd(), "cisbp_mm_meme/mouse_pwms_v2.meme")
write_meme(convert_motifs(mouse_pwms_v2), meme_out)

15g. Identify motifs of interest.

Grab top 20 unique motifs enriched in several pairwise comparisons for the JG populations.

15h. Run FIMO on the JG enriched motifs

Identify motif binding site occurrences at the Ren1 locus.

cd /sfs/lustre/bahamut/scratch/jps3dp/data/alexandre/ATACseq/processed/scATAC_analysis/R_analysis/sc_kidney/
module load gcc/9.2.0  openmpi/3.1.6 meme/5.3.3

export CISBPDIR=/sfs/lustre/bahamut/scratch/jps3dp/data/alexandre/ATACseq/processed/scATAC_analysis/R_analysis/sc_kidney/cisbp_mm_meme/

mkdir FIMO/

# JG only peaks around Ren1
fimo --motif Nfix --motif Nfic --motif Rfx1 --motif Rfx2 --motif Rfx4 --motif Rfx8 --motif Smarcc1 --motif Fli1 --motif Fos --motif Ebf1 --motif Bach1 --motif Bach2 --motif Creb1 --motif Creb3 --motif Creb3l3 --motif Elf1 --motif Elf3 --motif Elf5 --motif Erf --motif Erg --motif Ets1 --motif Etv2 --motif LINE5773 --motif LINE5785 --motif LINE5798 --motif LINE5816 --motif LINE5827 --motif LINE5828 --motif Maz --motif Mef2a --motif Mef2b --motif Mef2c --motif Mef2d --motif Nfia --motif Nfib --motif Pgr --motif Rarb --motif Srf --motif Stat3 --motif Thra --motif Thrb --motif Zfa --motif Zfy1 --motif Rbpj --bfile consensus_peakset.bg -oc FIMO/chromVAR_PWM_cisbp_JG_peaks/ ${CISBPDIR}mouse_pwms_v2.meme JG_ren1_peakset.fa

# All peaks around Ren1
fimo --motif Nfix --motif Nfic --motif Rfx1 --motif Rfx2 --motif Rfx4 --motif Rfx8 --motif Smarcc1 --motif Fli1 --motif Fos --motif Ebf1 --motif Bach1 --motif Bach2 --motif Creb1 --motif Creb3 --motif Creb3l3 --motif Elf1 --motif Elf3 --motif Elf5 --motif Erf --motif Erg --motif Ets1 --motif Etv2 --motif LINE5773 --motif LINE5785 --motif LINE5798 --motif LINE5816 --motif LINE5827 --motif LINE5828 --motif Maz --motif Mef2a --motif Mef2b --motif Mef2c --motif Mef2d --motif Nfia --motif Nfib --motif Pgr --motif Rarb --motif Srf --motif Stat3 --motif Thra --motif Thrb --motif Zfa --motif Zfy1 --motif Rbpj --bfile consensus_peakset.bg -oc FIMO/chromVAR_PWM_cisbp_ren1_all_marker_peaks/ ${CISBPDIR}mouse_pwms_v2.meme ren1_peakset.fa

Identify motif binding site occurrences at the Akr1b7 locus.

cd /sfs/lustre/bahamut/scratch/jps3dp/data/alexandre/ATACseq/processed/scATAC_analysis/R_analysis/sc_kidney/
module load gcc/9.2.0  openmpi/3.1.6 meme/5.3.3

export CISBPDIR=/sfs/lustre/bahamut/scratch/jps3dp/data/alexandre/ATACseq/processed/scATAC_analysis/R_analysis/sc_kidney/cisbp_mm_meme/

mkdir FIMO/

# JG peaks only around Akr1b7
fimo --motif Nfix --motif Nfic --motif Rfx1 --motif Rfx2 --motif Rfx4 --motif Rfx8 --motif Smarcc1 --motif Fli1 --motif Fos --motif Ebf1 --motif Bach1 --motif Bach2 --motif Creb1 --motif Creb3 --motif Creb3l3 --motif Elf1 --motif Elf3 --motif Elf5 --motif Erf --motif Erg --motif Ets1 --motif Etv2 --motif LINE5773 --motif LINE5785 --motif LINE5798 --motif LINE5816 --motif LINE5827 --motif LINE5828 --motif Maz --motif Mef2a --motif Mef2b --motif Mef2c --motif Mef2d --motif Nfia --motif Nfib --motif Pgr --motif Rarb --motif Srf --motif Stat3 --motif Thra --motif Thrb --motif Zfa --motif Zfy1 --motif Rbpj --bfile consensus_peakset.bg -oc FIMO/chromVAR_PWM_cisbp_akr1b7_JG_peaks/ ${CISBPDIR}mouse_pwms_v2.meme JG_akr1b7_peakset.fa

# All peaks around Akr1b7
fimo --motif Nfix --motif Nfic --motif Rfx1 --motif Rfx2 --motif Rfx4 --motif Rfx8 --motif Smarcc1 --motif Fli1 --motif Fos --motif Ebf1 --motif Bach1 --motif Bach2 --motif Creb1 --motif Creb3 --motif Creb3l3 --motif Elf1 --motif Elf3 --motif Elf5 --motif Erf --motif Erg --motif Ets1 --motif Etv2 --motif LINE5773 --motif LINE5785 --motif LINE5798 --motif LINE5816 --motif LINE5827 --motif LINE5828 --motif Maz --motif Mef2a --motif Mef2b --motif Mef2c --motif Mef2d --motif Nfia --motif Nfib --motif Pgr --motif Rarb --motif Srf --motif Stat3 --motif Thra --motif Thrb --motif Zfa --motif Zfy1 --motif Rbpj --bfile consensus_peakset.bg -oc FIMO/chromVAR_PWM_cisbp_akr1b7_all_marker_peaks/ ${CISBPDIR}mouse_pwms_v2.meme akr1b7_peakset.fa

15i. Run FIMO on literature reported Ren1 promoter binding TFs

fimo --motif Sp1 --motif Sp3 --motif Nfix --motif Pbx1 --motif Hoxa9 --motif Hoxa10 --motif Pknox1 --bfile consensus_peakset.bg -oc FIMO/chromVAR_PWM_cisbp_JGpeaks_promoter/ ${CISBPDIR}mouse_pwms_v2.meme JG_ren1_peakset.fa

fimo --motif Sp1 --motif Sp3 --motif Nfix --motif Pbx1 --motif Hoxa9 --motif Hoxa10 --motif Pknox1 --bfile consensus_peakset.bg -oc FIMO/chromVAR_PWM_cisbp_promoter/ ${CISBPDIR}mouse_pwms_v2.meme ren1_peakset.fa

15j. Run FIMO on literature reported Ren1 enhancer binding TFs

fimo --motif Sp1 --motif Sp3 --motif Nfix --motif Rxrb --motif Rara --motif Rxra --motif Rxrg --motif Nfya --motif Crem --motif Creb1 --motif Usf1 --motif Usf2 --bfile consensus_peakset.bg -oc FIMO/chromVAR_PWM_cisbp_JGpeaks_enhancer/ ${CISBPDIR}mouse_pwms_v2.meme JG_ren1_peakset.fa

fimo --motif Sp1 --motif Sp3 --motif Nfix --motif Rxrb --motif Rara --motif Rxra --motif Rxrg --motif Nfya --motif Crem --motif Creb1 --motif Usf1 --motif Usf2 --bfile consensus_peakset.bg -oc FIMO/chromVAR_PWM_cisbp_enhancer/ ${CISBPDIR}mouse_pwms_v2.meme ren1_peakset.fa