Integration with SCENICplus

[TingTingShao]

Dear,

As I can see from the paper that snapATAC2 can easily integrate with SCENICplus. But I see that in SCENICplus, the first step is to call consensus peaks then do quality control, which is very different from snapATAC2, if I am not wrong.

Could I please ask which step you were referring in the paper to best integrate snapATAC2 with SCENICplus?

Thanks
tingting

[TingTingShao]

It is easy to integrate with pycistpic!

[lf96abc]

Hi @TingTingShao
I'm unsure of how to best prepare the dataset before processing with pycistopic, would you be able to explain/provide code for this please?

I am trying to get the fragment count matrix in the correct format.

Many thanks

[TingTingShao]

Hi,

Here is the code, hope it helps

"""
consensus regions + filtered nuclei from snapATAC to create cisTopic object
"""
rule create_cisobject:
    input:
        frag_paths_pkl=os.path.join(config['work_dir'], "030results_cistopic/frag_paths.pkl"),
        tsv=os.path.join(config['work_dir'], "030results_cistopic/cell_data.tsv"),
        consensus_regions=os.path.join(config['work_dir'], "030results_cistopic/010consensus_peak_calling/consensus_regions.bed")
    output:
        obj_pkl=os.path.join(config['work_dir'], "030results_cistopic/cistopic_obj.pkl")
    run:
        cell_data = pd.read_csv(input.tsv, sep = '\t')
        sample_to_barcodes = cell_data.groupby('sample')['barcode'].agg(list).to_dict()
        print(sample_to_barcodes)
        # cell_data['barcode'] = cell_data['sample']+':'+ cell_data['barcode'] 
        fragments_dict=pickle.load(open(input.frag_paths_pkl, 'rb'))
        unique_samples = set(cell_data['sample'])
        path_to_regions= {sample: input.consensus_regions for sample in unique_samples} 
        cistopic_obj_list=[create_cistopic_object_from_fragments(path_to_fragments=fragments_dict[key],
                                                    path_to_regions=path_to_regions[key],
                                                    path_to_blacklist=config['path_to_blacklist'],
                                                    # metrics=metadata_bc[key],
                                                    valid_bc=sample_to_barcodes[key],
                                                    n_cpu=8,
                                                    # split_pattern= ":",
                                                    project=key) for key in sample_to_barcodes.keys()]   
        print(cistopic_obj_list)  
        cistopic_obj = merge(cistopic_obj_list)
        print(cistopic_obj.cell_data)
        pickle.dump(cistopic_obj,
                    open(output.obj_pkl, 'wb'))  

best,
tingting

[lf96abc]

Hi Tingting,

Thank you for this, I really appreciate it. Could you explain how you generated the consensus_regions.bed please?

Thanks again