Differential motif enrichment

[CarstenKnutsen]

I was wondering on everyone's opinion on an appropriate workflow for testing for differential TF motif enrichment.

In a scenario where you are comparing a healthy and diseased cell type, what do you think is the most appropriate way for finding differential TF motif enrichment?

Below is a snippet of code outlining how I have planned it out. Essentially, taking the differentially accessible peaks between treatments, then running the motif enrichment on the significantly different peaks, using all peaks detected in that cell type as a background signal.

b_cell = peak_mat[(peak_mat.obs['celltype']=='B cell')]
df = snap.tl.diff_test(peak_mat,
                  b_cell[b_cell.obs['treatment']=='Healthy'].obs_names,
                  b_cell[b_cell.obs['treatment']=='Disease'].obs_names,
                 )
df_pd = df.to_pandas()
df_pd_t = df_pd.loc[df_pd['adjusted p-value']<0.05]
up = df_pd_t[df_pd_t['log2(fold_change)']>0]['feature name'].values
dn = df_pd_t[df_pd_t['log2(fold_change)']<0]['feature name'].values
dict = {'Healthy':up,
       'Diseased':dn,
       }
peaks = peaks.to_pandas()
b_peaks = peaks[peaks['B cell']==True]
motifs = snap.tl.motif_enrichment(
    motifs=snap.datasets.Meuleman_2020(),
    regions=dict,
    background=b_peaks['Peaks'].values,
    genome_fasta=genome,
    method='hypergeometric'
)

Do you think this workflow is appropriate for finding different TF binding patterns between treatments?

[kaizhang]

The analysis you have done looks good to me.

[CarstenKnutsen]

Thanks so much for the response! I appreciate it.

What is the best way within SnapATAC2 to annotate our peaks with genomic information, similar to how Cell Ranger ARC annotates peaks. When going through these lists, it is nice to be able to see what genes they are nearest.

Should I go through SnapATAC2's network analysis?
Or would using something like HOMER where I just use genomic location to annotate the peaks make more sense?

[kaizhang]

This is indeed a very useful feature. This is on my to-do list and will be added to SnapATAC2 in the future. Network analysis can also give you such information but is something more advanced. Please keep in mind that the network analysis is still in beta.

[CarstenKnutsen]

Sounds good, thanks for the information. Been loving the package so far!

[yojetsharma]

The analysis you have done looks good to me.

I also got inspired by the author's post and did the following, can you please tell me if it is right? I basically want to see the enrichment in the patient samples wrt to the control. Had carried out multiome (10x) for these samples:
(Initial steps are from the multisample tutorial)

snap.tl.macs3(data, groupby='sample') #instead of grouping by 'leiden' I chose to do it as per sample
peaks = snap.tl.merge_peaks(data.uns['macs3'], snap.genome.hg38)
peaks.head()
peak_mat = snap.pp.make_peak_matrix(data, use_rep=peaks['Peaks'])
peak_mat
healthy_cells = peak_mat[peak_mat.obs['treatment'] == 'Healthy']
diseased_cells = peak_mat[peak_mat.obs['treatment'] == 'Disease']
df = snap.tl.diff_test(
    peak_mat,
    healthy_cells.obs_names,
    diseased_cells.obs_names,
)
df_pd = df.to_pandas()
df_pd_t = df_pd.loc[df_pd['adjusted p-value'] < 0.05]

up = df_pd_t[df_pd_t['log2(fold_change)'] > 0]['feature name'].values
dn = df_pd_t[df_pd_t['log2(fold_change)'] < 0]['feature name'].values

dict = {
    'Healthy': up,
    'Diseased': dn,
}
peaks = peaks.to_pandas()
motifs = snap.tl.motif_enrichment(
    motifs=snap.datasets.cis_bp(unique=True),
    regions=dict,
    background=peaks['Peaks'].values,  
    genome_fasta=snap.genome.hg38,
)
motifs_df = motifs.to_pandas()
significant_motifs = motifs_df[motifs_df['p_value'] < 0.05]
top_motifs = significant_motifs.sort_values(by='enrichment_score', ascending=False)

[yojetsharma]

I was wondering on everyone's opinion on an appropriate workflow for testing for differential TF motif enrichment.

In a scenario where you are comparing a healthy and diseased cell type, what do you think is the most appropriate way for finding differential TF motif enrichment?

Below is a snippet of code outlining how I have planned it out. Essentially, taking the differentially accessible peaks between treatments, then running the motif enrichment on the significantly different peaks, using all peaks detected in that cell type as a background signal.

b_cell = peak_mat[(peak_mat.obs['celltype']=='B cell')]
df = snap.tl.diff_test(peak_mat,
                  b_cell[b_cell.obs['treatment']=='Healthy'].obs_names,
                  b_cell[b_cell.obs['treatment']=='Disease'].obs_names,
                 )
df_pd = df.to_pandas()
df_pd_t = df_pd.loc[df_pd['adjusted p-value']<0.05]
up = df_pd_t[df_pd_t['log2(fold_change)']>0]['feature name'].values
dn = df_pd_t[df_pd_t['log2(fold_change)']<0]['feature name'].values
dict = {'Healthy':up,
       'Diseased':dn,
       }
peaks = peaks.to_pandas()
b_peaks = peaks[peaks['B cell']==True]
motifs = snap.tl.motif_enrichment(
    motifs=snap.datasets.Meuleman_2020(),
    regions=dict,
    background=b_peaks['Peaks'].values,
    genome_fasta=genome,
    method='hypergeometric'
)

Do you think this workflow is appropriate for finding different TF binding patterns between treatments?

I have a basic query. How can I load multiple atac frag files to do this analyses? I have a WT and diseased too. It would be really helpful if you could share your code?