Pseudobulk profiles not being generated.

[patricksis]

Hello!

I think this is the right place to post this question. I'm somewhat new to this, so any suggestions will be helpful!

I am trying to reproduce results from scenicplus tutorial here:
https://scenicplus.readthedocs.io/en/latest/pbmc_multiome_tutorial.html#Tutorial:-10x-multiome-pbmc

However I am encountering an issue when it comes to generating pseudobulk profiles.

My output looks like this:

2023-08-30 13:47:04,741 cisTopic     INFO     Reading fragments from /home/psisler/Projects/scenic/analysis/pbmc/data/pbmc_granulocyte_sorted_3k_atac_fragments.tsv.gz
2023-08-30 13:47:29,373 INFO worker.py:1621 -- Started a local Ray instance.
(export_pseudobulk_ray pid=2891529) 2023-08-30 13:47:32,349 cisTopic     INFO     Creating pseudobulk for B_cells
(export_pseudobulk_ray pid=2891526) 2023-08-30 13:47:32,655 cisTopic     INFO     Creating pseudobulk for B_cells_1
(export_pseudobulk_ray pid=2891525) 2023-08-30 13:47:33,868 cisTopic     INFO     Creating pseudobulk for NK_cells
(export_pseudobulk_ray pid=2891531) 2023-08-30 13:47:34,390 cisTopic     INFO     Creating pseudobulk for Dendritic_cells
(export_pseudobulk_ray pid=2891527) 2023-08-30 13:47:34,752 cisTopic     INFO     Creating pseudobulk for FCGR3A_Monocytes [repeated 5x across cluster] (Ray deduplicates logs by default. Set RAY_DEDUP_LOGS=0 to disable log deduplication, or see https://docs.ray.io/en/master/ray-observability/ray-logging.html#log-deduplication for more options.)

According to the tutorial is should look like this:

(export_pseudobulk_ray pid=13683) 2022-08-09 17:57:43,021 cisTopic     INFO     Creating pseudobulk for B_cells_1
(export_pseudobulk_ray pid=13679) 2022-08-09 17:57:43,404 cisTopic     INFO     Creating pseudobulk for B_cells_2
(export_pseudobulk_ray pid=13677) 2022-08-09 17:57:43,888 cisTopic     INFO     Creating pseudobulk for CD14_Monocytes
(export_pseudobulk_ray pid=13681) 2022-08-09 17:57:44,305 cisTopic     INFO     Creating pseudobulk for CD4_T_cells
(export_pseudobulk_ray pid=13682) 2022-08-09 17:57:44,706 cisTopic     INFO     Creating pseudobulk for CD8_T_cells
(export_pseudobulk_ray pid=13680) 2022-08-09 17:57:45,230 cisTopic     INFO     Creating pseudobulk for Dendritic_cells
(export_pseudobulk_ray pid=13684) 2022-08-09 17:57:45,726 cisTopic     INFO     Creating pseudobulk for FCGR3A_Monocytes
(export_pseudobulk_ray pid=13678) 2022-08-09 17:57:46,225 cisTopic     INFO     Creating pseudobulk for NK_cells
(export_pseudobulk_ray pid=13679) 2022-08-09 17:57:52,695 cisTopic     INFO     B_cells_2 done!
(export_pseudobulk_ray pid=13680) 2022-08-09 17:58:00,289 cisTopic     INFO     Dendritic_cells done!
(export_pseudobulk_ray pid=13678) 2022-08-09 17:58:01,561 cisTopic     INFO     NK_cells done!
(export_pseudobulk_ray pid=13684) 2022-08-09 17:58:08,569 cisTopic     INFO     FCGR3A_Monocytes done!
(export_pseudobulk_ray pid=13682) 2022-08-09 17:58:24,903 cisTopic     INFO     CD8_T_cells done!
(export_pseudobulk_ray pid=13683) 2022-08-09 17:58:29,227 cisTopic     INFO     B_cells_1 done!
(export_pseudobulk_ray pid=13677) 2022-08-09 18:00:49,939 cisTopic     INFO     CD14_Monocytes done!

Here is the relevant code

import os
import pycisTopic
import scanpy as sc
import pyranges as pr
import requests
import pandas as pd
import pickle

work_dir = '/home/psisler/Projects/scenic/analysis/pbmc'
tmp_dir = '/home/psisler/Projects/scenic/tmp'
#set some figure parameters for nice display inside jupyternotebooks.
#%matplotlib inline

#make a directory for to store the processed scRNA-seq data.
if not os.path.exists(os.path.join(work_dir, 'scATAC')):
        os.makedirs(os.path.join(work_dir, 'scATAC'))
        tmp_dir = '/home/psisler/Projects/scenic/tmp'

fragments_dict = {'10x_pbmc': os.path.join(work_dir, 'data/pbmc_granulocyte_sorted_3k_atac_fragments.tsv.gz')}

adata = sc.read_h5ad(os.path.join(work_dir, 'scRNA/adata.h5ad'))
cell_data = adata.obs
cell_data['sample_id'] = '10x_pbmc'
cell_data['celltype'] = cell_data['celltype'].astype(str) # set data type of the celltype column to str, otherwise the export_pseudobulk function will complain.
del(adata)

#Get chromosome sizes (for hg38 here)
target_url='http://hgdownload.cse.ucsc.edu/goldenPath/hg38/bigZips/hg38.chrom.sizes'
chromsizes=pd.read_csv(target_url, sep='\t', header=None)
chromsizes.columns=['Chromosome', 'End']
chromsizes['Start']=[0]*chromsizes.shape[0]
chromsizes=chromsizes.loc[:,['Chromosome', 'Start', 'End']]
# Exceptionally in this case, to agree with CellRangerARC annotations
chromsizes['Chromosome'] = [chromsizes['Chromosome'][x].replace('v', '.') for x in range(len(chromsizes['Chromosome']))]
chromsizes['Chromosome'] = [chromsizes['Chromosome'][x].split('_')[1] if len(chromsizes['Chromosome'][x].split('_')) > 1 else chromsizes['Chromosome'][x] for x in range(len(chromsizes['Chromosome']))]
chromsizes=pr.PyRanges(chromsizes)
from pycisTopic.pseudobulk_peak_calling import export_pseudobulk

bw_paths, bed_paths = export_pseudobulk(input_data = cell_data,
                variable = 'celltype',
                sample_id_col = 'sample_id',
                chromsizes = chromsizes,
                bed_path = os.path.join(work_dir, 'scATAC/consensus_peak_calling/pseudobulk_bed_files/'),
                bigwig_path = os.path.join(work_dir, 'scATAC/consensus_peak_calling/pseudobulk_bw_files/'),
                path_to_fragments = fragments_dict,
                n_cpu = 8,
                normalize_bigwig = True,
                remove_duplicates = True,
                _temp_dir = os.path.join(tmp_dir, 'ray_spill'),
                split_pattern = '-')

I saw a few issue related to this(#61) but it seems to be related to a different dataset.

Would anybody have any ideas on how to resolve this issue? I have an idea that it might be a dependency issue, but I am not sure.

Thank you

[patricksis]

Ah it seems I resolved the issue. If anyone has a similar issue please go here:
#138 (comment)