update

54d1f5c 5 months ago

5.67 kB

import argparse
import multiprocessing as mp
import os
import sys
import traceback

import loompy
import numpy as np
import pandas as pd
import scanpy as sc
import scrublet
from loguru import logger


def get_filename() -> str:
    current_file_name = os.path.basename(__file__)
    log_name = "{}.log"
    return log_name.format(current_file_name.split('.')[0])


def set_log() -> None:
    filename = get_filename()
    logger.add(f'../log/{filename}')


set_log()


gene_info_table_path = '/scDifformer/data/240412_test/material/gene_info_table.csv'
out_path = '/scDifformer/data/240412_test/looms'
home_dir_list = ['/scDifformer/data/240412_test/raw_h5ad']

gene_info = pd.read_csv(gene_info_table_path)

# Convert to dictionaries
gene_name_id_combine_dict = gene_info.set_index('gene_name')['ensembl_id'].to_dict()
gene_name_type_dict = gene_info.set_index('gene_name')['gene_type'].to_dict()
gene_id_type_dict = gene_info.set_index('ensembl_id')['gene_type'].to_dict()
func_gene_list = [i for i in
                  gene_info[(gene_info["gene_type"] == "protein_coding") | (gene_info["gene_type"] == "miRNA")][
                      "ensembl_id"]]


def preprocess(file_path: str) -> None:
    try:
        file_name = file_path.split(os.sep)[-1].split(".")[0]
        adata = sc.read_h5ad(file_path)
        adata.raw = None
        adata.X = np.round(adata.X, decimals=1)
        # 1. Filter mitochondrial ratio and total gene counts within 3 standard deviations of the mean
        adata.var['mt'] = adata.var_names.str.startswith('MT-')  # annotate the group of mitochondrial genes as 'mt'
        sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], percent_top=None, log1p=False, inplace=True)

        # Compute mean and standard deviation of mitochondrial ratio
        mean_ratio_mt = np.mean(adata.obs["pct_counts_mt"])
        std_ratio_mt = np.std(adata.obs["pct_counts_mt"])

        # Keep cells within three standard deviations of the mean
        min_ratio_mt = mean_ratio_mt - 3 * std_ratio_mt
        max_ratio_mt = mean_ratio_mt + 3 * std_ratio_mt
        within_range_mt = (adata.obs["pct_counts_mt"] >= min_ratio_mt) & (adata.obs["pct_counts_mt"] <= max_ratio_mt)
        # Remove cells outside the valid range
        adata = adata[within_range_mt, :]

        # Compute mean and standard deviation of total gene counts
        mean_ratio = np.mean(adata.obs["total_counts"])
        std_ratio = np.std(adata.obs["total_counts"])

        # Keep cells within three standard deviations of the mean
        min_ratio = mean_ratio - 3 * std_ratio
        max_ratio = mean_ratio + 3 * std_ratio
        within_range = (adata.obs["total_counts"] >= min_ratio) & (adata.obs["total_counts"] <= max_ratio)

        # Remove cells outside the valid range
        adata = adata[within_range]

        # 2. Create a DoubletFinder-like object and run doublet detection
        scrub = scrublet.Scrublet(adata.X)
        # Run doublet detection
        doublet_scores, predicted_doublets = scrub.scrub_doublets()
        # Filter the cells based on the doublet scores and predicted doublets
        if predicted_doublets is not None:
            adata = adata[~predicted_doublets]

        # 3. exclude cells with less than seven detected Ensembl-annotated protein-coding or miRNA genes
        retain_gene_set = set[str]()
        func_gene_set = set(func_gene_list)
        for gene_name in adata.var_names:
            if gene_name in func_gene_set:
                retain_gene_set.add(gene_name)
            elif gene_name in gene_name_id_combine_dict and gene_name_id_combine_dict[gene_name] in func_gene_set:
                retain_gene_set.add(gene_name)

        adata = adata[:, list(retain_gene_set)]
        sc.pp.filter_cells(adata, min_genes=7)
        gene_type_arr = []
        ensembl_id_arr = []
        for var_info in adata.var_names:
            if var_info in gene_name_type_dict:
                gene_type_arr.append(gene_name_type_dict[var_info])
                ensembl_id_arr.append(gene_name_id_combine_dict[var_info])
            elif str(var_info).startswith("ENSG"):
                gene_type_arr.append(gene_id_type_dict[var_info])
                ensembl_id_arr.append(var_info)
            else:
                raise Exception(f"can not find {var_info}")

        total_genes_per_row = adata.X.sum(axis=1)
        filename = f'{out_path}/{file_name}.loom'
        loompy.create(filename, adata.X.T, {"ensembl_id": np.array(ensembl_id_arr), "gene_type": gene_type_arr},
                      {"n_counts": total_genes_per_row})
        logger.info(f"{filename} is done")
        del adata
    except Exception as e:
        logger.error(f"====={file_path}")
        logger.error(f'{traceback.print_exc()}')


def argv_to_args(argv):
    """Parse argument from cmd"""
    parser = argparse.ArgumentParser()

    # start
    parser.add_argument("--start", '-s', type=int, help="start num for iter", required=True)
    parser.add_argument("--end", '-e', type=int, help="end num for iter", default=1)

    try:
        args = parser.parse_args(argv[1:])
    except Exception as e:
        logger.error(e)
        parser.print_help()
        sys.exit(1)

    return args


no_need_handle_set = set()

if __name__ == '__main__':
    for home_dir in home_dir_list:
        home_path = home_dir
        files = [os.path.join(home_path, file_name) for file_name in os.listdir(home_path) if
                 str(file_name).endswith(".h5ad") and file_name not in no_need_handle_set]
        logger.info(f'Number of files processed in this run: {len(files)}')
        with mp.Pool(processes=32) as pool:
            pool.map(preprocess, files)

	import argparse
	import multiprocessing as mp
	import os
	import sys
	import traceback

	import loompy
	import numpy as np
	import pandas as pd
	import scanpy as sc
	import scrublet
	from loguru import logger


	def get_filename() -> str:
	current_file_name = os.path.basename(__file__)
	log_name = "{}.log"
	return log_name.format(current_file_name.split('.')[0])


	def set_log() -> None:
	filename = get_filename()
	logger.add(f'../log/{filename}')


	set_log()


	gene_info_table_path = '/scDifformer/data/240412_test/material/gene_info_table.csv'
	out_path = '/scDifformer/data/240412_test/looms'
	home_dir_list = ['/scDifformer/data/240412_test/raw_h5ad']

	gene_info = pd.read_csv(gene_info_table_path)

	# Convert to dictionaries
	gene_name_id_combine_dict = gene_info.set_index('gene_name')['ensembl_id'].to_dict()
	gene_name_type_dict = gene_info.set_index('gene_name')['gene_type'].to_dict()
	gene_id_type_dict = gene_info.set_index('ensembl_id')['gene_type'].to_dict()
	func_gene_list = [i for i in
	gene_info[(gene_info["gene_type"] == "protein_coding") \| (gene_info["gene_type"] == "miRNA")][
	"ensembl_id"]]


	def preprocess(file_path: str) -> None:
	try:
	file_name = file_path.split(os.sep)[-1].split(".")[0]
	adata = sc.read_h5ad(file_path)
	adata.raw = None
	adata.X = np.round(adata.X, decimals=1)
	# 1. Filter mitochondrial ratio and total gene counts within 3 standard deviations of the mean
	adata.var['mt'] = adata.var_names.str.startswith('MT-') # annotate the group of mitochondrial genes as 'mt'
	sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], percent_top=None, log1p=False, inplace=True)

	# Compute mean and standard deviation of mitochondrial ratio
	mean_ratio_mt = np.mean(adata.obs["pct_counts_mt"])
	std_ratio_mt = np.std(adata.obs["pct_counts_mt"])

	# Keep cells within three standard deviations of the mean
	min_ratio_mt = mean_ratio_mt - 3 * std_ratio_mt
	max_ratio_mt = mean_ratio_mt + 3 * std_ratio_mt
	within_range_mt = (adata.obs["pct_counts_mt"] >= min_ratio_mt) & (adata.obs["pct_counts_mt"] <= max_ratio_mt)
	# Remove cells outside the valid range
	adata = adata[within_range_mt, :]

	# Compute mean and standard deviation of total gene counts
	mean_ratio = np.mean(adata.obs["total_counts"])
	std_ratio = np.std(adata.obs["total_counts"])

	# Keep cells within three standard deviations of the mean
	min_ratio = mean_ratio - 3 * std_ratio
	max_ratio = mean_ratio + 3 * std_ratio
	within_range = (adata.obs["total_counts"] >= min_ratio) & (adata.obs["total_counts"] <= max_ratio)

	# Remove cells outside the valid range
	adata = adata[within_range]

	# 2. Create a DoubletFinder-like object and run doublet detection
	scrub = scrublet.Scrublet(adata.X)
	# Run doublet detection
	doublet_scores, predicted_doublets = scrub.scrub_doublets()
	# Filter the cells based on the doublet scores and predicted doublets
	if predicted_doublets is not None:
	adata = adata[~predicted_doublets]

	# 3. exclude cells with less than seven detected Ensembl-annotated protein-coding or miRNA genes
	retain_gene_set = set[str]()
	func_gene_set = set(func_gene_list)
	for gene_name in adata.var_names:
	if gene_name in func_gene_set:
	retain_gene_set.add(gene_name)
	elif gene_name in gene_name_id_combine_dict and gene_name_id_combine_dict[gene_name] in func_gene_set:
	retain_gene_set.add(gene_name)

	adata = adata[:, list(retain_gene_set)]
	sc.pp.filter_cells(adata, min_genes=7)
	gene_type_arr = []
	ensembl_id_arr = []
	for var_info in adata.var_names:
	if var_info in gene_name_type_dict:
	gene_type_arr.append(gene_name_type_dict[var_info])
	ensembl_id_arr.append(gene_name_id_combine_dict[var_info])
	elif str(var_info).startswith("ENSG"):
	gene_type_arr.append(gene_id_type_dict[var_info])
	ensembl_id_arr.append(var_info)
	else:
	raise Exception(f"can not find {var_info}")

	total_genes_per_row = adata.X.sum(axis=1)
	filename = f'{out_path}/{file_name}.loom'
	loompy.create(filename, adata.X.T, {"ensembl_id": np.array(ensembl_id_arr), "gene_type": gene_type_arr},
	{"n_counts": total_genes_per_row})
	logger.info(f"{filename} is done")
	del adata
	except Exception as e:
	logger.error(f"====={file_path}")
	logger.error(f'{traceback.print_exc()}')


	def argv_to_args(argv):
	"""Parse argument from cmd"""
	parser = argparse.ArgumentParser()

	# start
	parser.add_argument("--start", '-s', type=int, help="start num for iter", required=True)
	parser.add_argument("--end", '-e', type=int, help="end num for iter", default=1)

	try:
	args = parser.parse_args(argv[1:])
	except Exception as e:
	logger.error(e)
	parser.print_help()
	sys.exit(1)

	return args


	no_need_handle_set = set()

	if __name__ == '__main__':
	for home_dir in home_dir_list:
	home_path = home_dir
	files = [os.path.join(home_path, file_name) for file_name in os.listdir(home_path) if
	str(file_name).endswith(".h5ad") and file_name not in no_need_handle_set]
	logger.info(f'Number of files processed in this run: {len(files)}')
	with mp.Pool(processes=32) as pool:
	pool.map(preprocess, files)