diff --git a/config/config.yaml b/config/config.yaml index dded8ec6..066931d3 100644 --- a/config/config.yaml +++ b/config/config.yaml @@ -52,6 +52,10 @@ methods: inferelator: modular: False ext: 500000 + linger: + modular: False + version: 'LINGER' + mode: 'parallel' pando: modular: True exclude_exons: 'TRUE' @@ -108,6 +112,7 @@ colors: granie: '#97de8a' hummus: '#d52727' inferelator: '#ff9796' + linger: '#17a589' pando: '#9467bc' scdori: '#c4afd4' scenicplus: '#8c554a' @@ -194,6 +199,7 @@ method_names: grnboost: GRNBoost2 hummus: HuMMuS inferelator: Inferelator3.0 + linger: LINGER pando: Pando pearson: Pearson random: Random diff --git a/workflow/Snakefile b/workflow/Snakefile index 5a30af0a..f7c980d5 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -97,6 +97,7 @@ include: 'rules/dbs/gst.smk' include: 'rules/dbs/c2g.smk' include: 'rules/dbs/prt.smk' include: 'rules/dbs/sim.smk' +include: 'rules/dbs/linger_prior.smk' # Datasets include: 'rules/dts/pbmc10k.smk' @@ -133,6 +134,7 @@ include: 'rules/mth/scdori.smk' include: 'rules/mth/scenic.smk' include: 'rules/mth/scgpt.smk' include: 'rules/mth/scmtni.smk' +include: 'rules/mth/linger.smk' # Analyses include: 'rules/anl/metrics/prior.smk' diff --git a/workflow/envs/linger.def b/workflow/envs/linger.def new file mode 100644 index 00000000..1ce434ec --- /dev/null +++ b/workflow/envs/linger.def @@ -0,0 +1,90 @@ +Bootstrap: docker +From: ubuntu:20.04 + + +%environment + export PATH=/opt/:$PATH + . "/opt/conda/etc/profile.d/conda.sh" + . "/opt/conda/etc/profile.d/mamba.sh" + conda activate env + +%post + + export DEBIAN_FRONTEND=noninteractive + apt update -y + + # basic packages + apt-get install -y --no-install-recommends tzdata + apt install -y build-essential \ + gcc \ + g++ \ + libstdc++6 \ + cmake \ + wget \ + curl \ + git \ + zlib1g-dev \ + libbz2-dev \ + liblzma-dev \ + libcurl4-openssl-dev \ + libssl-dev \ + libxml2-dev \ + libcairo2-dev \ + libxt-dev \ + libopenblas-dev \ + ca-certificates \ + bedtools + + # conda + wget -O Miniforge3.sh "https://github.com/conda-forge/miniforge/releases/latest/download/Miniforge3-$(uname)-$(uname -m).sh" + bash Miniforge3.sh -b -p "/opt/conda" + rm Miniforge3.sh + . "/opt/conda/etc/profile.d/conda.sh" + . "/opt/conda/etc/profile.d/mamba.sh" + conda activate + + # Create env + mamba create -y -n env -c pytorch -c bioconda -c conda-forge \ + python=3.10 \ + pip \ + setuptools \ + wheel \ + cython \ + pytorch=2.5.1 \ + torchvision \ + torchaudio \ + cpuonly \ + pybedtools==0.10.0 \ + bedtools \ + homer \ + r-base==4.3 \ + r-ggplot2 \ + r-tidyr \ + r-tidyverse \ + rpy2 \ + jupyterlab \ + notebook + + # Install LingerGRN and dependencies + conda activate env + pip install \ + scipy==1.11.3 \ + numpy==1.24.3 \ + pandas==2.0.3 \ + shap==0.42.0 \ + scikit-learn==1.3.0 \ + joblib==1.3.2 \ + matplotlib==3.8.0 \ + seaborn==0.13.0 \ + statsmodels==0.14.1 \ + umap-learn \ + scanpy==1.9.5 \ + anndata==0.9.2 \ + "mudata<0.2.4" \ + LingerGRN==1.106 + + # Remove cache for lighter containers + pip cache purge + conda clean -a -y + apt-get clean + rm -rf /var/lib/apt/lists/* \ No newline at end of file diff --git a/workflow/rules/dbs/linger_prior.smk b/workflow/rules/dbs/linger_prior.smk new file mode 100644 index 00000000..8611c5db --- /dev/null +++ b/workflow/rules/dbs/linger_prior.smk @@ -0,0 +1,20 @@ +import os + +rule linger_prior: + threads: 1 + singularity: 'workflow/envs/linger.sif' + input: 'workflow/envs/linger.sif' + output: + archive='dbs/lingerGRN/data_bulk.tar.gz', + dir=directory('dbs/lingerGRN/data_bulk') + params: + dir_name=lambda wildcards, output: os.path.dirname(output.archive) + shell: + """ + mkdir -p {params.dir_name} + wget -nv --load-cookies /tmp/cookies.txt \ + "https://drive.usercontent.google.com/download?export=download&confirm=$(wget --quiet --save-cookies /tmp/cookies.txt --keep-session-cookies --no-check-certificate 'https://drive.usercontent.google.com/download?id=1jwRgRHPJrKABOk7wImKONTtUupV7yJ9b' -O- | sed -rn 's/.*confirm=([0-9A-Za-z_]+).*/\1/p')&id=1jwRgRHPJrKABOk7wImKONTtUupV7yJ9b" \ + -O '{output.archive}' + rm -f /tmp/cookies.txt + tar -xzf '{output.archive}' -C '{params.dir_name}' + """ \ No newline at end of file diff --git a/workflow/rules/dbs/tss.smk b/workflow/rules/dbs/tss.smk index 9a40dee0..60ecff95 100644 --- a/workflow/rules/dbs/tss.smk +++ b/workflow/rules/dbs/tss.smk @@ -1,4 +1,4 @@ -localrules: gen_tss_celloracle, gen_tss_crema, gen_tss_dictys, gen_tss_figr, gen_tss_granie, gen_tss_pando, gen_tss_scenicplus, gen_tss_scdori, gen_tss_scmtni, gen_tss_promoters +localrules: gen_tss_celloracle, gen_tss_crema, gen_tss_dictys, gen_tss_figr, gen_tss_granie, gen_tss_linger, gen_tss_pando, gen_tss_scenicplus, gen_tss_scdori, gen_tss_scmtni, gen_tss_promoters rule gen_tss_celloracle: @@ -83,6 +83,21 @@ rule gen_tss_inferelator: cp {input} {output} """ +rule gen_tss_linger: + threads: 1 + singularity: 'workflow/envs/gretabench.sif' + input: + gtf=rules.gen_genome_inferelator.output.gtf, + genes='dbs/lingerGRN/data_bulk/bulk_gene_all.txt' + output: 'dbs/hg38/gen/tss/linger.bed.gz' + shell: + """ + python workflow/scripts/dbs/gen/tss/linger.py \ + -i {input.gtf} \ + -g {input.genes} \ + -o {output} + """ + rule gen_tss_pando: threads: 1 singularity: 'workflow/envs/pando.sif' diff --git a/workflow/rules/mth/linger.smk b/workflow/rules/mth/linger.smk new file mode 100644 index 00000000..a36ac067 --- /dev/null +++ b/workflow/rules/mth/linger.smk @@ -0,0 +1,39 @@ +rule mdl_o_linger: + threads: 24 + singularity: 'workflow/envs/linger.sif' + input: + img='workflow/envs/linger.sif', + linger_GRN=rules.linger_prior.output.dir, + mdata=rules.extract_case.output.mdata + output: + dir=directory('dts/{org}/{dat}/cases/{case}/runs/linger/'), + out='dts/{org}/{dat}/cases/{case}/runs/o_linger.o_linger.o_linger.o_linger.mdl.csv' + params: + version=config['methods']['linger']['version'], + mode=config['methods']['linger']['mode'], + organism=lambda w: config['dts'][w.dat]['organism'], + script='workflow/scripts/mth/linger/linger.sh' + resources: + mem_mb=lambda wildcards, attempt: restart_mem(wildcards, attempt) * 2, + runtime=( + 30 if config['methods']['linger']['version'] == 'baseline' + else 140 if config['methods']['linger']['mode'] == 'parallel' + else 360 + ) + shell: + """ + mkdir -p {output.dir} + set -e + timeout $(({resources.runtime}-20))m \ + bash {params.script} \ + --linger_GRN {input.linger_GRN} \ + --out_dir {output.dir} \ + --path_mdata {input.mdata} \ + --version {params.version} \ + --genome {params.organism} \ + --mode {params.mode} \ + --path_out {output.out} + if [ $? -eq 124 ]; then + awk 'BEGIN {{ print "source,target,score,pval" }}' > {output.out} + fi + """ \ No newline at end of file diff --git a/workflow/scripts/dbs/gen/tss/linger.py b/workflow/scripts/dbs/gen/tss/linger.py new file mode 100644 index 00000000..53f35f5f --- /dev/null +++ b/workflow/scripts/dbs/gen/tss/linger.py @@ -0,0 +1,42 @@ +import argparse +import gzip +import re +import pandas as pd + +parser = argparse.ArgumentParser() +parser.add_argument('-g', '--path_genes', required=True) +parser.add_argument('-i', '--path_gtf', required=True) +parser.add_argument('-o', '--path_out', required=True) +args = parser.parse_args() + +gene_set = set(pd.read_csv(args.path_genes, header=None)[0].tolist()) +print(f"LINGER gene universe: {len(gene_set)} genes") + +records = [] +with gzip.open(args.path_gtf, 'rt') as f: + for line in f: + fields = line.rstrip('\n').split('\t') + if fields[2] != 'transcript': + continue + gene_name = re.search(r'gene_name "([^"]+)"', fields[8]).group(1) + records.append([fields[0], int(fields[3]), gene_name]) + +# mirror `get_TSS_ensembl` in LINGER_tr_fast.py +df = pd.DataFrame(records, columns=['chr', 'start', 'symbol']) +df = df.groupby(['chr', 'symbol'], as_index=False)['start'].min() +df = df[df['symbol'] != ''] +df['start'] -= 1 + +std_chroms = [f'chr{i}' for i in range(1, 23)] + ['chrX'] +df = df[df['symbol'].isin(gene_set) & df['chr'].isin(std_chroms)] +print(f"Matched {df['symbol'].nunique()}/{len(gene_set)} genes") + +bed = pd.DataFrame({ + 'chr': df['chr'], + 'start': df['start'], + 'end': df['start'], + 'name': df['symbol'] +}) +bed = bed.sort_values(['chr', 'start']) +bed.to_csv(args.path_out, sep='\t', index=False, header=False, compression='gzip') +print(f"Written {len(bed)} records to {args.path_out}") \ No newline at end of file diff --git a/workflow/scripts/mth/linger/LINGER_tr_fast.py b/workflow/scripts/mth/linger/LINGER_tr_fast.py new file mode 100644 index 00000000..0960a0f4 --- /dev/null +++ b/workflow/scripts/mth/linger/LINGER_tr_fast.py @@ -0,0 +1,547 @@ +import torch +import torch.nn as nn +import torch.optim as optim +from torch.nn import functional as F +from scipy.stats import pearsonr +from scipy.stats import spearmanr +#load data +import numpy as np +import pandas as pd +import random +from torch.optim import Adam +import os +from sklearn.linear_model import ElasticNet +from sklearn.datasets import make_regression +from sklearn.model_selection import KFold +import shap +hidden_size = 64 +hidden_size2 = 16 +output_size = 1 +from joblib import Parallel,delayed +seed_value = 42 + +class Net(nn.Module): + def __init__(self,input_size,activef): + super(Net, self).__init__() + self.fc1 = nn.Linear(input_size, 64) + self.fc2 = nn.Linear(64, 16) + self.fc3 = nn.Linear(16, output_size) + self.activef=activef + def forward(self, x): + #x = torch.sigmoid(self.fc1(x)) + if self.activef=='ReLU': + x = F.relu(self.fc1(x)) + x = F.relu(self.fc2(x)) + if self.activef=='sigmoid': + x = F.sigmoid(self.fc1(x)) + x = F.sigmoid(self.fc2(x)) + if self.activef=='tanh': + x = F.tanh(self.fc1(x)) + x = F.tanh(self.fc2(x)) + x = self.fc3(x) + return x + +#EWC +def EWC(fisher,params,net): + params_n = list(net.parameters()) + EWC=0 + i=0 + p=params_n[0] + cost=(p-params[i])*fisher*(p-params[i]) + EWC=EWC+cost.sum() + return EWC + + +def sc_nn(ii,gene_chr,TFindex,TFindex_bulk,REindex,REindex_bulk,REindex_bulk_match,Target,netall,adj_matrix_all,Exp,TF_match,input_size_all,fisherall,Opn,l1_lambda,fisher_w,activef): + warnings.filterwarnings("ignore") + alpha = 1 + eps=1e-12 + alpha = torch.tensor(alpha,dtype=torch.float32) + gene_idx=gene_chr['id_s'].values[ii]-1 + gene_idx_b=int(gene_chr['id_b'].values[ii])-1 + TFidxtemp=TFindex[gene_idx] + TFidxtemp=TFidxtemp.split('_') + TFidxtemp=[int(TFidxtemp[k])+1 for k in range(len(TFidxtemp))] + TFidxtemp_b=TFindex_bulk[gene_idx_b] + TFidxtemp_b=TFidxtemp_b.split('_') + TFidxtemp_b=[int(TFidxtemp_b[k]) for k in range(len(TFidxtemp_b))] + TFtemp=Exp[np.array(TFidxtemp)-1,:] + REidxtemp=REindex[gene_idx] + REidxtemp_b_m=REindex_bulk_match[gene_idx] + REidxtemp_b=REindex_bulk[gene_idx_b] + REidxtemp=str(REidxtemp).split('_') + REidxtemp_b_m=str(REidxtemp_b_m).split('_') + REidxtemp_b=str(REidxtemp_b).split('_') + if (len(REidxtemp)==1)&(REidxtemp[0]=='nan'): + REidxtemp=[] + REidxtemp_b_m=[] + inputs=TFtemp+1-1 + L=np.zeros([len(TFidxtemp)+len(REidxtemp),len(TFidxtemp)+len(REidxtemp)]) + L=torch.tensor(L, dtype=torch.float32) + else: + REidxtemp=[int(REidxtemp[k])+1 for k in range(len(REidxtemp))] + REidxtemp_b_m=[int(REidxtemp_b_m[k])+1 for k in range(len(REidxtemp_b_m))] + REtemp=Opn[np.array(REidxtemp)-1,:] + inputs=np.vstack((TFtemp, REtemp)) + adj_matrix=np.zeros([len(TFidxtemp)+len(REidxtemp),len(TFidxtemp)+len(REidxtemp)]) + AA=adj_matrix_all[np.array(REidxtemp)-1,:] + AA=AA[:,np.array(TFidxtemp)-1] + adj_matrix[:len(TFidxtemp),-len(REidxtemp):]=AA.T + adj_matrix[-len(REidxtemp):,:len(TFidxtemp)]=AA + A = torch.tensor(adj_matrix, dtype=torch.float32) + D = torch.diag(A.sum(1)) + degree = A.sum(dim=1) + degree += eps + D_sqrt_inv = 1 / degree.sqrt() + D_sqrt_inv = torch.diag(D_sqrt_inv) + L = D_sqrt_inv@(D - A)@D_sqrt_inv + if (len(REidxtemp_b)==1)&(REidxtemp_b[0]=='nan'): + REidxtemp_b=[] + else: + REidxtemp_b=[int(REidxtemp_b[k]) for k in range(len(REidxtemp_b))] + targets = torch.tensor(Target[gene_idx,:]) + inputs = torch.tensor(inputs,dtype=torch.float32) + targets = targets.type(torch.float32) + mean = inputs.mean(dim=1) + std = inputs.std(dim=1) + inputs = (inputs.T - mean) / (std+eps) + inputs=inputs.T + num_nodes=inputs.shape[0] + y=targets.reshape(len(targets),1) + #trainData testData + input_size=int(input_size_all[gene_idx_b]) + loaded_net = Net(input_size,activef) + loaded_net.load_state_dict(netall[gene_idx_b]) + params = list(loaded_net.parameters()) + fisher0=fisherall[gene_idx_b][0].data.clone() + data0=pd.DataFrame(TFidxtemp) + data1=pd.DataFrame(TFidxtemp_b) + data0.columns=['TF'] + data1.columns=['TF'] + A=TF_match.loc[data0['TF'].values-1]['id_b'] + data0=pd.DataFrame(A) + data0.columns=['TF'] + data1['id_b']=data1.index + data0['id_s']=range(0,len(A)) + merge_TF=pd.merge(data0,data1,how='left',on='TF') + if (len(REidxtemp)>0)&(len(REidxtemp_b)>0): + data0=pd.DataFrame(REidxtemp_b_m) + data1=pd.DataFrame(REidxtemp_b) + data0.columns=['RE'] + data1.columns=['RE'] + data0['id_s']=data0.index + data1['id_b']=data1.index + merge_RE=pd.merge(data0,data1,how='left',on='RE') + if merge_RE['id_b'].isna().sum()==0: + good=1 + indexall=merge_TF['id_b'].values.tolist()+(merge_RE['id_b'].values+merge_TF.shape[0]).tolist() + else: + good=0 + else: + indexall=merge_TF['id_b'].values.tolist() + good=1 + if good==1: + fisher=fisher0[:,np.array(indexall,dtype=int)] + params_bulk = params[0][:,np.array(indexall,dtype=int)] + with torch.no_grad(): + params_bulk = params_bulk.detach() + num_nodes=inputs.shape[0] + n_folds = 5 + kf = KFold(n_splits=n_folds,shuffle=True,random_state=0) + fold_size = len(inputs.T) // n_folds + input_size = num_nodes + mse_loss = nn.MSELoss() + y_pred_all=0*(y+1-1) + y_pred_all1=0*(y+1-1) + y_pred_all1=y_pred_all1.numpy().reshape(-1) + X_tr = inputs.T + y_tr = y + torch.manual_seed(seed_value) + net = Net(input_size,activef) + optimizer = Adam(net.parameters(),lr=0.01,weight_decay=l1_lambda) + #optimizer = Adam(net.parameters(),weight_decay=1) + # Perform backpropagation + Loss0=np.zeros([100,1]) + for i in range(100): + # Perform forward pass + y_pred = net(X_tr) + # Calculate loss + l1_norm = sum(torch.linalg.norm(p, 1) for p in net.parameters()) + #loss_EWC=EWC(fisher,params_bulk,net); + l2_bulk = -1* fisher_w* sum(sum(torch.mul(params_bulk,net.fc1.weight))) + lap_reg = alpha * torch.trace(torch.mm(torch.mm(net.fc1.weight, L), net.fc1.weight.t())) + loss = mse_loss(y_pred, y_tr) +l1_norm*l1_lambda+l2_bulk+lap_reg + Loss0[i,0]=loss.detach().numpy() + # Perform backpropagation + optimizer.zero_grad() + loss.backward() + optimizer.step() + np.random.seed(42) + background = X_tr[np.random.choice(X_tr.shape[0], 50, replace=False)] + explainer = shap.DeepExplainer(net,background) + shap_values = explainer.shap_values(X_tr) + warnings.resetwarnings() + return net,shap_values,0.5,0.5,1,Loss0 + else: + warnings.resetwarnings() + return 0,0,0,0,0,0 + +def get_TSS(GRNdir,genome,TSS_dis): + #import pyensembl +# Initialize Ensembl database for the desired genome assembly + #ensembl = pyensembl.EnsemblRelease(release=release, species=species) # For hg19 +# ensembl = pyensembl.EnsemblRelease(release=104, species='mouse') # For mm10 +# Get all genes in the genome + #genes = ensembl.genes() +# Retrieve TSS positions for each gene and store them in a list + #tss_positions = [] + #strand=[] + #chrom=[] + #genesymbol=[] + #for gene in genes: + #tss_positions.append(gene.transcripts[0].start) + #strand.append(gene.strand) + #chrom.append('chr'+gene.contig) + #genesymbol.append(gene.name) + import pandas as pd + Tssdf = pd.read_csv(GRNdir+'TSS_'+genome+'.txt',sep='\t',header=None) + Tssdf.columns=['chr','TSS','symbol','strand'] + Tssdf['1M-']=Tssdf['TSS']-TSS_dis + Tssdf['1M+']=Tssdf['TSS']+TSS_dis + temp=Tssdf['1M-'].values + temp[temp<1]=1 + Tssdf['1M-']=temp + Tssdf=Tssdf[Tssdf['symbol']!=''] + Tssdf[['chr','1M-','1M+','symbol','TSS', 'strand']].to_csv('data/TSS_extend_1M.txt',sep='\t',index=None) + +def load_data(GRNdir,outdir): + gene_all=pd.DataFrame([]) + for i in range(22): + chr='chr'+str(i+1) + gene_file=GRNdir+chr+'_gene.txt' + data0=pd.read_csv(gene_file,sep='\t',header=None) + data0['chr']=chr + data0['id_b']=data0.index+1 + gene_all=pd.concat([gene_all,data0]) + chr='chrX' + gene_file=GRNdir+chr+'_gene.txt' + data0=pd.read_csv(gene_file,sep='\t',header=None) + data0['chr']=chr + data0['id_b']=data0.index+1 + gene_all=pd.concat([gene_all,data0]) + gene_file=outdir+'Symbol.txt' + data0=pd.read_csv(gene_file,sep='\t',header=None) + data0.columns=['Symbol'] + data0['id_s']=data0.index+1 + gene_all.columns=['Symbol','chr','id_b'] + data_merge=pd.merge(data0,gene_all,how='left',on='Symbol') + TFName_b=pd.read_csv(GRNdir+'TFName.txt',header=None,sep='\t') + TFName_s=pd.read_csv(outdir+'TFName.txt',header=None,sep='\t') + TFName_b.columns=['TF'] + TFName_s.columns=['TF'] + TFName_b['id_b']=TFName_b.index+1# index from 1 + TFName_s['id_s']=TFName_s.index+1# index from 1 + TF_match=pd.merge(TFName_s,TFName_b,how='left',on='TF') + Opn_file=outdir+'Openness.txt' + idx_file=outdir+'index.txt' + geneexp_file=outdir+'Exp.txt' + Target=pd.read_csv(geneexp_file,header=None,sep='\t') + Target=Target.values + #def sc_NN(gene_file,Opn_file,idx_file,geneexp_file,out_PCC,out_net): + #alpha = torch.tensor(alpha,dtype=torch.float32) + bind_file=outdir+'TF_binding.txt' + adj_matrix_all=pd.read_csv(bind_file,header=None,sep='\t') + adj_matrix_all=adj_matrix_all.values + TFExp_file=outdir+'TFexp.txt' + Opn=pd.read_csv(Opn_file,header=None,sep='\t') + Opn=Opn.values + idx=pd.read_csv(idx_file,header=None,sep='\t') + Exp=pd.read_csv(TFExp_file,header=None,sep='\t') + Exp=Exp.values + return Exp,idx,Opn,adj_matrix_all,Target,data_merge,TF_match +def sc_nn_NN(ii,RE_TGlink_temp,Target,Exp,Opn,l1_lambda,activef): + warnings.filterwarnings("ignore") + alpha = 1 + eps=1e-12 + alpha = torch.tensor(alpha,dtype=torch.float32) + if RE_TGlink_temp[0] in Exp.index: + TFtemp = Exp.drop([RE_TGlink_temp[0]]).values + else: + TFtemp=Exp.values + REtemp=Opn.loc[RE_TGlink_temp[1]].values + inputs=np.vstack((TFtemp, REtemp)) + targets = torch.tensor(Target.loc[RE_TGlink_temp[0],:]) + inputs = torch.tensor(inputs,dtype=torch.float32) + targets = targets.type(torch.float32) + mean = inputs.mean(dim=1) + std = inputs.std(dim=1) + inputs = (inputs.T - mean) / (std+eps) + inputs=inputs.T + num_nodes=inputs.shape[0] + y=targets.reshape(len(targets),1) + #trainData testData + input_size=int(num_nodes) + mse_loss = nn.MSELoss() + y_pred_all=0*(y+1-1) + y_pred_all1=0*(y+1-1) + y_pred_all1=y_pred_all1.numpy().reshape(-1) + X_tr = inputs.T + y_tr = y + torch.manual_seed(seed_value) + net = Net(input_size,activef) + optimizer = Adam(net.parameters(),lr=0.01,weight_decay=l1_lambda) + #optimizer = Adam(net.parameters(),weight_decay=1) + # Perform backpropagation + Loss0=np.zeros([100,1]) + for i in range(100): + # Perform forward pass + y_pred = net(X_tr) + # Calculate loss + l1_norm = sum(torch.linalg.norm(p, 1) for p in net.parameters()) + #loss_EWC=EWC(fisher,params_bulk,net); + #l2_bulk = -1* fisher_w* sum(sum(torch.mul(params_bulk,net.fc1.weight))) + #lap_reg = alpha * torch.trace(torch.mm(torch.mm(net.fc1.weight, L), net.fc1.weight.t())) + loss = mse_loss(y_pred, y_tr) +l1_norm*l1_lambda#+l2_bulk+lap_reg + Loss0[i,0]=loss.detach().numpy() + # Perform backpropagation + optimizer.zero_grad() + loss.backward() + optimizer.step() + np.random.seed(42) + background = X_tr[np.random.choice(X_tr.shape[0], 50, replace=False)] + explainer = shap.DeepExplainer(net,background) + shap_values = explainer.shap_values(X_tr) + warnings.resetwarnings() + return net,shap_values,Loss0 + + +def load_data_scNN(GRNdir,species): + import pandas as pd + if species=='New': + Match2=pd.read_csv(GRNdir+'MotifMatch.txt',header=0,sep='\t') + else: + Match2=pd.read_csv(GRNdir+'Match_TF_motif_'+species+'.txt',header=None,sep='\t') + Match2.columns = ['Motif','TF'] + TFName = pd.DataFrame(Match2['TF'].unique()) + Target=pd.read_csv('data/TG_pseudobulk.tsv',sep=',',header=0,index_col=0) + TFlist=list(set(Target.index)&set(TFName[0].values)) + Exp=Target.loc[TFlist] + Opn=pd.read_csv('data/RE_pseudobulk.tsv',sep=',',header=0,index_col=0) + RE_TGlink=pd.read_csv('data/RE_gene_distance.txt',sep='\t',header=0) + RE_TGlink = RE_TGlink.groupby('gene').apply(lambda x: x['RE'].values.tolist()).reset_index() + geneoverlap=list(set(Target.index)&set(RE_TGlink['gene'])) + RE_TGlink.index=RE_TGlink['gene'] + RE_TGlink=RE_TGlink.loc[geneoverlap] + RE_TGlink=RE_TGlink.reset_index(drop=True) + return Exp,Opn,Target,RE_TGlink + +def RE_TG_dis(outdir): + import pandas as pd + import pybedtools + import numpy as np + print('Overlap the regions with gene loc ...') + import os# Create the directory + current_directory = os.getcwd() + os.makedirs(outdir, exist_ok=True) + import pandas as pd + peakList=pd.read_csv(current_directory+'/data/Peaks.txt',index_col=None,header=None) + peakList1=[temp.split(':')[0] for temp in peakList[0].values.tolist()] + peakList2=[temp.split(':')[1].split('-')[0] for temp in peakList[0].values.tolist()] + peakList3=[temp.split(':')[1].split('-')[1] for temp in peakList[0].values.tolist()] + peakList['chr']=peakList1 + peakList['start']=peakList2 + peakList['end']=peakList3 + peakList[['chr','start','end']].to_csv(current_directory+'/data/Peaks.bed',sep='\t',header=None,index=None) + TSS_1M=pd.read_csv(current_directory+'/data/TSS_extend_1M.txt',sep='\t',header=0) + TSS_1M.to_csv(current_directory+'/data/TSS_extend_1M.bed',sep='\t',header=None,index=None) + a = pybedtools.example_bedtool(current_directory+'/data/Peaks.bed') + b = pybedtools.example_bedtool(current_directory+'/data/TSS_extend_1M.bed') + a_with_b = a.intersect(b, wa=True,wb=True) + a_with_b.saveas(outdir+'temp.bed') + a_with_b=pd.read_csv(outdir+'temp.bed',sep='\t',header=None) + a_with_b['RE']=a_with_b[0].astype(str) + ':' + a_with_b[1].astype(str) + '-' + a_with_b[2].astype(str) + temp=a_with_b[['RE',6]] + temp.columns=[['RE','gene']] + temp['distance']=np.abs(a_with_b[7]-a_with_b[1]) + temp.to_csv(current_directory+'/data/RE_gene_distance.txt',sep='\t',index=None) + +def get_system_resources(): + import psutil + cpus = psutil.cpu_count(logical=True) + ram = psutil.virtual_memory().total + + # Override with SLURM if available + cpus = int(os.environ.get('SLURM_CPUS_PER_TASK', cpus)) + ram_mb = int(os.environ.get('SLURM_MEM_PER_NODE', ram // (1024*1024))) + ram_gb = ram_mb / 1024 + + return cpus, ram_gb + +def _process_chr(chr, GRNdir, outdir, data_merge, idx, Target, adj_matrix_all, Exp, TF_match, Opn, l1_lambda, fisher_w, activef): + import warnings + import torch + import pandas as pd + import numpy as np + from tqdm import tqdm + + torch.set_num_threads(1) + + netall_s = {} + shapall_s = {} + result = np.zeros([data_merge.shape[0], 2]) + Lossall = np.zeros([data_merge.shape[0], 100]) + + #print(chr, flush=True) + + idx_file1 = GRNdir + chr + '_index.txt' + idx_file_all = GRNdir + chr + '_index_all.txt' + idx_bulk = pd.read_csv(idx_file1, header=None, sep='\t') + idxRE_all = pd.read_csv(idx_file_all, header=None, sep='\t') + + gene_chr = data_merge[data_merge['chr'] == chr] + N = len(gene_chr) + + TFindex = idx.values[:, 2] + REindex = idx.values[:, 1] + REindex_bulk_match = idx.values[:, 3] + REindex_bulk = idxRE_all.values[:, 0] + TFindex_bulk = idx_bulk.values[:, 2] + input_size_all = idx_bulk.values[:, 3] + + fisherall = torch.load(GRNdir + 'fisher_' + chr + '.pt', weights_only=False) + netall = torch.load(GRNdir + 'all_models_' + chr + '.pt', weights_only=False) + + for ii in tqdm(range(N), desc=chr): + warnings.filterwarnings("ignore") + res = sc_nn( + ii, gene_chr, TFindex, TFindex_bulk, REindex, REindex_bulk, + REindex_bulk_match, Target, netall, adj_matrix_all, Exp, + TF_match, input_size_all, fisherall, Opn, l1_lambda, fisher_w, activef + ) + warnings.resetwarnings() + index_all = gene_chr.index[ii] + if res[4] == 1: + result[index_all, 0] = res[2] + result[index_all, 1] = res[3] + netall_s[index_all] = res[0] + shapall_s[index_all] = res[1] + Lossall[index_all, :] = res[5].T + else: + result[index_all, 0] = -100 + + result = pd.DataFrame(result) + result.index = data_merge['Symbol'].values + genetemp = data_merge[data_merge['chr'] == chr]['Symbol'].values + result = result.loc[genetemp] + result.to_csv(outdir + 'result_' + chr + '.txt', sep='\t') + torch.save(netall_s, outdir + 'net_' + chr + '.pt') + torch.save(shapall_s, outdir + 'shap_' + chr + '.pt') + Lossall = pd.DataFrame(Lossall) + Lossall.index = data_merge['Symbol'].values + Lossall = Lossall.loc[genetemp] + Lossall.to_csv(outdir + 'Loss_' + chr + '.txt', sep='\t') + +from tqdm import tqdm +import warnings +import time +import pandas as pd +import numpy as np +def training(GRNdir,method,outdir,activef,species): + import psutil + if method=='LINGER': + hidden_size = 64 + hidden_size2 = 16 + output_size = 1 + l1_lambda = 0.01 + alpha_l = 0.01#elastic net parameter + lambda0 = 0.00 #bulk + fisher_w=0.1 + + n_cpus, ram_av = get_system_resources() + base_ram = psutil.Process(os.getpid()).memory_info().rss / 1e9 + #print(f"Base ram : {base_ram_gb:.2f}") + + # available RAM (GB) beyond the base RAM usage and safety RAM (10GB) + ram_free = ram_av - base_ram - 10 + + # per worker RAM (load fisherfisher_{chr}.pt, net_{chr}.pt and create shap_{chr}.pt) + ram_worker = 5 + + n_jobs = max(1, min(int(ram_free / ram_worker), n_cpus, 23)) + print(f"With {n_cpus} CPUs, {ram_av:.2f} GB RAM (avail. : {ram_free:.2f} GB), start {n_jobs} workers") + + Exp,idx,Opn,adj_matrix_all,Target,data_merge,TF_match=load_data(GRNdir,outdir) + data_merge.to_csv(outdir+'data_merge.txt',sep='\t') + chrall=['chr'+str(i+1) for i in range(22)] + chrall.append('chrX') + Parallel(n_jobs=n_jobs, backend='loky', verbose=10)( + delayed(_process_chr)( + chr, GRNdir, outdir, data_merge, idx, Target, + adj_matrix_all, Exp, TF_match, Opn, + l1_lambda, fisher_w, activef + ) + for chr in chrall + ) + if method=='scNN': + hidden_size = 64 + hidden_size2 = 16 + output_size = 1 + l1_lambda = 0.01 + alpha_l = 0.01#elastic net parameter + lambda0 = 0.00 #bulk + fisher_w=0.1 + n_jobs=16 + Exp,Opn,Target,RE_TGlink=load_data_scNN(GRNdir,species) + import warnings + import time + from tqdm import tqdm + netall_s={} + shapall_s={} + #result=np.zeros([data_merge.shape[0],2]) + chrall=[RE_TGlink[0][i][0].split(':')[0] for i in range(RE_TGlink.shape[0])] + RE_TGlink['chr']=chrall + chrlist=RE_TGlink['chr'].unique() + for jj in tqdm(range(len(chrlist))): + chrtemp=chrlist[jj] + RE_TGlink1=RE_TGlink[RE_TGlink['chr']==chrtemp] + Lossall=np.zeros([RE_TGlink1.shape[0],100]) + for ii in range(RE_TGlink1.shape[0]): + warnings.filterwarnings("ignore") + #res = Parallel(n_jobs=n_jobs)(delayed(sc_nn_NN)(ii,RE_TGlink_temp,Target,netall,Exp,Opn,l1_lambda,activef) for ii in tqdm(range(RE_TGlink.shape[0])) + RE_TGlink_temp=RE_TGlink1.values[ii,:] + res=sc_nn_NN(ii,RE_TGlink_temp,Target,Exp,Opn,l1_lambda,activef) + warnings.resetwarnings() + netall_s[ii]=res[0] + shapall_s[ii]=res[1] + Lossall[ii,:]=res[2].T + torch.save(netall_s,outdir+chrtemp+'_net.pt') + torch.save(shapall_s,outdir+chrtemp+'_shap.pt') + Lossall=pd.DataFrame(Lossall) + Lossall.index=RE_TGlink1['gene'].values + Lossall.to_csv(outdir+chrtemp+'_Loss.txt',sep='\t') + RE_TGlink.to_csv(outdir+'RE_TGlink.txt',sep='\t',index=None) + + +def get_TSS_ensembl(genome_short,gtf_file,GRNdir): + import pyensembl + import subprocess + from pyensembl import Genome + ensembl = Genome( + reference_name=genome_short, + annotation_name="My_annotation", + gtf_path_or_url=gtf_file) + ensembl.index() + genes = ensembl.genes() +# Retrieve TSS positions for each gene and store them in a list + tss_positions = [] + strand=[] + chrom=[] + genesymbol=[] + for gene in genes: + tss_positions.append(gene.transcripts[0].start) + strand.append(gene.strand) + chrom.append('chr'+gene.contig) + genesymbol.append(gene.name) + import pandas as pd + Tssdf = pd.DataFrame({'chr': chrom, 'TSS': tss_positions, 'symbol': genesymbol,'strand': strand}) + Tssdf.to_csv(GRNdir+'TSS_'+genome_short+'.txt',sep='\t',index=None,header=0) \ No newline at end of file diff --git a/workflow/scripts/mth/linger/LL_net_fast.py b/workflow/scripts/mth/linger/LL_net_fast.py new file mode 100644 index 00000000..de58f9a1 --- /dev/null +++ b/workflow/scripts/mth/linger/LL_net_fast.py @@ -0,0 +1,1160 @@ +import numpy as np +import pandas as pd +from scipy.sparse import coo_matrix +from scipy.sparse import csc_matrix +from tqdm import tqdm +import torch +import csv +import torch.nn as nn +import torch.optim as optim +from torch.nn import functional as F +from scipy.stats import pearsonr +from scipy.stats import spearmanr +#load data +import random +from torch.optim import Adam +import os +from sklearn.linear_model import ElasticNet +from sklearn.datasets import make_regression +from sklearn.model_selection import KFold +hidden_size = 64 +hidden_size2 = 16 +output_size = 1 +seed_value = 42 +torch.manual_seed(seed_value) +class Net(nn.Module): + def __init__(self,input_size,activef): + super(Net, self).__init__() + self.fc1 = nn.Linear(input_size, 64) + self.fc2 = nn.Linear(64, 16) + self.fc3 = nn.Linear(16, output_size) + self.activef=activef + def forward(self, x): + #x = torch.sigmoid(self.fc1(x)) + if self.activef=='ReLU': + x = F.relu(self.fc1(x)) + x = F.relu(self.fc2(x)) + if self.activef=='sigmoid': + x = F.sigmoid(self.fc1(x)) + x = F.sigmoid(self.fc2(x)) + if self.activef=='tanh': + x = F.tanh(self.fc1(x)) + x = F.tanh(self.fc2(x)) + x = self.fc3(x) + return x + +def cosine_similarity_0(X): + A=X.T/((X**2).sum(axis=1)**(1/2)+((X**2).sum(axis=1)**(1/2)).mean()/1000000) + return np.dot(A.T,A) + + +def list2mat(df,i_n,j_n,x_n): + TFs = df[j_n].unique() + REs = df[i_n].unique() +#Initialize matrix as numpy array +#Map row and col indices for lookup + row_map = {r:i for i,r in enumerate(REs)} + col_map = {c:i for i,c in enumerate(TFs)} + row_indices = np.array([row_map[row] for row in df[i_n]]) + col_indices = np.array([col_map[col] for col in df[j_n]]) + from scipy.sparse import coo_matrix + matrix = coo_matrix((df[x_n], (row_indices, col_indices)), shape=(len(REs), len(TFs))) + mat=coo_matrix.toarray(matrix) + return mat,REs,TFs + + + +def list2mat_s(df,REs,TFs,i_n,j_n,x_n): +#Initialize matrix as numpy array +#Map row and col indices for lookup + row_map = {r:i for i,r in enumerate(REs)} + col_map = {c:i for i,c in enumerate(TFs)} + row_indices = np.array([row_map[row] for row in df[i_n]]) + col_indices = np.array([col_map[col] for col in df[j_n]]) + from scipy.sparse import coo_matrix + import scipy.sparse as sp + matrix = sp.csr_matrix((df[x_n], (row_indices, col_indices)), shape=(len(REs), len(TFs))) + return matrix,REs,TFs + +def merge_columns_in_bed_file(file_path,startcol): + merged_values = [] + with open(file_path, 'r') as file: + for line in file: + columns = line.strip().split('\t') + col1 = columns[-1+startcol] + col2 = columns[startcol] + col3 = columns[1+startcol] + merged_value = f"{col1}:{col2}-{col3}" + merged_values.append(merged_value) + return merged_values +def merge_columns_in_bed_file2(file_path,startcol): + merged_values = [] + with open(file_path, 'r') as file: + for line in file: + columns = line.strip().split('\t') + col1 = columns[-1+startcol] + col2 = columns[startcol] + col3 = columns[1+startcol] + merged_value = f"{col1}_{col2}_{col3}" + merged_values.append(merged_value) + return merged_values +def format_RE_tran12(region): + chr, range_ = region.split(":") + start, end = range_.split("-") + return "_".join([chr, start, end]) +def get_TF_RE(data_merge_temp,j,net_all,TFindex,TFName,REindex,REName): + index_all=data_merge_temp[j] + result={'TF':[],'RE':[],'score':[]} + result = pd.DataFrame(result) + #for ii in range(1): + temps=list(net_all[index_all].parameters())[0] + TFidxtemp=TFindex[index_all] + TFidxtemp=TFidxtemp.split('_') + TFidxtemp=[int(TFidxtemp[i]) for i in range(len(TFidxtemp))] + TFName_temp=TFName[np.array(TFidxtemp)] + REidxtemp=REindex[index_all] + if REidxtemp=='': + REidxtemp=[] + else: + REidxtemp=REidxtemp.split('_') + REidxtemp=[int(REidxtemp[i]) for i in range(len(REidxtemp))] #146 RE idx + if len(REidxtemp)>0: + corr_matrix = cosine_similarity_0(temps.detach().numpy().T) + REName_temp=REName[np.array(REidxtemp)] + corr_matrix=corr_matrix[:len(TFidxtemp),len(TFidxtemp):] + for k in range(len(REidxtemp)): + datatemp=pd.DataFrame({'score':corr_matrix[:,k].tolist()}) + datatemp['TF']=TFName_temp.tolist() + datatemp['RE']=REName_temp[k] + result=pd.concat([result,datatemp]) + return result +def load_TFbinding(GRNdir,O_overlap,O_overlap_u,O_overlap_hg19_u,chrN): + TFbinding=pd.read_csv(GRNdir+'TF_binding_'+chrN+'.txt',sep='\t',index_col=0) + TFbinding1=np.zeros([len(O_overlap_u),TFbinding.shape[1]]) + TFbinding1=np.zeros([len(O_overlap_u),TFbinding.shape[1]]) + O_overlap1=list(set(O_overlap_hg19_u)&set(TFbinding.index)) + List=pd.DataFrame(range(len(TFbinding.index)), index=TFbinding.index) + index0=List.loc[O_overlap1][0].values + #O_overlap_df=pd.DataFrame(range(len(O_overlap)), index=O_overlap) + O_overlap_hg19_u_df=pd.DataFrame(range(len(O_overlap_hg19_u)), index=O_overlap_hg19_u) + #hg19_38=pd.DataFrame(O_overlap_u, index=O_overlap_hg19_u) + index1=O_overlap_hg19_u_df.loc[O_overlap1][0].values + #index1=O_overlap_df.loc[index1][0].values + TFbinding1[index1,:]=TFbinding.iloc[index0,:].values + #TFbinding=pd.DataFrame(TFbinding1,index=O_overlap_u,columns=TFbinding.columns) + O_overlap_u_df=pd.DataFrame(range(len(O_overlap_u)), index=O_overlap_u) + hg19_38=pd.DataFrame(O_overlap_u, index=O_overlap_hg19_u) + TFbinding2=np.zeros([len(O_overlap),TFbinding.shape[1]]) + index=O_overlap_u_df.loc[O_overlap][0].values + TFbinding2=TFbinding1[index,:] + TFbinding=pd.DataFrame(TFbinding2,index=O_overlap,columns=TFbinding.columns) + return TFbinding + +def load_region(GRNdir,genome,chrN,outdir): + O_overlap=merge_columns_in_bed_file(outdir+'Region_overlap_'+chrN+'.bed',1) + N_overlap=merge_columns_in_bed_file(outdir+'Region_overlap_'+chrN+'.bed',4) + O_overlap_u=list(set(O_overlap)) + N_overlap_u=list(set(N_overlap)) + #O_all=merge_columns_in_bed_file(GRNdir+'Peaks_'+chrN+'.bed',1) + hg19_region=merge_columns_in_bed_file(GRNdir+'hg19_Peaks_'+chrN+'.bed',1) + hg19_region=pd.DataFrame(range(len(hg19_region)),index=hg19_region) + hg38_region=merge_columns_in_bed_file(GRNdir+'hg38_Peaks_'+chrN+'.bed',1) + hg38_region=pd.DataFrame(range(len(hg38_region)),index=hg38_region) + if genome=='hg19': + idx=hg19_region.loc[O_overlap_u][0].values + O_overlap_u=hg38_region.index[idx].tolist() + O_overlap_hg19_u=hg19_region.index[idx].tolist() + if genome=='hg38': + idx=hg38_region.loc[O_overlap_u][0].values + O_overlap_hg19_u=hg19_region.index[idx].tolist() + return O_overlap, N_overlap,O_overlap_u,N_overlap_u,O_overlap_hg19_u + +def load_TF_RE(GRNdir,chrN,O_overlap,O_overlap_u,O_overlap_hg19_u): + #print('load prior TF-RE for '+chrN+'...') + mat=pd.read_csv(GRNdir+'Primary_TF_RE_'+chrN+'.txt',sep='\t',index_col=0) + mat1=np.zeros([len(O_overlap_u),mat.shape[1]]) + O_overlap1=list(set(O_overlap_u)&set(mat.index)) + List=pd.DataFrame(range(len(mat.index)), index=mat.index) + index0 = List.loc[O_overlap1][0].values + O_overlap_u_df=pd.DataFrame(range(len(O_overlap_u)), index=O_overlap_u) + index1 = O_overlap_u_df.loc[O_overlap1][0].values + mat1[index1,:] = mat.iloc[index0,:].values + #mat = pd.DataFrame(mat1,index=O_overlap_u,columns=mat.columns) + O_overlap_u_df=pd.DataFrame(range(len(O_overlap_u)), index=O_overlap_u) + hg19_38=pd.DataFrame(O_overlap_u, index=O_overlap_hg19_u) + mat2=np.zeros([len(O_overlap),mat.shape[1]]) + index=O_overlap_u_df.loc[O_overlap][0].values + mat2=mat1[index,:] + mat=pd.DataFrame(mat2,index=O_overlap,columns=mat.columns) + return mat +def TF_RE_LINGER_chr(chr,outdir): + REName = 'data/Peaks.txt' +# Open the file in read mode + with open(REName, "r") as file: + # Create a CSV reader + reader = csv.reader(file, delimiter='\t') + # Read the first column and store it in a list + first_column = [row[0] for row in reader] + REName=np.array(first_column) + idx_file=outdir+'index.txt' + from scipy.stats import zscore + data0=pd.read_csv(outdir+'result_'+chr+'.txt',sep='\t') + data0.columns=['gene','x','y'] + idx_file=outdir+'index.txt' + idx=pd.read_csv(idx_file,sep='\t',header=None) + idx.columns=['gene','REid','TF_id','REid_b'] + idx.fillna('', inplace=True) + TFName=outdir+'TFName.txt' + TFName=pd.read_csv(TFName,sep='\t',header=None) + TFName.columns=['Name'] + TFName=TFName['Name'].values + TFindex=idx['TF_id'].values + REindex=idx['REid'].values + geneName=idx['gene'].values + net_all=torch.load(outdir+"net_"+chr+".pt", weights_only=False) + data_merge=pd.read_csv(outdir+'data_merge.txt',sep='\t',header=0,index_col=0) + data_merge_temp=data_merge[data_merge['chr']==chr].index + batchsize=50 + AAA=np.abs(data0[['x']].values) + N=data_merge_temp.shape[0] + times=int(np.floor(N/batchsize)) + resultlist=[0 for i in range(times+1)] + for ii in tqdm(range(times), desc=chr): + result_all=pd.DataFrame([]) + for j in range(ii*batchsize,(ii+1)*batchsize): + if (AAA[j]>0)&(AAA[j]<10): + result=get_TF_RE(data_merge_temp,j,net_all,TFindex,TFName,REindex,REName) + result_all=pd.concat([result_all,result],axis=0) + result_all=result_all.groupby(['TF', 'RE'])['score'].max().reset_index() + resultlist[ii]=result_all + result_all=pd.DataFrame([]) + ii=ii+1 + for j in range(ii*batchsize,N): + if (AAA[j]>0)&(AAA[j]<10): + result=get_TF_RE(data_merge_temp,j,net_all,TFindex,TFName,REindex,REName) + result_all=pd.concat([result_all,result],axis=0) + if result_all.shape[0]>0: + result_all=result_all.groupby(['TF', 'RE'])['score'].max().reset_index() + resultlist[ii]=result_all + result_all1=pd.concat(resultlist,axis=0) + A=result_all1.groupby(['TF', 'RE'])['score'].max().reset_index() + mat,REs,TFs=list2mat(A,'RE','TF','score') + mat=pd.DataFrame(mat,index=REs,columns=TFs) + return mat + +def _process_TF_RE_LINGER_chr(chrN, outdir, TG_index): + import pandas as pd + mat = TF_RE_LINGER_chr(chrN, outdir) + TFoverlap = list(set(mat.columns) & set(TG_index)) + mat = mat[TFoverlap] + mat.to_csv(outdir + chrN + '_cell_population_TF_RE_binding.txt', sep='\t') + return mat + +def TF_RE_binding_chr(adata_RNA,adata_ATAC,GRNdir,chrN,genome,outdir): + ## the regions + O_overlap, N_overlap,O_overlap_u,N_overlap_u,O_overlap_hg19_u=load_region(GRNdir,genome,chrN,outdir) + import numpy as np + import pandas as pd +## read the count file. + #RE=pd.DataFrame(adata_ATAC.raw.X.toarray().T,index=adata_ATAC.raw.var['gene_ids'].values,columns=adata_ATAC.obs['barcode'].values) + TG=pd.DataFrame(adata_RNA.X.toarray().T,index=adata_RNA.var['gene_ids'].values,columns=adata_RNA.obs['barcode'].values) +## cell annotation +## extact the overlapped peaks. + #RE=RE.loc[N_overlap] + #TFbinding=load_TFbinding(GRNdir,O_overlap,O_overlap_u,O_overlap_hg19_u,chrN) + mat=load_TF_RE(GRNdir,chrN,O_overlap,O_overlap_u,O_overlap_hg19_u) + TFs = mat.columns + TFoverlap = list(set(TFs) & set(TG.index)) + mat = mat[TFoverlap] + #TFbinding = TFbinding[TFoverlap] + #TF = TG.loc[TFoverlap] + #mat_m=np.mean(mat.values[mat>0]) + #mat = mat / mat_m + mat.values[mat.values<0]=0 + #TFbinding = TFbinding / TFbinding.mean(axis=1).mean() + #TF_cluster = TF.values.mean(axis=1) + #TF_cluster = TF_cluster[None,:] + #RE_cluster = RE.values.mean(axis=1) + #RE_cluster = RE_cluster[:,None] + #S = np.log(RE_cluster+0.1) + np.log(mat+TFbinding+0.1) + np.log(TF_cluster+0.1) + S = mat#+TFbinding + #S = np.exp(S) + S.index=N_overlap + mean_S = S.groupby(S.index).max() + return mean_S +import ast +def TF_RE_scNN(TFName,geneName,net_all,RE_TGlink,REName): + batchsize=50 + REName=pd.DataFrame(range(len(REName)),index=REName) + N=RE_TGlink.shape[0] + times=int(np.floor(N/batchsize)) + resultlist=[0 for i in range(times+1)] + for ii in range(times): + result_all=pd.DataFrame([]) + for j in range(ii*batchsize,(ii+1)*batchsize): + RE_TGlink_temp=RE_TGlink.values[j,:] + temps=list(net_all[j].parameters())[0] + actual_list = ast.literal_eval(RE_TGlink_temp[1]) + REidxtemp=REName.loc[actual_list].index + TFidxtemp=np.array(range(len(TFName))) + TFidxtemp=TFidxtemp[TFName!=RE_TGlink_temp[0]] + if len(REidxtemp)>0: + corr_matrix = cosine_similarity_0(temps.detach().numpy().T) + corr_matrix=corr_matrix[:len(TFidxtemp),len(TFidxtemp):] + result={'TF':[],'RE':[],'score':[]} + result = pd.DataFrame(result) + for k in range(len(REidxtemp)): + datatemp=pd.DataFrame({'score':corr_matrix[:,k].tolist()}) + datatemp['TF']=TFName[TFidxtemp].tolist() + datatemp['RE']=REidxtemp[k] + result=pd.concat([result,datatemp]) + result_all=pd.concat([result_all,result],axis=0) + result_all=result_all.groupby(['TF', 'RE'])['score'].max().reset_index() + #print(result_all) + resultlist[ii]=result_all + result_all=pd.DataFrame([]) + ii=times + if N>ii*batchsize: + for j in range(ii*batchsize,N): + RE_TGlink_temp=RE_TGlink.values[j,:] + temps=list(net_all[j].parameters())[0] + actual_list = ast.literal_eval(RE_TGlink_temp[1]) + REidxtemp=REName.loc[actual_list].index + TFidxtemp=np.array(range(len(TFName))) + TFidxtemp=TFidxtemp[TFName!=RE_TGlink_temp[0]] + if len(REidxtemp)>0: + corr_matrix = cosine_similarity_0(temps.detach().numpy().T) + corr_matrix=corr_matrix[:len(TFidxtemp),len(TFidxtemp):] + result={'TF':[],'RE':[],'score':[]} + result = pd.DataFrame(result) + for k in range(len(REidxtemp)): + datatemp=pd.DataFrame({'score':corr_matrix[:,k].tolist()}) + datatemp['TF']=TFName[TFidxtemp].tolist() + datatemp['RE']=REidxtemp[k] + result=pd.concat([result,datatemp]) + result_all=pd.concat([result_all,result],axis=0) + result_all=result_all.groupby(['TF', 'RE'])['score'].max().reset_index() + #print(result_all) + resultlist[ii]=result_all + result_all=pd.concat(resultlist,axis=0) + result_all=result_all.groupby(['TF', 'RE'])['score'].max().reset_index() + return result_all + +def load_data_scNN(GRNdir,genome): + import pandas as pd + genome_map=pd.read_csv(GRNdir+'genome_map_homer.txt',sep='\t',header=0) + genome_map.index=genome_map['genome_short'].values + if genome in genome_map.index: + Match2=pd.read_csv(GRNdir+'Match_TF_motif_'+genome_map.loc[genome]['species_ensembl']+'.txt',sep='\t',header=0) + else: + Match2=pd.read_csv(GRNdir+'MotifMatch.txt',sep='\t',header=0) + TFName = pd.DataFrame(Match2['TF'].unique()) + Target=pd.read_csv('data/TG_pseudobulk.tsv',sep=',',header=0,index_col=0) + TFlist=list(set(Target.index)&set(TFName[0].values)) + Exp=Target.loc[TFlist] + Opn=pd.read_csv('data/RE_pseudobulk.tsv',sep=',',header=0,index_col=0) + RE_TGlink=pd.read_csv('data/RE_gene_distance.txt',sep='\t',header=0) + RE_TGlink = RE_TGlink.groupby('gene').apply(lambda x: x['RE'].values.tolist()).reset_index() + geneoverlap=list(set(Target.index)&set(RE_TGlink['gene'])) + RE_TGlink.index=RE_TGlink['gene'] + RE_TGlink=RE_TGlink.loc[geneoverlap] + RE_TGlink=RE_TGlink.reset_index(drop=True) + return Exp,Opn,Target,RE_TGlink + +def TF_RE_binding(GRNdir,adata_RNA,adata_ATAC,genome,method,outdir): + from tqdm import tqdm + import numpy as np + import pandas as pd + print('Generating cellular population TF binding strength ...') + chrom = ['chr'+str(i+1) for i in range(22)] + chrom.append('chrX') + if method=='baseline': + result=pd.DataFrame() + for i in tqdm(range(23)): + chrN=chrom[i] + out=TF_RE_binding_chr(adata_RNA,adata_ATAC,GRNdir,chrN,genome,outdir) + out.to_csv(outdir+chrN+'_cell_population_TF_RE_binding.txt',sep='\t') + #result=pd.concat([result,out],axis=1).fillna(0) + result = pd.concat([result, out], join='outer', axis=0) + + if method=='LINGER': + import psutil, os + from joblib import Parallel, delayed + TG_index = adata_RNA.var['gene_ids'].values.tolist() + n_cpus = int(os.environ.get('SLURM_CPUS_PER_TASK', os.cpu_count())) + n_jobs = min(n_cpus, 23) # assuming n_cpus <= n_gb_ram (1GB/worker) + results = Parallel(n_jobs=n_jobs, backend='loky', verbose=5)( + delayed(_process_TF_RE_LINGER_chr)(chrN, outdir, TG_index) + for chrN in chrom + ) + result = pd.concat(results, join='outer', axis=0) + + if method=='scNN': + Exp,Opn,Target,RE_TGlink=load_data_scNN(GRNdir,genome) + RE_TGlink=pd.read_csv(outdir+'RE_TGlink.txt',sep='\t',header=0) + RE_TGlink.columns=[0,1,'chr'] + #chrall=[RE_TGlink[0][i][0].split(':')[0] for i in range(RE_TGlink.shape[0])] + chrlist=RE_TGlink['chr'].unique() + REName=Opn.index + geneName=Target.index + TFName=Exp.index + result_all=pd.DataFrame([]) + for jj in tqdm(range(0,len(chrlist))): + chrtemp=chrlist[jj] + RE_TGlink1=RE_TGlink[RE_TGlink['chr']==chrtemp] + net_all=torch.load(outdir+chrtemp+'_net.pt') + result=TF_RE_scNN(TFName,geneName,net_all,RE_TGlink1,REName) + result.to_csv(outdir+chrtemp+'_cell_population_TF_RE_binding.txt',sep='\t') + result_all=pd.concat([result_all,result],axis=0) + result=result_all.copy() + result.to_csv(outdir+'cell_population_TF_RE_binding.txt',sep='\t') + +def load_TFbinding_scNN(GRNdir,outdir,genome): + import pandas as pd + import numpy as np + genome_map=pd.read_csv(GRNdir+'genome_map_homer.txt',sep='\t',header=0) + genome_map.index=genome_map['genome_short'].values + A=pd.read_csv(outdir+'MotifTarget.bed',sep='\t',header=0,index_col=None) + #Motif_binding,REs1,motifs=list2mat(A,'PositionID','Motif Name','MotifScore') + A['MotifScore']=np.log(1+A['MotifScore']); + if genome in genome_map.index: + Match2=pd.read_csv(GRNdir+'Match_TF_motif_'+genome_map.loc[genome]['species_ensembl']+'.txt',sep='\t',header=0) + else: + Match2=pd.read_csv(GRNdir+'MotifMatch.txt',sep='\t',header=0) + TF_binding,REs1,motifs=list2mat(A,'PositionID','Motif Name','MotifScore') + TF_binding1=pd.DataFrame(TF_binding.T,index=motifs,columns=REs1) + TF_binding1['motif']=motifs + TF_binding1=TF_binding1.merge(Match2,how='inner',left_on='motif',right_on='Motif') + TF_binding=TF_binding1.groupby(['TF'])[REs1].max() + TF_binding=TF_binding.reset_index() + TF_binding.index=TF_binding['TF'] + TF_binding=TF_binding[REs1] + return TF_binding.T + +def cell_type_specific_TF_RE_binding_chr(adata_RNA,adata_ATAC,GRNdir,chrN,genome,celltype,outdir,method,mat): + ## the regions + ## the regions + O_overlap, N_overlap,O_overlap_u,N_overlap_u,O_overlap_hg19_u=load_region(GRNdir,genome,chrN,outdir) + import numpy as np + import pandas as pd + label=adata_RNA.obs['label'].values.tolist() + labelset=list(set(label)) + temp=adata_ATAC.X[np.array(label)==celltype,:].mean(axis=0) + RE=pd.DataFrame(temp.T,index=adata_ATAC.var['gene_ids'].values,columns=['values']) + temp=adata_RNA.X[np.array(label)==celltype,:].mean(axis=0) + TG=pd.DataFrame(temp.T,index=adata_RNA.var['gene_ids'].values,columns=['values']) + del temp +## cell annotation +## extact the overlapped peaks. + RE=RE.loc[N_overlap] + TFbinding=load_TFbinding(GRNdir,O_overlap,O_overlap_u,O_overlap_hg19_u,chrN) + if method=='LINGER': + other_RE=list(set(N_overlap)-set(mat.index)) + if len(other_RE)>0: + B_arr = pd.DataFrame(np.zeros((len(other_RE), mat.shape[1])), columns=mat.columns, index=other_RE) + mat = pd.concat([mat, B_arr]) + mat = mat.loc[N_overlap] + if method=='baseline': + mat=load_TF_RE(GRNdir,chrN,O_overlap,O_overlap_u,O_overlap_hg19_u) + mat.index=N_overlap + TFs = mat.columns + TFoverlap = list(set(TFs) & set(TG.index)) + mat = mat[TFoverlap] + TFbinding = TFbinding[TFoverlap] + TFbinding.index=N_overlap + TF = TG.loc[TFoverlap] + mat_m=np.mean(mat.values[mat>0]) + mat = mat / mat_m + mat.values[mat.values<0]=0 + TFbinding = TFbinding / TFbinding.mean(axis=1).mean() + TF_cluster = TF.values#[:,np.array(label)==celltype].mean(axis=1) + TF_cluster=TF_cluster/TF_cluster.mean() + #TF_cluster = TF_cluster[None,:] + RE_cluster = RE.values#[:,np.array(label)==celltype].mean(axis=1) + RE_cluster=RE_cluster/RE_cluster.mean() + #RE_cluster = RE_cluster[:,None] + S = (np.log(RE_cluster+0.1) + np.log(mat+TFbinding+0.1)).T + np.log(TF_cluster+0.1) + S = np.exp(S.T) + S.index=N_overlap + S_all = S.groupby(S.index).max() + return S_all + + + +def cell_type_specific_TF_RE_binding_score_scNN(mat,TFbinding,RE,TG,TFoverlap): + TF = TG.loc[TFoverlap] + mat_m=np.mean(mat.values[mat>0]) + mat = mat / mat_m + mat.values[mat.values<0]=0 + TFbinding = TFbinding / TFbinding.mean(axis=1).mean() + TF_cluster = TF.values#[:,np.array(label)==celltype].mean(axis=1) + TF_cluster=TF_cluster/TF_cluster.mean() + #TF_cluster = TF_cluster[None,:] + RE_cluster = RE.values#[:,np.array(label)==celltype].mean(axis=1) + RE_cluster=RE_cluster/RE_cluster.mean() + #RE_cluster = RE_cluster[:,None] + S = (np.log(RE_cluster+0.1) + np.log(mat+TFbinding+0.1)).T + np.log(TF_cluster+0.1) + S = np.exp(S.T) + S.index=mat.index + return S + +def cell_type_specific_TF_RE_binding(GRNdir,adata_RNA,adata_ATAC,genome,celltype,outdir,method): + label=adata_RNA.obs['label'].values.tolist() + labelset=list(set(label)) + if (celltype == 'all')&(method!='scNN'): + for label0 in labelset: + print('Generate cell type specitic TF binding potential for cell type '+ str(label0)+'...') + result=pd.DataFrame() + from tqdm import tqdm + for i in tqdm(range(22)): + chrN='chr'+str(i+1) + mat=pd.read_csv(outdir+chrN+'_cell_population_TF_RE_binding.txt',sep='\t',index_col=0,header=0) + out=cell_type_specific_TF_RE_binding_chr(adata_RNA,adata_ATAC,GRNdir,chrN,genome,label0,outdir,method,mat) + #result=pd.concat([result,out],axis=1).fillna(0) + result = pd.concat([result, out], join='outer', axis=0) + chrN='chrX' + mat=pd.read_csv(outdir+chrN+'_cell_population_TF_RE_binding.txt',sep='\t',index_col=0,header=0) + out=cell_type_specific_TF_RE_binding_chr(adata_RNA,adata_ATAC,GRNdir,chrN,genome,label0,outdir,method,mat) + result = pd.concat([result, out], join='outer', axis=0).fillna(0) + result.to_csv(outdir+'cell_type_specific_TF_RE_binding_'+str(label0)+'.txt', sep='\t') + elif method!='scNN': + result=pd.DataFrame() + from tqdm import tqdm + chrom=['chr'+str(i+1) for i in range(22)] + chrom.append('chrX') + for i in tqdm(range(23)): + chrN=chrom[i] + mat=pd.read_csv(outdir+chrN+'_cell_population_TF_RE_binding.txt',sep='\t',index_col=0,header=0) + out=cell_type_specific_TF_RE_binding_chr(adata_RNA,adata_ATAC,GRNdir,chrN,genome,celltype,outdir,method,mat) + #result=pd.concat([result,out],axis=1).fillna(0) + result = pd.concat([result, out], join='outer', axis=0) + result.to_csv(outdir+'cell_type_specific_TF_RE_binding_'+str(celltype)+'.txt', sep='\t') + elif (celltype == 'all')&(method=='scNN'): + A=pd.read_csv(outdir+'cell_population_TF_RE_binding.txt',sep='\t',header=0,index_col=0) + mat,REs,TFs=list2mat(A,'RE','TF','score') + mat=pd.DataFrame(mat,index=REs,columns=TFs) + TFs = mat.columns + TFbinding=load_TFbinding_scNN(GRNdir,outdir,genome) + TG=pd.DataFrame([],index=adata_RNA.var['gene_ids'].values) + TFoverlap = list(set(TFs) & set(TG.index)) + TFoverlap=list(set(TFoverlap) & set(TFbinding.columns)) + mat = mat[TFoverlap] + TFbinding = TFbinding[TFoverlap] + REoverlap=list(set(TFbinding.index)&set(mat.index)) + TFbinding=TFbinding.loc[REoverlap] + TFbinding1=np.zeros((mat.shape[0],len(TFoverlap))) + REidx=pd.DataFrame(range(mat.shape[0]),index=mat.index) + TFbinding1[REidx.loc[TFbinding.index][0].values,:]=TFbinding.values + TFbinding1 = pd.DataFrame(TFbinding1,index=mat.index,columns=TFoverlap) + TFbinding=TFbinding1.copy() + for label0 in labelset: + print('Generate cell type specitic TF binding potential for cell type '+ str(label0)+'...') + from tqdm import tqdm + temp=adata_ATAC.X[np.array(label)==label0,:].mean(axis=0).T + RE=pd.DataFrame(temp,index=adata_ATAC.var['gene_ids'].values,columns=['values']) + temp=adata_RNA.X[np.array(label)==label0,:].mean(axis=0).T + TG=pd.DataFrame(temp,index=adata_RNA.var['gene_ids'].values,columns=['values']) + RE=RE.loc[REs] + result=cell_type_specific_TF_RE_binding_score_scNN(mat,TFbinding,RE,TG,TFoverlap) + result.to_csv(outdir+'cell_type_specific_TF_RE_binding_'+str(label0)+'.txt', sep='\t') + else: + label0=celltype + A=pd.read_csv(outdir+'cell_population_TF_RE_binding.txt',sep='\t',header=0,index_col=0) + mat,REs,TFs=list2mat(A,'RE','TF','score') + mat=pd.DataFrame(mat,index=REs,columns=TFs) + TFs = mat.columns + TFbinding=load_TFbinding_scNN(GRNdir,outdir,genome) + TG=pd.DataFrame([],index=adata_RNA.var['gene_ids'].values) + TFoverlap = list(set(TFs) & set(TG.index)) + TFoverlap=list(set(TFoverlap) & set(TFbinding.columns)) + mat = mat[TFoverlap] + TFbinding = TFbinding[TFoverlap] + REoverlap=list(set(TFbinding.index)&set(RE.index)) + TFbinding=TFbinding.loc[REoverlap] + TFbinding1=np.zeros((mat.shape[0],len(TFoverlap))) + REidx=pd.DataFrame(range(mat.shape[0]),index=mat.index) + TFbinding1[REidx.loc[TFbinding.index][0].values,:]=TFbinding.values + TFbinding1 = pd.DataFrame(TFbinding1,index=mat.index,columns=TFoverlap) + print('Generate cell type specitic TF binding potential for cell type '+ str(label0)+'...') + temp=adata_ATAC.X[np.array(label)==label0,:].mean(axis=0).T + RE=pd.DataFrame(temp,index=adata_ATAC.var['gene_ids'].values,columns=['values']) + temp=adata_RNA.X[np.array(label)==label0,:].mean(axis=0).T + TG=pd.DataFrame(temp,index=adata_RNA.var['gene_ids'].values,columns=['values']) + RE=RE.loc[REs] + result=cell_type_specific_TF_RE_binding_score_scNN(mat,TFbinding,RE,TG,TFoverlap) + result.to_csv(outdir+'cell_type_specific_TF_RE_binding_'+str(label0)+'.txt', sep='\t') + + +def load_shap(chr,outdir): + import torch + import pandas as pd + import numpy as np + import csv + #print('loading shapley value '+chr+' ...') + shap_all=torch.load(outdir+"shap_"+chr+".pt", weights_only=False) + import pandas as pd + idx_file=outdir+'index.txt' + TFName=outdir+'TFName.txt' + #TFE=Input_dir+'TFexp.txt' + REName = 'data/Peaks.txt' +# Open the file in read mode + with open(REName, "r") as file: + # Create a CSV reader + reader = csv.reader(file, delimiter='\t') + # Read the first column and store it in a list + first_column = [row[0] for row in reader] + REName=np.array(first_column) + idx=pd.read_csv(idx_file,sep='\t',header=None) + idx.fillna('', inplace=True) + TFName=pd.read_csv(TFName,sep='\t',header=None) + import numpy as np + idx.columns=['gene','REid','TF_id','REid_b'] + TFName.columns=['Name'] + TFName=TFName['Name'].values + #TFE=pd.read_csv(TFE,header=None,sep='\t') + from scipy.stats import zscore + TFindex=idx['TF_id'].values + REindex=idx['REid'].values + geneName=idx['gene'].values + data_merge=pd.read_csv(outdir+'data_merge.txt',sep='\t',header=0,index_col=0) + data_merge_temp=data_merge[data_merge['chr']==chr] + return data_merge_temp,geneName,REindex,TFindex,shap_all,TFName,REName +def cis_shap(chr,outdir): + RE_2=[] + TG_2=[] + score_2=[] + data_merge_temp,geneName,REindex,TFindex,shap_all,TFName,REName=load_shap(chr,outdir) + from tqdm import tqdm + for j in tqdm(range(data_merge_temp.shape[0]), desc=chr): + ii=data_merge_temp.index[j] + if ii in shap_all.keys(): + AA0=shap_all[ii] + REidxtemp=REindex[ii] + REidxtemp=str(REidxtemp).split('_') + #AA0[:,0:len(TFidxtemp)]=np.multiply(AA0[:,0:len(TFidxtemp)],TFE.values[np.array(TFidxtemp),:].T) + temps=np.abs(AA0).mean(axis=0) + #zscored_arr = zscore(temps) + zscored_arr = np.nan_to_num(temps, nan=0.0) + if (REidxtemp[0]=='') : + REidxtemp=[] + else: + REidxtemp=[int(REidxtemp[i]) for i in range(len(REidxtemp))] + if len(REidxtemp)>0: + REName_temp=REName[np.array(REidxtemp)] + for k in range(len(REidxtemp)): + TG_2.append(geneName[ii]) + RE_2.append(REName_temp[k]) + score_2.append(zscored_arr[k+len(zscored_arr)-len(REidxtemp)]) + RE_TG=pd.DataFrame(TG_2) + RE_TG.columns=['TG'] + RE_TG['RE']=RE_2 + RE_TG['score']=score_2 + RE_TG=RE_TG.groupby(['RE', 'TG'])['score'].max().reset_index() + return RE_TG +def trans_shap(chr,outdir): + TG_1=[] + TF_1=[] + score_1=[] + data_merge_temp,geneName,REindex,TFindex,shap_all,TFName,REName=load_shap(chr,outdir) + from tqdm import tqdm + for j in tqdm(range(data_merge_temp.shape[0]), desc=chr): + ii=data_merge_temp.index[j] + if ii in shap_all.keys(): + AA0=shap_all[ii] + TFidxtemp=TFindex[ii] + TFidxtemp=TFidxtemp.split('_') + TFidxtemp=[int(TFidxtemp[i]) for i in range(len(TFidxtemp))] + TFName_temp=TFName[np.array(TFidxtemp)] + #AA0[:,0:len(TFidxtemp)]=np.multiply(AA0[:,0:len(TFidxtemp)],TFE.values[np.array(TFidxtemp),:].T) + temps=np.abs(AA0).mean(axis=0) + #zscored_arr = zscore(temps) + zscored_arr = np.nan_to_num(temps, nan=0.0) + for k in range(len(TFidxtemp)): + TG_1.append(geneName[ii]) + TF_1.append(TFName_temp[k]) + score_1.append(zscored_arr[k]) + TF_TG=pd.DataFrame(TG_1) + TF_TG.columns=['TG'] + TF_TG['TF']=TF_1 + TF_TG['score']=score_1 + mat,TGs,TFs=list2mat(TF_TG,'TG','TF','score') + mat=pd.DataFrame(mat,index=TGs,columns=TFs) + mat.fillna(0, inplace=True) + return mat + +def load_RE_TG(GRNdir,chrN,O_overlap_u,O_overlap_hg19_u,O_overlap): + #print('load prior RE-TG ...') + from scipy.sparse import coo_matrix + primary_s=pd.read_csv(GRNdir+'Primary_RE_TG_'+chrN+'.txt',sep='\t') + primary_s["RE"] = primary_s["RE"].apply(lambda x: x.split('_')[0]+':'+x.split('_')[1]+'-'+x.split('_')[2]) + primary_s = primary_s[primary_s["RE"].isin(O_overlap_u)] + TGset=primary_s["TG"].unique() + REset=O_overlap_u + # Create a dictionary mapping column names and row names to integer indices + col_dict = {col: i for i, col in enumerate(TGset)} + row_dict = {row: i for i, row in enumerate(REset)} +# Map the column names and row names to integer indices in the DataFrame + primary_s.loc[:,"col_index"] = primary_s["TG"].map(col_dict) + primary_s.loc[:,"row_index"] = primary_s["RE"].map(row_dict) + # Extract the column indices, row indices, and values from the DataFrame + col_indices = primary_s["col_index"].tolist() + row_indices = primary_s["row_index"].tolist() + values = primary_s["score"].tolist() + # Create the sparse matrix using coo_matrix + sparse_S = coo_matrix((values, (row_indices, col_indices))) + sparse_S.colnames = TGset + sparse_S.rownames = REset + array = sparse_S.toarray() + O_overlap_u_df=pd.DataFrame(range(len(O_overlap_u)), index=O_overlap_u) + hg19_38=pd.DataFrame(O_overlap_u, index=O_overlap_hg19_u) + array2=np.zeros([len(O_overlap),array.shape[1]]) + index=O_overlap_u_df.loc[O_overlap][0].values + array2=array[index,:] + array=pd.DataFrame(array2,index=O_overlap,columns=TGset) + return array,TGset +def load_RE_TG_distance(GRNdir,chrN,O_overlap_hg19_u,O_overlap_u,O_overlap,TGoverlap): + #print('load RE-TG distance for '+chrN+'...') + from scipy.sparse import coo_matrix + Dis=pd.read_csv(GRNdir+'RE_TG_distance_'+chrN+'.txt',sep='\t',header=None) + Dis.columns=['RE','TG','dis'] + Dis["RE"] = Dis["RE"].apply(lambda x: x.split('_')[0]+':'+x.split('_')[1]+'-'+x.split('_')[2]) + Dis = Dis[Dis["RE"].isin(O_overlap_hg19_u)] + Dis = Dis[Dis['TG'].isin(TGoverlap)] + col_dict = {col: i for i, col in enumerate(TGoverlap)} + row_dict = {row: i for i, row in enumerate(O_overlap_hg19_u)} +# Map the column names and row names to integer indices in the DataFrame + Dis.loc[:,"col_index"] = Dis["TG"].map(col_dict) + Dis.loc[:,"row_index"] = Dis["RE"].map(row_dict) + col_indices = Dis["col_index"].tolist() + row_indices = Dis["row_index"].tolist() + values = Dis["dis"].tolist() +# Create the sparse matrix using coo_matrix + sparse_dis = coo_matrix((values, (row_indices, col_indices)),shape=(len(O_overlap_u), len(TGoverlap))) + sparse_dis.colnames = TGoverlap + sparse_dis.rownames = O_overlap_u + sparse_dis = sparse_dis.tocsc() + A=sparse_dis.multiply(1 / 25000) + A.data +=0.5 + A.data = np.exp(-A.data) + sparse_dis=A + array = sparse_dis.toarray() + O_overlap_u_df=pd.DataFrame(range(len(O_overlap_u)), index=O_overlap_u) + hg19_38=pd.DataFrame(O_overlap_u, index=O_overlap_hg19_u) + array2=np.zeros([len(O_overlap),array.shape[1]]) + index=O_overlap_u_df.loc[O_overlap][0].values + array2=array[index,:] + array=pd.DataFrame(array2,index=O_overlap,columns=TGoverlap) + return array + + +def load_RE_TG_scNN(outdir): + #print('load prior RE-TG ...') + from scipy.sparse import coo_matrix + import pandas as pd + import numpy as np + dis=pd.read_csv('data/RE_gene_distance.txt',sep='\t',header=0) + dis['distance']=np.exp(-(0.5+dis['distance']/25000)) + REs=dis['RE'].unique() + TGs=dis['gene'].unique() + cis=pd.read_csv(outdir+'cell_population_cis_regulatory.txt',sep='\t',header=None) + cis.columns=['RE','TG','score'] + REs2=cis['RE'].unique() + TGs2=cis['TG'].unique() + REoverlap=list(set(REs2)&set(REs)) + TGoverlap=list(set(TGs)&set(TGs2)) + cisGRN,REs2,TGs2=list2mat_s(cis,REoverlap,TGoverlap,'RE','TG','score') + dis=dis[dis['RE'].isin(REoverlap)] + dis=dis[dis['gene'].isin(TGoverlap)] + distance,REs,TGs=list2mat_s(dis,REoverlap,TGoverlap,'RE','gene','distance') + return distance,cisGRN,REoverlap,TGoverlap + +def cis_reg_chr(GRNdir,adata_RNA,adata_ATAC,genome,chrN,outdir): + import numpy as np + import pandas as pd + from scipy.sparse import csc_matrix + from scipy.sparse import coo_matrix + O_overlap, N_overlap,O_overlap_u,N_overlap_u,O_overlap_hg19_u=load_region(GRNdir,genome,chrN,outdir) + sparse_S,TGset=load_RE_TG(GRNdir,chrN,O_overlap_u,O_overlap_hg19_u,O_overlap) + RE=pd.DataFrame(adata_ATAC.X.toarray().T,index=adata_ATAC.var['gene_ids'].values,columns=adata_ATAC.obs['barcode'].values) +## cell annotation +## extact the overlapped peaks. + RE=RE.loc[N_overlap] + RE=RE.mean(axis=1) + RE=RE/RE.mean()+0.1 + ## select the genes + TG=pd.DataFrame(adata_RNA.X.toarray().T,index=adata_RNA.var['gene_ids'].values,columns=adata_RNA.obs['barcode'].values) + TGoverlap=list(set(TGset)&set(TG.index)) + #target_col_indices = [col_dict[col] for col in TGoverlap] + sparse_S = sparse_S[TGoverlap] + TG=TG.loc[TGoverlap] + TG=TG.mean(axis=1) + TG=TG/TG.mean()+0.1 + sparse_dis=load_RE_TG_distance(GRNdir,chrN,O_overlap_hg19_u,O_overlap_u,O_overlap,TGoverlap) + sparse_S+=0.1 + #Score=csc_matrix(RE).T.multiply(sparse_S.values).multiply(sparse_dis.values).multiply(csc_matrix(TG)).toarray() + Score=np.multiply(sparse_S.values,sparse_dis.values) + Score=pd.DataFrame(Score,index=N_overlap,columns=TGoverlap) + Score=Score.groupby(Score.index).max() + data = Score.values[Score.values!=0] + rows, cols = np.nonzero(Score.values) + coo = coo_matrix((data,(rows,cols)),shape=Score.shape) + combined = np.zeros([len(data),3], dtype=object) + combined[:,0]=Score.index[coo.row] + combined[:,1]=np.array(TGoverlap)[coo.col] + combined[:,2]=coo.data + combined=pd.DataFrame(combined) + return combined + +def cis_shap_scNN(chrtemp,outdir,RE_TGlink1,REName,TFName): + import ast + REName=pd.DataFrame(range(len(REName)),index=REName) + RE_2=[] + TG_2=[] + score_2=[] + shap_all=torch.load(outdir+chrtemp+"_shap"+".pt") + N=RE_TGlink1.shape[0] + for ii in tqdm(range(N)): + AA0=shap_all[ii] + RE_TGlink_temp=RE_TGlink1.values[ii,:] + actual_list = ast.literal_eval(RE_TGlink_temp[1]) + REidxtemp=REName.loc[actual_list].index + TFidxtemp=np.array(range(len(TFName))) + TFidxtemp=TFidxtemp[TFName!=RE_TGlink_temp[0]] + if len(REidxtemp)>0: + temps=np.abs(AA0).mean(axis=0) + zscored_arr = np.nan_to_num(temps, nan=0.0) + for k in range(len(REidxtemp)): + TG_2.append(RE_TGlink_temp[0]) + RE_2.append(REidxtemp[k]) + score_2.append(zscored_arr[k+len(zscored_arr)-len(REidxtemp)]) + RE_TG=pd.DataFrame(TG_2) + RE_TG.columns=['TG'] + RE_TG['RE']=RE_2 + RE_TG['score']=score_2 + RE_TG=RE_TG.groupby(['RE', 'TG'])['score'].max().reset_index() + return RE_TG + + +def cis_reg(GRNdir,adata_RNA,adata_ATAC,genome,method,outdir): + from tqdm import tqdm + chrom=['chr'+str(i+1) for i in range(22)] + chrom.append('chrX') + if method=='baseline': + result=pd.DataFrame([]) + for i in tqdm(range(23)): + chrN=chrom[i] + temp=cis_reg_chr(GRNdir,adata_RNA,adata_ATAC,genome,chrN,outdir) + temp.columns=['RE','TG','Score'] + result=pd.concat([result,temp],axis=0,join='outer') + if method=='LINGER': + import os + from joblib import Parallel, delayed + n_cpus = int(os.environ.get('SLURM_CPUS_PER_TASK', os.cpu_count())) + n_jobs = min(n_cpus, 23) + results = Parallel(n_jobs=n_jobs, backend='loky')( + delayed(cis_shap)(chrN, outdir) + for chrN in chrom + ) + result = pd.concat(results, join='outer', axis=0) + if method=='scNN': + Exp,Opn,Target,RE_TGlink=load_data_scNN(GRNdir,genome) + RE_TGlink=pd.read_csv(outdir+'RE_TGlink.txt',sep='\t',header=0) + RE_TGlink.columns=[0,1,'chr'] + #chrall=[RE_TGlink[0][i][0].split(':')[0] for i in range(RE_TGlink.shape[0])] + chrlist=RE_TGlink['chr'].unique() + REName=Opn.index + geneName=Target.index + TFName=Exp.index + result=pd.DataFrame([]) + for i in tqdm(range(len(chrlist))): + chrN=chrlist[i] + RE_TGlink1=RE_TGlink[RE_TGlink['chr']==chrN] + temp=cis_shap_scNN(chrN,outdir,RE_TGlink1,REName,TFName) + result=pd.concat([result,temp],axis=0,join='outer') + result.to_csv(outdir+'cell_population_cis_regulatory.txt',sep='\t',header=None,index=None) + + +def cell_type_specific_cis_reg_chr(GRNdir,adata_RNA,adata_ATAC,genome,chrN,celltype,outdir): + import numpy as np + import pandas as pd + from scipy.sparse import csc_matrix + from scipy.sparse import coo_matrix + O_overlap, N_overlap,O_overlap_u,N_overlap_u,O_overlap_hg19_u=load_region(GRNdir,genome,chrN,outdir) + sparse_S,TGset=load_RE_TG(GRNdir,chrN,O_overlap_u,O_overlap_hg19_u,O_overlap) + label=adata_RNA.obs['label'].values.tolist() + labelset=list(set(label)) + temp=adata_ATAC.X[np.array(label)==celltype,:].mean(axis=0).T + RE=pd.DataFrame(temp,index=adata_ATAC.var['gene_ids'].values,columns=['values']) + temp=adata_RNA.X[np.array(label)==celltype,:].mean(axis=0).T + TG=pd.DataFrame(temp,index=adata_RNA.var['gene_ids'].values,columns=['values']) + del temp +## cell annotation +## extact the overlapped peaks. + RE=RE.loc[N_overlap] + ## select the genes + TGoverlap=list(set(TGset)&set(TG.index)) + #target_col_indices = [col_dict[col] for col in TGoverlap] + sparse_S = sparse_S[TGoverlap] + TG=TG.loc[TGoverlap] + sparse_dis=load_RE_TG_distance(GRNdir,chrN,O_overlap_hg19_u,O_overlap_u,O_overlap,TGoverlap) + sparse_S+=0.1 +## cell annotation + TG_temp=TG.values#[:,np.array(label)==celltype].mean(axis=1) + TG_temp=TG_temp/TG_temp.mean()+0.1 + RE_temp=RE.values#[:,np.array(label)==celltype].mean(axis=1) + RE_temp=RE_temp/RE_temp.mean()+0.1 + Score=csc_matrix(RE_temp).multiply(sparse_S.values).multiply(sparse_dis.values).multiply(csc_matrix(TG_temp.T)).toarray() + Score=pd.DataFrame(Score,index=N_overlap,columns=TGoverlap) + Score=Score.groupby(Score.index).max() + data = Score.values[Score.values!=0] + rows, cols = np.nonzero(Score.values) + coo = coo_matrix((data,(rows,cols)),shape=Score.shape) + combined = np.zeros([len(data),3], dtype=object) + combined[:,0]=Score.index[coo.row] + combined[:,1]=np.array(TGoverlap)[coo.col] + combined[:,2]=coo.data + resultall=pd.DataFrame(combined) + return resultall +def cell_type_specific_cis_reg_scNN(distance,cisGRN,RE,TG,REs,TGs): + import numpy as np + import pandas as pd + from scipy.sparse import csr_matrix + from scipy.sparse import coo_matrix + RE=RE.loc[REs] + ## select the genes + #target_col_indices = [col_dict[col] for col in TGoverlap] + TG=TG.loc[TGs] + ## cell annotation + TG_temp=TG.values#[:,np.array(label)==celltype].mean(axis=1) + TG_temp=TG_temp/TG_temp.mean()+0.1 + RE_temp=RE.values#[:,np.array(label)==celltype].mean(axis=1) + RE_temp=RE_temp/RE_temp.mean()+0.1 + Score=(cisGRN.multiply(csr_matrix(RE_temp))).multiply(distance).multiply(csr_matrix(TG_temp.T)) + row_indices, col_indices = Score.nonzero() + row_indices=np.array(REs)[row_indices] + col_indices = np.array(TGs)[col_indices] + values = Score.data + combined = np.zeros([len(row_indices),3], dtype=object) + combined[:,0]=row_indices + combined[:,1]=col_indices + combined[:,2]=values + resultall=pd.DataFrame(combined) + return resultall + +def cell_type_specific_cis_reg(GRNdir,adata_RNA,adata_ATAC,genome,celltype,outdir,method): + import pandas as pd + import numpy as np + label=adata_RNA.obs['label'].values.tolist() + labelset=list(set(label)) + chrom=['chr'+str(i+1) for i in range(22)] + chrom.append('chrX') + from tqdm import tqdm + if (celltype=='all')&(method!='scNN'): + for label0 in labelset: + label0=str(label0) + result=pd.DataFrame([]) + for i in tqdm(range(23)): + chrN=chrom[i] + temp=cell_type_specific_cis_reg_chr(GRNdir,adata_RNA,adata_ATAC,genome,chrN,label0,outdir) + result=pd.concat([result,temp],axis=0,join='outer') + chrN='chrX' + temp=cell_type_specific_cis_reg_chr(GRNdir,adata_RNA,adata_ATAC,genome,chrN,label0,outdir) + result=pd.concat([result,temp],axis=0,join='outer') + result.to_csv(outdir+'cell_type_specific_cis_regulatory_'+str(label0)+'.txt',sep='\t',header=None,index=None) + elif (method!='scNN'): + result=pd.DataFrame([]) + for i in tqdm(range(23)): + chrN=chrom[i] + temp=cell_type_specific_cis_reg_chr(GRNdir,adata_RNA,adata_ATAC,genome,chrN,celltype,outdir) + result=pd.concat([result,temp],axis=0,join='outer') + chrN='chrX' + temp=cell_type_specific_cis_reg_chr(GRNdir,adata_RNA,adata_ATAC,genome,chrN,celltype,outdir) + result=pd.concat([result,temp],axis=0,join='outer') + result.to_csv(outdir+'cell_type_specific_cis_regulatory_'+celltype+'.txt',sep='\t',header=None,index=None) + elif (celltype=='all')&(method=='scNN'): + distance,cisGRN,REs,TGs=load_RE_TG_scNN(outdir) + for label0 in labelset: + label0=str(label0) + temp=adata_ATAC.X[np.array(label)==label0,:].mean(axis=0).T + RE=pd.DataFrame(temp,index=adata_ATAC.var['gene_ids'].values,columns=['values']) + temp=adata_RNA.X[np.array(label)==label0,:].mean(axis=0).T + TG=pd.DataFrame(temp,index=adata_RNA.var['gene_ids'].values,columns=['values']) + del temp + + result=cell_type_specific_cis_reg_scNN(distance,cisGRN,RE,TG,REs,TGs) + result.to_csv(outdir+'cell_type_specific_cis_regulatory_'+label0+'.txt',sep='\t',header=None,index=None) + else: + label0=celltype + label0=str(label0) + temp=adata_ATAC.X[np.array(label)==label0,:].mean(axis=0).T + RE=pd.DataFrame(temp,index=adata_ATAC.var['gene_ids'].values,columns=['values']) + temp=adata_RNA.X[np.array(label)==label0,:].mean(axis=0).T + TG=pd.DataFrame(temp,index=adata_RNA.var['gene_ids'].values,columns=['values']) + del temp + result=cell_type_specific_cis_reg_scNN(distance,cisGRN,RE,TG,REs,TGs) + result.to_csv(outdir+'cell_type_specific_cis_regulatory_'+label0+'.txt',sep='\t',header=None,index=None) + +def trans_shap_scNN(chrtemp,outdir,RE_TGlink1,REName,TFName): + import ast + TG_1=[] + TF_1=[] + score_1=[] + REName=pd.DataFrame(range(len(REName)),index=REName) + shap_all=torch.load(outdir+chrtemp+"_shap"+".pt") + N=RE_TGlink1.shape[0] + from tqdm import tqdm + for ii in range(N): + AA0=shap_all[ii] + RE_TGlink_temp=RE_TGlink1.values[ii,:] + actual_list = ast.literal_eval(RE_TGlink_temp[1]) + REidxtemp=REName.loc[actual_list].index + TFidxtemp=np.array(range(len(TFName))) + TFidxtemp=TFidxtemp[TFName!=RE_TGlink_temp[0]] + temps=np.abs(AA0).mean(axis=0) + #zscored_arr = zscore(temps) + zscored_arr = np.nan_to_num(temps, nan=0.0) + for k in range(len(TFidxtemp)): + TG_1.append(RE_TGlink_temp[0]) + TF_1.append(TFName[TFidxtemp[k]]) + score_1.append(zscored_arr[k]) + TF_TG=pd.DataFrame(TG_1) + TF_TG.columns=['TG'] + TF_TG['TF']=TF_1 + TF_TG['score']=score_1 + mat,TGs,TFs=list2mat(TF_TG,'TG','TF','score') + mat=pd.DataFrame(mat,index=TGs,columns=TFs) + mat.fillna(0, inplace=True) + return mat + + +def load_cis(Binding,celltype,outdir): + from scipy.sparse import coo_matrix + import pandas as pd + import numpy as np + if celltype=='': + cis=pd.read_csv(outdir+'cell_population_cis_regulatory.txt',sep='\t',header=None) + else: + cis=pd.read_csv(outdir+'cell_type_specific_cis_regulatory_'+celltype+'.txt',sep='\t',header=None) + cis.columns=['RE','TG','Score'] + TGset=cis['TG'].unique() + REset=Binding.index + TFset=Binding.columns + col_dict = {col: i for i, col in enumerate(TGset)} + row_dict = {row: i for i, row in enumerate(REset)} + cis=cis[cis["RE"].isin(REset)] +# Map the column names and row names to integer indices in the DataFrame + cis["col_index"] = cis["TG"].map(col_dict) + cis["row_index"] = cis["RE"].map(row_dict) + # Extract the column indices, row indices, and values from the DataFrame + col_indices = cis["col_index"].tolist() + row_indices = cis["row_index"].tolist() + values = cis["Score"].tolist() + # Create the sparse matrix using coo_matrix + sparse_S = coo_matrix((values, (row_indices, col_indices)),shape=(len(REset), len(TGset))) + sparse_S.colnames = TGset + sparse_S.rownames = REset + cis=sparse_S.toarray() + cis=pd.DataFrame(cis,index=REset,columns=TGset) + return cis + +def load_TF_TG( GRNdir, TFset,TGset): + TF_TG_all=np.zeros([len(TGset),len(TFset)]) + a=list(range(1,23)) + a.append('X') + for i in a: + chrN='chr'+str(i) + TF_TG = pd.read_csv(GRNdir+'Primary_TF_TG_'+chrN+'.txt',sep='\t') + TF_TG = TF_TG[TF_TG['TF'].isin(TFset)] + TF_TG = TF_TG[TF_TG['TG'].isin(TGset)] + col_dict = {col: i for i, col in enumerate(TFset)} + row_dict = {row: i for i, row in enumerate(TGset)} + TF_TG["col_index"] = TF_TG["TF"].map(col_dict) + TF_TG["row_index"] = TF_TG["TG"].map(row_dict) + col_indices = TF_TG["col_index"].tolist() + row_indices = TF_TG["row_index"].tolist() + values = TF_TG["score"].tolist() + sparse_S = coo_matrix((values, (row_indices, col_indices)),shape=( len(TGset),len(TFset))) + idx=list(set(row_indices)) + TGset1=TGset[idx] + TF_TG=sparse_S.toarray() + TF_TG=pd.DataFrame(TF_TG,index=TGset,columns=TFset) + TF_TG=TF_TG.loc[TGset1] + TF_TG_all[idx,:]=TF_TG.values + TF_TG_all=pd.DataFrame(TF_TG_all,index=TGset,columns=TFset) + return TF_TG_all + +def trans_reg(GRNdir,method,outdir,genome): + import ast + import pandas as pd + from scipy.sparse import coo_matrix + import numpy as np + import pandas as pd + from scipy.sparse import csc_matrix + from scipy.sparse import coo_matrix + print('Generate trans-regulatory netowrk ...') + if method=='baseline': + Binding=pd.read_csv(outdir+'cell_population_TF_RE_binding.txt',sep='\t',index_col=0) + cis=load_cis(Binding,'',outdir) + TFset=Binding.columns + TGset=cis.columns + TF_TG=load_TF_TG(GRNdir, TFset,TGset) + S=np.matmul(Binding.values.T, cis.values).T*(TF_TG.values.T).T + S=pd.DataFrame(S, index=TGset,columns=TFset) + elif method=='LINGER': + chrom=['chr'+str(i+1) for i in range(22)] + chrom.append('chrX') + import os + from joblib import Parallel, delayed + n_cpus = int(os.environ.get('SLURM_CPUS_PER_TASK', os.cpu_count())) + n_jobs = min(n_cpus, 23) + results = Parallel(n_jobs=n_jobs, backend='loky')( + delayed(trans_shap)(chrN, outdir) + for chrN in chrom + ) + S = pd.concat(results, join='outer', axis=0) + elif method=='scNN': + Exp,Opn,Target,RE_TGlink=load_data_scNN(GRNdir,genome) + RE_TGlink=pd.read_csv(outdir+'RE_TGlink.txt',sep='\t',header=0) + RE_TGlink.columns=[0,1,'chr'] + #chrall=[RE_TGlink[0][i][0].split(':')[0] for i in range(RE_TGlink.shape[0])] + chrlist=RE_TGlink['chr'].unique() + REName=Opn.index + geneName=Target.index + TFName=Exp.index + result=pd.DataFrame([]) + S=pd.DataFrame([]) + for i in tqdm(range(len(chrlist))): + chrN=chrlist[i] + RE_TGlink1=RE_TGlink[RE_TGlink['chr']==chrN] + temp=trans_shap_scNN(chrN,outdir,RE_TGlink1,REName,TFName) + S=pd.concat([S,temp],axis=0,join='outer') + print('Save trans-regulatory netowrk ...') + S.to_csv(outdir+'cell_population_trans_regulatory.txt',sep='\t') + +def cell_type_specific_trans_reg(GRNdir,adata_RNA,celltype,outdir): + import pandas as pd + import numpy as np + from scipy.sparse import csc_matrix + from scipy.sparse import coo_matrix + label=adata_RNA.obs['label'].values.tolist() + labelset=list(set(label)) + if celltype=='all': + for label0 in labelset: + Binding=pd.read_csv(outdir+'cell_type_specific_TF_RE_binding_'+str(label0)+'.txt',sep='\t',index_col=0) + label0=str(label0) + cis=load_cis(Binding,label0,outdir) + TFset=Binding.columns + TGset=cis.columns + #TF_TG=load_TF_TG(GRNdir, TFset,TGset) + S=np.matmul(Binding.values.T, cis.values).T#*(TF_TG.values.T).T + S=pd.DataFrame(S, index=TGset,columns=TFset) + S.to_csv(outdir+'cell_type_specific_trans_regulatory_'+str(label0)+'.txt',sep='\t') + else: + Binding=pd.read_csv(outdir+'cell_type_specific_TF_RE_binding_'+celltype+'.txt',sep='\t',index_col=0) + cis=load_cis(Binding,celltype,outdir) + TFset=Binding.columns + TGset=cis.columns + #TF_TG=load_TF_TG(GRNdir, TFset,TGset) + S=np.matmul(Binding.values.T, cis.values).T#*(TF_TG.values.T).T + S=pd.DataFrame(S, index=TGset,columns=TFset) + S.to_csv(outdir+'cell_type_specific_trans_regulatory_'+celltype+'.txt',sep='\t') \ No newline at end of file diff --git a/workflow/scripts/mth/linger/grn.py b/workflow/scripts/mth/linger/grn.py new file mode 100644 index 00000000..c3cd47da --- /dev/null +++ b/workflow/scripts/mth/linger/grn.py @@ -0,0 +1,41 @@ +import argparse +import pandas as pd + +parser = argparse.ArgumentParser() +parser.add_argument('-d','--out_dir', required=True) +parser.add_argument('-o','--path_out', required=True) +args = vars(parser.parse_args()) + +out_dir = args['out_dir'] +path_out = args['path_out'] + +# Load three cell population GRNs +alpha = pd.read_csv(f'{out_dir}/cell_population_TF_RE_binding.txt', sep='\t', index_col=0) +beta = pd.read_csv(f'{out_dir}/cell_population_cis_regulatory.txt', sep='\t', header=None) +gamma = pd.read_csv(f'{out_dir}/cell_population_trans_regulatory.txt', sep='\t', index_col=0) +beta.columns = ['RE', 'TG', 'beta'] + +alpha.index.name = 'RE' +alpha_long = alpha.reset_index().melt(id_vars='RE', var_name='TF', value_name='alpha') + +gamma.index.name = 'TG' +gamma_long = gamma.reset_index().melt(id_vars='TG', var_name='TF', value_name='gamma') + +# Filter to keep top 10% (reduce the combinatorial explosion of later steps) +alpha_long = alpha_long[alpha_long['alpha'] > alpha_long['alpha'].quantile(0.90)] +beta = beta[beta['beta'] > beta['beta'].quantile(0.90)] +gamma_long = gamma_long[gamma_long['gamma'] > gamma_long['gamma'].quantile(0.90)] + +# Build triplets +triplets = pd.merge(alpha_long, beta, on='RE') +triplets = pd.merge(triplets, gamma_long, on=['TF', 'TG']) + +# Compute score +triplets['score'] = triplets['alpha'] * triplets['beta'] * triplets['gamma'] +triplets = triplets[['TF', 'RE', 'TG', 'score']].sort_values('score', ascending=False) + +# Reformat +triplets.columns = ['source', 'cre', 'target', 'score'] +triplets['cre'] = triplets['cre'].str.replace(':', '-') + +triplets.to_csv(path_out, index=False) \ No newline at end of file diff --git a/workflow/scripts/mth/linger/linger.py b/workflow/scripts/mth/linger/linger.py new file mode 100644 index 00000000..41e1d870 --- /dev/null +++ b/workflow/scripts/mth/linger/linger.py @@ -0,0 +1,151 @@ +import argparse +import mudata as mu +import pandas as pd +import time +import os + +from LingerGRN.pseudo_bulk import * + +def log(msg): + print(f"[{time.strftime('%H:%M:%S')}] {msg}", flush=True) + +# Init args +parser = argparse.ArgumentParser() +parser.add_argument('-g','--linger_GRN', required=True) +parser.add_argument('-d','--out_dir', required=True) +parser.add_argument('-m','--path_mdata', required=True) +parser.add_argument('-v','--version', required=True) +parser.add_argument('-gen','--genome', required=True) +parser.add_argument('-md','--mode', required=True) +args = vars(parser.parse_args()) + +GRNdir = args['linger_GRN'] + "/" +out_dir = args['out_dir'] + "/" +path_mdata = args['path_mdata'] +version = args['version'] +genome = args['genome'] +mode = args['mode'] + +if mode == 'parallel': + from preprocess_fast import * + import LINGER_tr_fast as LINGER_tr + import LL_net_fast as LL_net +else: + from LingerGRN.preprocess import * + import LingerGRN.LINGER_tr as LINGER_tr + import LingerGRN.LL_net as LL_net + +log(f"Using {mode} mode") + +# Read mdata +mdata = mu.read(path_mdata) +adata_RNA = mdata["rna"].copy() +adata_ATAC = mdata["atac"].copy() + +# Restore raw counts (pseudo_bulk expects raw counts matrices) +adata_RNA.X = adata_RNA.layers['counts'].copy() +adata_ATAC.X = adata_ATAC.layers['counts'].copy() +log(f"adata_RNA loaded:\n{adata_RNA}\n") +log(f"adata_ATAC loaded:\n{adata_ATAC}\n") + +# Add barcode, gene_ids, label columns +adata_RNA.obs['barcode'] = adata_RNA.obs_names +adata_ATAC.obs['barcode'] = adata_ATAC.obs_names +adata_RNA.var['gene_ids'] = adata_RNA.var_names +adata_ATAC.var['gene_ids'] = adata_ATAC.var_names.str.replace('-', ':', n=1) +adata_RNA.obs['label'] = mdata.obs.loc[adata_RNA.obs_names, 'celltype'].values +adata_ATAC.obs['label'] = mdata.obs.loc[adata_ATAC.obs_names, 'celltype'].values + +# Add sample column +bc = adata_RNA.obs['barcode'].values +if len(bc[0].split("-")) == 2: + adata_RNA.obs['sample'] = [int(s.split("-")[1]) for s in bc] + adata_ATAC.obs['sample'] = [int(s.split("-")[1]) for s in bc] +else: + adata_RNA.obs['sample'] = 1 + adata_ATAC.obs['sample'] = 1 + +adata_RNA.var_names_make_unique() +adata_RNA.var['gene_ids'] = adata_RNA.var_names.values + +# Intersect barcodes +shared = adata_RNA.obs_names.intersection(adata_ATAC.obs_names) +adata_RNA = adata_RNA[shared].copy() +adata_ATAC = adata_ATAC[shared].copy() +log(f"adata_RNA for Linger:\n{adata_RNA}\n") +log(f"adata_ATAC for Linger:\n{adata_ATAC}\n") + + +# Generate pseudo-bulk +log("Generating pseudo-bulk / metacells...") +samplelist = list(set(adata_ATAC.obs['sample'].values)) +TG_pseudobulk = pd.DataFrame([]) # TG x metacells +RE_pseudobulk = pd.DataFrame([]) # RE x metacells +singlepseudobulk = adata_RNA.obs['sample'].nunique() > 10 + +for tempsample in samplelist: + adata_RNAtemp = adata_RNA[adata_RNA.obs['sample'] == tempsample] + adata_ATACtemp = adata_ATAC[adata_ATAC.obs['sample'] == tempsample] + TG_temp, RE_temp = pseudo_bulk(adata_RNAtemp, adata_ATACtemp, singlepseudobulk) + TG_pseudobulk = pd.concat([TG_pseudobulk, TG_temp], axis=1) + RE_pseudobulk = pd.concat([RE_pseudobulk, RE_temp], axis=1) + RE_pseudobulk[RE_pseudobulk > 100] = 100 + + +# Save pseudobulk and adata matrices +log("Saving pseudobulk data...") +os.makedirs(out_dir + 'data/', exist_ok=True) +adata_ATAC.write(out_dir + 'data/adata_ATAC.h5ad') +adata_RNA.write(out_dir + 'data/adata_RNA.h5ad') +TG_pseudobulk.fillna(0).to_csv(out_dir + 'data/TG_pseudobulk.tsv') +RE_pseudobulk.fillna(0).to_csv(out_dir + 'data/RE_pseudobulk.tsv') +pd.DataFrame(adata_ATAC.var['gene_ids']).to_csv(out_dir + 'data/Peaks.txt', header=None, index=None) + + +# Train the model +log("Preprocessing and training LINGER model...") + + +"""" +# LINGER assumes following dir strucutre: +. +├── LINGER_data/ +│ └── data_bulk/ +│ +├── LINGER_output/ +│ ├── cell_population_TF_RE_binding.txt +│ ├── cell_population_cis_regulatory.txt +│ ├── cell_population_trans_regulatory.txt +│ └── ... +│ +└── data/ + ├── Peaks.txt + ├── TG_pseudobulk.tsv + ├── RE_pseudobulk.tsv + ├── adata_RNA.h5ad + └── adata_ATAC.h5ad +""" + +# LINGER_data/ is dbs/lingerGRN/data_bulk (GRNdir) +# LINGER_output/ is dts/{org}/{dat}/cases/{case}/runs/linger/ (out_dir) +# ├── data/ is out_dir + data/ (to keep files together) + +# since data/ is not at top level anymore, we need to chdir to out_dir +GRNdir = os.path.abspath(GRNdir) + "/" +out_dir = os.path.abspath(out_dir) + "/" +os.chdir(out_dir) + +print(GRNdir) +print(out_dir) + +preprocess(TG_pseudobulk, RE_pseudobulk, GRNdir, genome, version, out_dir) + +activef='ReLU' +LINGER_tr.training(GRNdir, version, out_dir, activef, species='Human') + +# Generate regulatory networks +log("Generating cell population GRNs...") +LL_net.TF_RE_binding(GRNdir, adata_RNA, adata_ATAC, genome, version, out_dir) +LL_net.cis_reg(GRNdir, adata_RNA, adata_ATAC, genome, version, out_dir) +LL_net.trans_reg(GRNdir, version, out_dir, genome) +log("GRNs generation done") diff --git a/workflow/scripts/mth/linger/linger.sh b/workflow/scripts/mth/linger/linger.sh new file mode 100644 index 00000000..a50738a0 --- /dev/null +++ b/workflow/scripts/mth/linger/linger.sh @@ -0,0 +1,27 @@ +#!/bin/bash + +while [[ "$#" -gt 0 ]]; do + case $1 in + --linger_GRN) linger_GRN="$2"; shift ;; + --out_dir) out_dir="$2"; shift ;; + --path_mdata) path_mdata="$2"; shift ;; + --version) version="$2"; shift ;; + --genome) genome="$2"; shift ;; + --mode) mode="$2"; shift ;; + --path_out) path_out="$2"; shift ;; + *) echo echo "Unknown parameter passed: $1"; exit 1 ;; + esac + shift +done + +python workflow/scripts/mth/linger/linger.py \ +-g $linger_GRN \ +-d $out_dir \ +-m $path_mdata \ +-v $version \ +-gen $genome \ +-md $mode + +python workflow/scripts/mth/linger/grn.py \ +-d $out_dir \ +-o $path_out \ No newline at end of file diff --git a/workflow/scripts/mth/linger/preprocess_fast.py b/workflow/scripts/mth/linger/preprocess_fast.py new file mode 100644 index 00000000..60531245 --- /dev/null +++ b/workflow/scripts/mth/linger/preprocess_fast.py @@ -0,0 +1,349 @@ +import os +import numpy as np +import pandas as pd +#from LingerGRN.immupute_dis import immupute_dis +#import LingerGRN.pseudo_bulk as pseudo_bulk +import subprocess +from tqdm import tqdm +from joblib import Parallel, delayed + +import warnings +warnings.filterwarnings("ignore", category=FutureWarning, module="anndata") + +def list2mat(df,i_n,j_n,x_n): + TFs = df[j_n].unique() + REs = df[i_n].unique() +#Initialize matrix as numpy array +#Map row and col indices for lookup + row_map = {r:i for i,r in enumerate(REs)} + col_map = {c:i for i,c in enumerate(TFs)} + row_indices = np.array([row_map[row] for row in df[i_n]]) + col_indices = np.array([col_map[col] for col in df[j_n]]) + from scipy.sparse import coo_matrix + matrix = coo_matrix((df[x_n], (row_indices, col_indices)), shape=(len(REs), len(TFs))) + mat=coo_matrix.toarray(matrix) + return mat,REs,TFs + + + +def gene_expression(GRNdir,TG_pseudobulk,outdir): + gene = pd.read_csv(GRNdir+'bulk_gene_all.txt') + gene.columns=['gene'] + #gene=gene['gene'] + d1 = np.isin(TG_pseudobulk.index, gene['gene'].values) + List = TG_pseudobulk.index[d1] + A = np.log2(1 + TG_pseudobulk.loc[List]) + #Write Exp.txt and Symbol.txt + pd.DataFrame(A).to_csv(outdir+'Exp.txt',sep='\t',index=False,header=False) + pd.DataFrame(List).to_csv(outdir+'Symbol.txt', sep='\t', header=False, index=False) + pd.DataFrame(A.columns).to_csv(outdir+'Col.txt',sep='\t',index=False,header=False) + return List,A + + +def TF_expression(TFName,List,Match2,A,outdir): + d= np.isin(TFName,List) + TFName = TFName[d] + List_idx=pd.DataFrame(range(len(List)),index=List) + f=List_idx.loc[TFName][0].values + TF = A.values[f, :] + Match2 = Match2[np.isin(Match2[:, 1], TFName)] + d = np.isin(TFName, Match2[:, 1]) + TFName = TFName[d] + TF = TF[d, :] + pd.DataFrame(TF).to_csv(outdir+'TFexp.txt', sep='\t', header=False, index=False) + pd.DataFrame(TFName).to_csv(outdir+'TFName.txt', sep='\t', header=False, index=False) + return TFName + +def index_generate(choosL_i,merged_s,merged_b,TFName): + if choosL_i in merged_s.index: + REid = merged_s.loc[choosL_i]['id_s'] + REid_b = merged_b.loc[choosL_i]['id_b'] + else: + REid='' + REid_b='' + TFName_1 = np.delete(TFName, np.where(TFName == choosL_i)) + TFid = np.where(np.isin(TFName, TFName_1))[0] + RE_s = '_'.join(map(str, REid)) + TF_s = '_'.join(map(str, TFid)) + RE_b = '_'.join(map(str, REid_b)) + return choosL_i, RE_s, TF_s, RE_b + +def load_corr_RE_TG(List,Element_name,Element_name_bulk,outdir): + Element_gene = pd.read_csv(outdir+"hg19_Peak_hg19_gene_u.txt", delimiter="\t", header=None) + choosL = List + index_ElementName=pd.DataFrame(np.arange(0, len(Element_name)),index=Element_name) + index_Element_name_bulk=pd.DataFrame(np.arange(0, len(Element_name_bulk)),index=Element_name_bulk) + index_Element_name_bulk=index_Element_name_bulk.groupby(index_Element_name_bulk.index).min() + Element_gene.columns=['Element_name_b','Element_name_s','TG'] + Element_gene['value']=1 + #RE_all_s=Element_gene['Element_name_s'].unique() + #RE_all_b=Element_gene['Element_name_b'].unique() + #index_RE_all_s=pd.DataFrame(np.arange(0, len(RE_all_s)),index=RE_all_s) + #index_RE_all_b=pd.DataFrame(np.arange(0, len(RE_all_b)),index=RE_all_b) + Element_gene['id_s']=index_ElementName.loc[Element_gene['Element_name_s']][0].values + Element_gene['id_b']=index_Element_name_bulk.loc[Element_gene['Element_name_b']][0].values + merged_s = Element_gene.groupby('TG')['id_s'].agg(list).reset_index() + merged_b = Element_gene.groupby('TG')['id_b'].agg(list).reset_index() + merged_s = merged_s.set_index('TG') + merged_b =merged_b.set_index('TG') + #index_ElementName1=index_ElementName.loc[RE_all_s][0].values + #index_Element_name_bulk1=index_Element_name_bulk.loc[RE_all_b][0].values + return merged_s,merged_b + +def load_motifbinding_chr(chrN,GRNdir,motifWeight,outdir): + Motif_binding_temp=pd.read_csv(GRNdir+'MotifTarget_Matrix_'+chrN+'.txt',sep='\t',index_col=0) + REs=Motif_binding_temp.index + march_hg19_Regrion=pd.read_csv(outdir+'MotifTarget_hg19_hg38_'+chrN+'.txt',sep='\t',header=None) + REoverlap=list(set(march_hg19_Regrion[1].values)) + Motif_binding_temp1=Motif_binding_temp.loc[REoverlap] + REs=Motif_binding_temp1.index + Motif_binding_temp=np.zeros([march_hg19_Regrion.shape[0],Motif_binding_temp.shape[1]]) + Motif_binding_temp=Motif_binding_temp1.loc[march_hg19_Regrion[1].values].values + Motif_binding_temp=pd.DataFrame(Motif_binding_temp,index=march_hg19_Regrion[0].values,columns=Motif_binding_temp1.columns) + Motif_binding_temp1=Motif_binding_temp.groupby(Motif_binding_temp.index).max() + motifoverlap=list(set(Motif_binding_temp1.columns)&set(motifWeight.index)) + Motif_binding_temp1=Motif_binding_temp1[motifoverlap] + motifWeight=motifWeight.loc[Motif_binding_temp1.columns] + Motif_binding = np.diag(1.0 / (motifWeight.T + 0.1)) * Motif_binding_temp1.values.T + Motif_binding = np.log1p(Motif_binding) + return Motif_binding_temp1 + #return Motif_binding # why not ? + +def load_TFbinding(GRNdir,motifWeight,Match2,TFName,Element_name,outdir): + chrall=['chr'+str(i+1) for i in range(22)] + chrall.append('chrX') + n_jobs = 4 + results = Parallel(n_jobs=n_jobs, backend='loky', verbose=10)( + delayed(load_motifbinding_chr)(chrN, GRNdir, motifWeight, outdir) + for chrN in chrall + ) + motif_binding = pd.concat(results, join='outer', axis=0) + motif_binding=motif_binding.fillna(0) + motif_binding=motif_binding.groupby(motif_binding.index).max() + motifoverlap=list(set(motif_binding.columns)&set(motifWeight.index)) + Match2=Match2[np.isin(Match2[:, 0],motifoverlap), :] + TF_binding_temp = np.zeros((len(TFName), len(Element_name))) + #Motif_binding=np.zeros((motif_binding.shape[1], len(Element_name))) + Element_name_idx=pd.DataFrame(range(len(Element_name)),index=Element_name) + idx=Element_name_idx.loc[motif_binding.index][0].values + Motif_binding=np.zeros((motif_binding.shape[1], len(Element_name))) + Motif_binding[:,idx]=motif_binding.loc[Element_name[idx]].values.T + Motif_binding=pd.DataFrame(Motif_binding,index=motif_binding.columns,columns=Element_name) + Match2=Match2[np.isin(Match2[:, 1],TFName), :] + Motif_binding=Motif_binding.loc[Match2[:, 0]] + Motif_binding.index=Match2[:, 1] + TF_binding=Motif_binding.groupby(Motif_binding.index).sum() + a = np.sum(TF_binding.values, axis=1) + a[a == 0] =1 + TF_binding_n = np.diag(1.0 / a) @TF_binding.values + TF_binding_n=pd.DataFrame(TF_binding_n.T,index=Element_name,columns=TF_binding.index) + TF_binding=np.zeros((len(Element_name),len(TFName))) + idx=np.isin(TFName,TF_binding_n.columns) + TF_binding[:,idx]=TF_binding_n[TFName[idx]].values + TF_binding=pd.DataFrame(TF_binding,index=Element_name,columns=TFName) + TF_binding.to_csv(outdir+'TF_binding.txt',sep='\t',index=None,header=None) + +def _process_chr_overlap(chrtemp, GRNdir, genome, outdir): + import pybedtools, os + a = pybedtools.example_bedtool(outdir+'match_hg19_peak.bed') + b = pybedtools.example_bedtool(GRNdir+'MotifTarget_matrix_'+chrtemp+'.bed') + a_with_b = a.intersect(b, wa=True,wb=True) + a_with_b.saveas(outdir+'temp_'+chrtemp+'.bed') + a_with_b = pd.read_csv(outdir+'temp_'+chrtemp+'.bed', sep='\t', header=None) + a_with_b=a_with_b[(a_with_b[1].values==a_with_b[7].values)&(a_with_b[2].values==a_with_b[8].values)] + a_with_b_n = pd.DataFrame({ + 'column1': a_with_b[3] + ':' + a_with_b[4].astype(str) + '-' + a_with_b[5].astype(str), + 'column2': a_with_b[6] + ':' + a_with_b[7].astype(str) + '-' + a_with_b[8].astype(str)}) + a_with_b_n=a_with_b_n.drop_duplicates() + a_with_b_n.to_csv(outdir+'MotifTarget_hg19_hg38_'+chrtemp+'.txt',sep='\t',header=None,index=None) + a = pybedtools.example_bedtool(GRNdir+genome+'_Peaks_'+chrtemp+'.bed') + b = pybedtools.example_bedtool(outdir+'Region.bed') + a_with_b = a.intersect(b, wa=True,wb=True) + a_with_b.saveas(outdir+'Region_overlap_'+chrtemp+'.bed') + os.remove(outdir+'temp_'+chrtemp+'.bed') + +def extract_overlap_regions(genome,GRNdir,outdir,method): + import pybedtools + import pandas as pd + import os + os.makedirs(outdir, exist_ok=True) + input_file = 'data/Peaks.txt' + output_file = outdir+'Region.bed' +# Read the input file + df = pd.read_csv(input_file, sep='\t',header=None) + chromosomes = [item.split(':')[0] for item in df[0].values] +# Drop the first row +# Replace ':' and '-' with tabs + df = df.replace({':': '\t', '-': '\t'}, regex=True) + chrall=['chr'+str(i+1) for i in range(23)]+['chrX'] + df=df[pd.DataFrame(chromosomes)[0].isin(chrall).values] + df.to_csv(output_file, index=None, header=None) + if method=='LINGER': + if genome=='hg38': + a = pybedtools.example_bedtool(outdir+'Region.bed') + b = pybedtools.example_bedtool(GRNdir+'hg38_hg19_pair.bed') + a_with_b = a.intersect(b, wa=True,wb=True) + a_with_b.saveas(outdir+'temp.bed') + a_with_b=pd.read_csv(outdir+'temp.bed',sep='\t',header=None) + a_with_b[[6,7,8,0,1,2]].to_csv(outdir+'match_hg19_peak.bed',sep='\t',header=None,index=None) + if genome=='hg19': + a = pybedtools.example_bedtool(outdir+'Region.bed') + b = pybedtools.example_bedtool(GRNdir+'hg19_hg38_pair.bed') + a_with_b = a.intersect(b, wa=True,wb=True) + a_with_b.saveas(outdir+'temp.bed') + a_with_b=pd.read_csv(outdir+'temp.bed',sep='\t',header=None) + a_with_b[[6,7,8,0,1,2]].to_csv(outdir+'match_hg19_peak.bed',sep='\t',header=None,index=None) + a = pybedtools.example_bedtool(outdir+'match_hg19_peak.bed') + b = pybedtools.example_bedtool(GRNdir+'RE_gene_corr_hg19.bed') + a_with_b = a.intersect(b, wa=True,wb=True) + a_with_b.saveas(outdir+'temp.bed') + a_with_b=pd.read_csv(outdir+'temp.bed',sep='\t',header=None) + a_with_b=a_with_b[(a_with_b[1].values==a_with_b[7].values)&(a_with_b[2].values==a_with_b[8].values)] + a_with_b_n = pd.DataFrame({ + 'column1': a_with_b[0] + ':' + a_with_b[1].astype(str) + '-' + a_with_b[2].astype(str), + 'column2': a_with_b[3] + ':' + a_with_b[4].astype(str) + '-' + a_with_b[5].astype(str), + 'column3': a_with_b[9]}) + a_with_b_n=a_with_b_n.drop_duplicates() + a_with_b_n.to_csv(outdir+'hg19_Peak_hg19_gene_u.txt',sep='\t',header=None,index=None) + chr_all=['chr'+str(i+1) for i in range(22)] + chr_all.append('chrX') + + n_jobs = 8 + Parallel(n_jobs=n_jobs, backend='loky', verbose=10)( + delayed(_process_chr_overlap)(chrtemp, GRNdir, genome, outdir) + for chrtemp in chr_all + ) + + # cleanup files in /tmp + pybedtools.cleanup(remove_all=True) + + if method=='baseline': + chr_all=['chr'+str(i+1) for i in range(22)] + chr_all.append('chrX') + for chrtemp in chr_all: + a = pybedtools.example_bedtool(GRNdir+genome+'_Peaks_'+chrtemp+'.bed') + b = pybedtools.example_bedtool(outdir+'Region.bed') + a_with_b = a.intersect(b, wa=True,wb=True) + a_with_b.saveas(outdir+'Region_overlap_'+chrtemp+'.bed') + + +def preprocess(TG_pseudobulk,RE_pseudobulk,GRNdir,genome,method,outdir): + #package_dir = os.path.dirname(os.path.abspath(__file__)) + if method=='LINGER': + print('Overlapping regions...') + extract_overlap_regions(genome,GRNdir,outdir,method) + print('Mapping gene expression...') + TFName = pd.read_csv(GRNdir+'TFName.txt',header=None) + TFName.columns=['TFName'] + TFName=TFName['TFName'].values + Match2=pd.read_csv(GRNdir+'Match2.txt',sep='\t') + Match2=Match2.values + List,A=gene_expression(GRNdir,TG_pseudobulk,outdir) + print('Generate TF expression...') + TFName=TF_expression(TFName,List,Match2,A,outdir) + print('Generate RE chromatin accessibility...') + RE_pseudobulk.to_csv(outdir+'Openness.txt',sep='\t',header=None,index=None) + print('Generate TF binding...') + Element_name_bulk = pd.read_csv(GRNdir+'all_hg19.txt', delimiter="\t", header=None) + Element_name_bulk=Element_name_bulk[0].values + Element_name = RE_pseudobulk.index + motifWeight=pd.read_csv(GRNdir+'motifWeight.txt',index_col=0,sep='\t') + load_TFbinding(GRNdir,motifWeight,Match2,TFName,Element_name,outdir) + print('Generate Index...') + #Read hg19_Peak_hg19_gene_u.txt + merged_s,merged_b=load_corr_RE_TG(List,Element_name,Element_name_bulk,outdir) + from tqdm import tqdm + #Assuming you have imported the necessary libraries and defined the variables + #Create a progress bar for the loop + choosL=List + out=np.empty([len(choosL),4], dtype=object) + for i in tqdm(range(len(choosL))): + choosL_i = choosL[i] + out[i, :] = index_generate(choosL_i,merged_s,merged_b,TFName) + pd.DataFrame(out).to_csv(outdir+'index.txt', sep='\t', header=None, index=None) + elif method=='baseline': + print('Overlap the regions with bulk data ...') + #script_path = os.path.join( "extract_overlap_regions_baseline.sh") + #subprocess.run(["sh", script_path, GRNdir, genome,outdir,workdir]) + extract_overlap_regions(genome,GRNdir,outdir,method) + else: + print('Method:' +method+ 'is not found! Please set method as baseline or LINGER') + +import scanpy as sc +#set some figure parameters for nice display inside jupyternotebooks. +import scipy +import pandas as pd +import anndata +import numpy as np +from scipy.sparse import coo_matrix +from scipy.sparse import csc_matrix +def get_adata(matrix,features,barcodes,label): + ### generate the anndata + matrix.data=matrix.data.astype(np.float32) + adata=anndata.AnnData(X= csc_matrix(matrix.T)) + adata.var['gene_ids']=features[1].values + adata.obs['barcode']=barcodes[0].values + if len(barcodes[0].values[0].split("-"))==2: + adata.obs['sample'] = [int(string.split("-")[1]) for string in barcodes[0].values] + else: + adata.obs['sample'] = 1 + rows_to_select=features[features[2]=='Gene Expression'].index + adata_RNA = adata[:,rows_to_select] + rows_to_select=features[features[2]=='Peaks'].index + adata_ATAC = adata[:,rows_to_select] +### if you have the label (cell type annotation) + idx=adata_RNA.obs['barcode'].isin(label['barcode_use'].values) + adata_RNA=adata_RNA[idx] + adata_ATAC=adata_ATAC[idx] + label.index=label['barcode_use'] + adata_RNA.obs['label']=label.loc[adata_RNA.obs['barcode']]['label'].values + #barcode_indices = np.where(np.isin(adata_RNA.obs['barcode'].values, label['barcode_use'].values))[0] + #adata_ATAC = adata_ATAC[barcode_indices, :] + adata_ATAC.obs['label']=label.loc[adata_ATAC.obs['barcode']]['label'].values + adata_RNA.var["mt"] = adata_RNA.var_names.str.startswith("MT-") + sc.pp.calculate_qc_metrics( + adata_RNA, qc_vars=["mt"], percent_top=None, log1p=False, inplace=True +) + adata_RNA = adata_RNA[adata_RNA.obs.pct_counts_mt < 5, :].copy() + adata_RNA.var.index=adata_RNA.var['gene_ids'].values + adata_RNA.var_names_make_unique() + adata_RNA.var['gene_ids']=adata_RNA.var.index + selected_barcode=list(set(adata_RNA.obs['barcode'].values)&set(adata_ATAC.obs['barcode'].values)) + barcode_idx=pd.DataFrame(range(adata_RNA.shape[0]), index=adata_RNA.obs['barcode'].values) + adata_RNA = adata_RNA[barcode_idx.loc[selected_barcode][0]] + barcode_idx=pd.DataFrame(range(adata_ATAC.shape[0]), index=adata_ATAC.obs['barcode'].values) + adata_ATAC = adata_ATAC[barcode_idx.loc[selected_barcode][0]] + return adata_RNA,adata_ATAC + + +def get_adata_h5(adata_RNA,adata_ATAC,label): + ### generate the anndata + if len(adata_RNA.obs['barcode'].values[0].split("-"))==2: + adata_RNA.obs['sample'] = [int(string.split("-")[1]) for string in adata_RNA.obs['barcode'].values] + adata_ATAC.obs['sample'] = [int(string.split("-")[1]) for string in adata_ATAC.obs['barcode'].values] + else: + adata_RNA.obs['sample'] = 1 + adata_ATAC.obs['sample'] = 1 +### if you have the label (cell type annotation) + idx=adata_RNA.obs['barcode'].isin(label['barcode_use'].values) + adata_RNA=adata_RNA[idx] + adata_ATAC=adata_ATAC[idx] + label.index=label['barcode_use'] + adata_RNA.obs['label']=label.loc[adata_RNA.obs['barcode']]['label'].values + #barcode_indices = np.where(np.isin(adata_RNA.obs['barcode'].values, label['barcode_use'].values))[0] + #adata_ATAC = adata_ATAC[barcode_indices, :] + adata_ATAC.obs['label']=label.loc[adata_ATAC.obs['barcode']]['label'].values + adata_RNA.var["mt"] = adata_RNA.var_names.str.startswith("MT-") + sc.pp.calculate_qc_metrics( + adata_RNA, qc_vars=["mt"], percent_top=None, log1p=False, inplace=True +) + adata_RNA = adata_RNA[adata_RNA.obs.pct_counts_mt < 5, :].copy() + adata_RNA.var.index=adata_RNA.var['gene_ids'].values + adata_RNA.var_names_make_unique() + adata_RNA.var['gene_ids']=adata_RNA.var.index + selected_barcode=list(set(adata_RNA.obs['barcode'].values)&set(adata_ATAC.obs['barcode'].values)) + barcode_idx=pd.DataFrame(range(adata_RNA.shape[0]), index=adata_RNA.obs['barcode'].values) + adata_RNA = adata_RNA[barcode_idx.loc[selected_barcode][0]] + barcode_idx=pd.DataFrame(range(adata_ATAC.shape[0]), index=adata_ATAC.obs['barcode'].values) + adata_ATAC = adata_ATAC[barcode_idx.loc[selected_barcode][0]] + return adata_RNA,adata_ATAC \ No newline at end of file