import scanpy as sc
import pandas as pd
import anndata as ad
import yaml
import socket
import os
import warnings
outdir = "../data/output"
if not os.path.exists(outdir):
# Create a new directory because it does not exist
os.makedirs(outdir)
with open("../data/resources/iPSC_lines_map.yaml", 'r') as f:
iPSC_lines_map = yaml.load(f, Loader=yaml.FullLoader)["lines"]
indir = "../data"
UpD50Path = "../data/Sample_S20272_157/filtered_feature_bc_matrix/"
DownD50Path = "../data/Sample_S20273_158/filtered_feature_bc_matrix/"
UpD100_1Path = "../data/Sample_S20812_258/filtered_feature_bc_matrix/"
UpD100_2Path = "../data/Sample_S20813_259/filtered_feature_bc_matrix/"
DownD100Path = "../data/Sample_S31807_MET6_GEX/filtered_feature_bc_matrix/"
UpD300Path = "../data/Sample_S33846_C_GEX/filtered_feature_bc_matrix/"
DownD250Path = "../data/Sample_S20814_260/filtered_feature_bc_matrix/"
BranchBCs = sc.read_h5ad(outdir+"/adatas/adataPaga.h5ad")
BranchBCs = BranchBCs.obs_names[BranchBCs.obs["endpoint_GlutamatergicNeurons_both"] == 1].tolist()
Branchtag = "Exc"
pd.Series(["_".join(ds) for ds in [bc.split("_")[1:] for bc in BranchBCs]]).value_counts()
UpD300 3584 UpD50 2258 DownD50 1851 DownD100 699 DownD250 682 UpD100_2 546 UpD100_1 281 dtype: int64
#D50 HVG <------------------------------------------------------------------------------
#D50 HVG <------------------------------------------------------------------------------
#UpD50-----------------------------------------
#UpD50-----------------------------------------
DSname="UpD50"
DSnameDirName="Sample_S20272_157"
adata=sc.read_10x_mtx(indir+'/'+DSnameDirName+"/filtered_feature_bc_matrix/",var_names='gene_symbols', cache=True)
adata.var_names_make_unique()
cellID = pd.read_csv(indir+'/'+DSnameDirName+'/aggregatedCall/aggregatedCall.tsv', sep ="\t",index_col = 0)
adata.obs['cellID'] = cellID.loc[adata.obs_names, 'Consensus']
adata = adata[adata.obs['cellID'].isin(["3391B","KOLF","MIFF1","809","LowQuality","doublet"])]
adata = adata[pd.Series(list(set(adata.obs_names +"_"+DSname).intersection(set(BranchBCs)))).str.split("_", expand=True)[0]]
sc.pp.filter_genes(adata, min_cells=3)
UpD50 = adata.copy()
#DownD50-----------------------------------------
#DownD50-----------------------------------------
DSname="DownD50"
DSnameDirName="Sample_S20273_158"
adata=sc.read_10x_mtx(indir+'/'+DSnameDirName+"/filtered_feature_bc_matrix/",var_names='gene_symbols', cache=True)
adata.var_names_make_unique()
cellID = pd.read_csv(indir+'/'+DSnameDirName+'/aggregatedCall/aggregatedCall.tsv', sep ="\t",index_col = 0)
adata.obs['cellID'] = cellID.loc[adata.obs_names, 'Consensus']
adata = adata[adata.obs['cellID'].isin(["3391B","KOLF","MIFF1","809","LowQuality","doublet"])]
adata = adata[pd.Series(list(set(adata.obs_names +"_"+DSname).intersection(set(BranchBCs)))).str.split("_", expand=True)[0]]
sc.pp.filter_genes(adata, min_cells=3)
DownD50 = adata.copy()
del adata
UpD50.obs['dataset'] = "UpD50"
UpD50.obs_names = [i + "_" + j for i, j in zip(UpD50.obs_names.tolist(), UpD50.obs["dataset"].tolist())]
DownD50.obs['dataset'] = "DownD50"
DownD50.obs_names = [i + "_" + j for i, j in zip(DownD50.obs_names.tolist(), DownD50.obs["dataset"].tolist())]
#Merge datasets
adata = ad.concat([UpD50,DownD50], merge="same")
sc.pp.normalize_total(adata, target_sum=2e4, layers = "all")
sc.pp.log1p(adata)
#scv.pp.filter_genes(adata, min_shared_counts = 30)
#HVG detection
sc.pp.highly_variable_genes(adata, n_top_genes=2000, batch_key = "dataset")
#adata = adata[:, adata.var.highly_variable]
HVGd50 = adata.var[adata.var.highly_variable_nbatches >= 2].index.tolist()
HVGd50Union = adata.var[adata.var.highly_variable_nbatches >= 1].index.tolist()
/usr/local/lib/python3.8/dist-packages/scanpy/preprocessing/_simple.py:251: ImplicitModificationWarning: Trying to modify attribute `.var` of view, initializing view as actual. adata.var['n_cells'] = number /usr/local/lib/python3.8/dist-packages/scanpy/preprocessing/_normalization.py:134: FutureWarning: The `layers` argument is deprecated. Instead, specify individual layers to normalize with `layer`. warn( /usr/local/lib/python3.8/dist-packages/pandas/core/indexing.py:1637: SettingWithCopyWarning: A value is trying to be set on a copy of a slice from a DataFrame See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy self._setitem_single_block(indexer, value, name)
#D100 HVG <--------------------------
#D100 HVG <--------------------------
#D100 HVG <--------------------------
#UpD100_1-----------------------------------------
#UpD100_1-----------------------------------------
DSname="UpD100_1"
DSnameDirName="Sample_S20812_258"
adata=sc.read_10x_mtx(indir+'/'+DSnameDirName+"/filtered_feature_bc_matrix/",var_names='gene_symbols', cache=True)
adata.var_names_make_unique()
cellID = pd.read_csv(indir+'/'+DSnameDirName+'/aggregatedCall/aggregatedCall.tsv', sep ="\t",index_col = 0)
adata.obs['cellID'] = cellID.loc[adata.obs_names, 'Consensus']
adata = adata[adata.obs['cellID'].isin(["3391B","KOLF","MIFF1","809","LowQuality","doublet"])]
adata = adata[pd.Series(list(set(adata.obs_names +"_"+DSname).intersection(set(BranchBCs)))).str.split("_", expand=True)[0]]
sc.pp.filter_genes(adata, min_cells=3)
UpD100_1 = adata.copy()
#UpD100_2-----------------------------------------
#UpD100_2-----------------------------------------
DSname="UpD100_2"
DSnameDirName="Sample_S20813_259"
adata=sc.read_10x_mtx(indir+'/'+DSnameDirName+"/filtered_feature_bc_matrix/",var_names='gene_symbols', cache=True)
adata.var_names_make_unique()
cellID = pd.read_csv(indir+'/'+DSnameDirName+'/aggregatedCall/aggregatedCall.tsv', sep ="\t",index_col = 0)
adata.obs['cellID'] = cellID.loc[adata.obs_names, 'Consensus']
adata = adata[adata.obs['cellID'].isin(["3391B","KOLF","MIFF1","809","LowQuality","doublet"])]
adata = adata[pd.Series(list(set(adata.obs_names +"_"+DSname).intersection(set(BranchBCs)))).str.split("_", expand=True)[0]]
sc.pp.filter_genes(adata, min_cells=3)
UpD100_2 = adata.copy()
#DownD100-----------------------------------------
#DownD100-----------------------------------------
DSname="DownD100"
DSnameDirName="Sample_S31807_MET6"
adata=sc.read_10x_mtx(indir+'/'+DSnameDirName+"/filtered_feature_bc_matrix/",var_names='gene_symbols', cache=True)
adata.var_names_make_unique()
cellID = pd.read_csv(indir+'/'+DSnameDirName+'/aggregatedCall/aggregatedCall.tsv', sep ="\t",index_col = 0)
adata.obs['cellID'] = cellID.loc[adata.obs_names, 'Consensus']
adata = adata[adata.obs['cellID'].isin(["3391B","KOLF","MIFF1","809","LowQuality","doublet"])]
adata = adata[pd.Series(list(set(adata.obs_names +"_"+DSname).intersection(set(BranchBCs)))).str.split("_", expand=True)[0]]
sc.pp.filter_genes(adata, min_cells=3)
DownD100 = adata.copy()
del adata
UpD100_1.obs['dataset'] = "UpD100_1"
UpD100_1.obs_names = [i + "_" + j for i, j in zip(UpD100_1.obs_names.tolist(), UpD100_1.obs["dataset"].tolist())]
UpD100_2.obs['dataset'] = "UpD100_2"
UpD100_2.obs_names = [i + "_" + j for i, j in zip(UpD100_2.obs_names.tolist(), UpD100_2.obs["dataset"].tolist())]
DownD100.obs['dataset'] = "DownD100"
DownD100.obs_names = [i + "_" + j for i, j in zip(DownD100.obs_names.tolist(), DownD100.obs["dataset"].tolist())]
#Merge datasets
adata = ad.concat([UpD100_1,UpD100_2, DownD100], merge="same")
sc.pp.normalize_total(adata, target_sum=2e4, layers = "all")
sc.pp.log1p(adata)
#scv.pp.filter_genes(adata, min_shared_counts = 30)
#HVG detection
sc.pp.highly_variable_genes(adata, n_top_genes=2000, batch_key = "dataset")
#adata = adata[:, adata.var.highly_variable]
HVGd100 = adata.var[adata.var.highly_variable_nbatches >= 2].index.tolist()
HVGd100Union = adata.var[adata.var.highly_variable_nbatches >= 1].index.tolist()
/usr/local/lib/python3.8/dist-packages/scanpy/preprocessing/_simple.py:251: ImplicitModificationWarning: Trying to modify attribute `.var` of view, initializing view as actual. adata.var['n_cells'] = number /usr/local/lib/python3.8/dist-packages/scanpy/preprocessing/_normalization.py:134: FutureWarning: The `layers` argument is deprecated. Instead, specify individual layers to normalize with `layer`. warn( /usr/local/lib/python3.8/dist-packages/pandas/core/indexing.py:1637: SettingWithCopyWarning: A value is trying to be set on a copy of a slice from a DataFrame See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy self._setitem_single_block(indexer, value, name)
#D250 HVG <--------------------------
#D250 HVG <--outdir----------------------
#D250 HVG <--------------------------
#UpD300-----------------------------------------
#UpD300-----------------------------------------
DSname="UpD300"
DSnameDirName="Sample_S33846_C_GEX"
adata=sc.read_10x_mtx(indir+'/'+DSnameDirName+"/filtered_feature_bc_matrix/",var_names='gene_symbols', cache=True)
adata.var_names_make_unique()
cellID = pd.read_csv(indir+'/'+DSnameDirName+'/aggregatedCall/aggregatedCall.tsv', sep ="\t",index_col = 0)
adata.obs['cellID'] = cellID.loc[adata.obs_names, 'Consensus']
adata = adata[adata.obs['cellID'].isin(["3391B","KOLF","MIFF1","809","LowQuality","doublet"])]
adata = adata[pd.Series(list(set(adata.obs_names +"_"+DSname).intersection(set(BranchBCs)))).str.split("_", expand=True)[0]]
sc.pp.filter_genes(adata, min_cells=3)
UpD300 = adata.copy()
#DownD250-----------------------------------------
#DownD250-----------------------------------------
DSname="DownD250"
DSnameDirName="Sample_S20814_260"
adata=sc.read_10x_mtx(indir+'/'+DSnameDirName+"/filtered_feature_bc_matrix/",var_names='gene_symbols', cache=True)
adata.var_names_make_unique()
cellID = pd.read_csv(indir+'/'+DSnameDirName+'/aggregatedCall/aggregatedCall.tsv', sep ="\t",index_col = 0)
adata.obs['cellID'] = cellID.loc[adata.obs_names, 'Consensus']
adata = adata[adata.obs['cellID'].isin(["3391B","KOLF","MIFF1","809","LowQuality","doublet"])]
adata = adata[pd.Series(list(set(adata.obs_names +"_"+DSname).intersection(set(BranchBCs)))).str.split("_", expand=True)[0]]
sc.pp.filter_genes(adata, min_cells=3)
DownD250 = adata.copy()
del adata
UpD300.obs['dataset'] = "UpD300"
UpD300.obs_names = [i + "_" + j for i, j in zip(UpD300.obs_names.tolist(), UpD300.obs["dataset"].tolist())]
DownD250.obs['dataset'] = "DownD250"
DownD250.obs_names = [i + "_" + j for i, j in zip(DownD250.obs_names.tolist(), DownD250.obs["dataset"].tolist())]
#Merge datasets
adata = ad.concat([UpD300,DownD250], merge="same")
sc.pp.normalize_total(adata, target_sum=2e4, layers = "all")
sc.pp.log1p(adata)
#scv.pp.filter_genes(adata, min_shared_counts = 30)
#HVG detection
sc.pp.highly_variable_genes(adata, n_top_genes=2000, batch_key = "dataset")
#adata = adata[:, adata.var.highly_variable]
HVGd250 = adata.var[adata.var.highly_variable_nbatches >= 2].index.tolist()
HVGd250Union = adata.var[adata.var.highly_variable_nbatches >= 1].index.tolist()
HVGoverlap = list(set(HVGd100+HVGd50+HVGd250))
HVUnionsGoverlap = list(set(HVGd250Union+HVGd50Union+HVGd100Union))
/usr/local/lib/python3.8/dist-packages/scanpy/preprocessing/_simple.py:251: ImplicitModificationWarning: Trying to modify attribute `.var` of view, initializing view as actual. adata.var['n_cells'] = number /usr/local/lib/python3.8/dist-packages/scanpy/preprocessing/_normalization.py:134: FutureWarning: The `layers` argument is deprecated. Instead, specify individual layers to normalize with `layer`. warn( /usr/local/lib/python3.8/dist-packages/pandas/core/indexing.py:1637: SettingWithCopyWarning: A value is trying to be set on a copy of a slice from a DataFrame See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy self._setitem_single_block(indexer, value, name)
External = []
HVGoverlap_curated = list(set(HVGoverlap+External))
pd.DataFrame(HVGoverlap_curated, columns = ["HVG"]).to_csv(outdir+"/HVG_list_intersection_Curated_"+Branchtag+"_PagaTest.txt", sep="\t")
len(HVGoverlap_curated)
1856
len(HVGoverlap_curated)
1856