<jupyter_start><jupyter_text>Pseudobulk Differential Expression Analysis Dependencies:Python:    scanpy    pandasR:    DESeq2Results related to Figure 1C, Nyquist et al 2022Generates the contents of Supplemental Dataset 2Overview:1. Data is saved as raw, pseudobulk gene expression counts by celltype and sample2. DESeq2 is run in R to identify marker genes blocking for donor3. DESeq2 results are read back in and filtered based on percent of cells in each group expressing significant differential genes<jupyter_code>import scanpy as sc
import pandas as pd
import os
import sys
sys.path.insert(0,"../helper_functions")
import pseudobulk_helpers as ph
print("Scanpy version " + sc.__version__)
print("Pandas version " +pd.__version__)

pseudobulk_markers_folder = "../../Results/tables/General_Celltype_Pseudobulk_Marker_Genes/"
# TODO: set this up to download from some cloud storage instead of directly to the objects
adata_no_doublets=sc.read_h5ad("../../Data/processed_data/all_timepoints/adata_no_doublets_FINAL.h5ad")
adata_raw = sc.read_h5ad("../../Data/processed_data/all_timepoints/adata_raw_counts_with_doublets_FINAL.h5ad")
# raw data contains all cells including doublets, so subset to include only non-doublet cells
raw_filtered = adata_raw[adata_no_doublets.obs_names]

# Prepare pseudobulk data

processed_meta = adata_no_doublets.obs[['sample','time_post_partum_days', 'time_post_partum_weeks', 'milk_stage', 'donor',
                      'infant_sick_YN', 'weaning_YN', 'mastisis_YN', 'breast_soreness_YN', 'directly_breastfeeding_YN', 'any_formula_YN',
                      'mother_medications_YN',    'hormonal_birthcontrol_YN', 'daycare_YN', 'vaccines_reported_YN', 'vaccines_list', 'solid_foods_YN', 'reported_infant_medical_events_YN', 'reported_infant_medical_events_description', 
                      'milk stage ordered',"simplified celltypes"]]
pseudo_counts = ph.save_pseudobulk_counts(raw_filtered, processed_meta,celltype_col="simplified celltypes")
pseudo_counts.to_csv("../../Data/processed_data/cell_subset_expression/general_celltype_pseudobulk_counts.csv")

ph.save_metadata_for_pseudobulk_deseq(processed_meta, celltype_col='simplified celltypes',cell_counts=True).to_csv("../../Data/processed_data/cell_subset_expression/general_celltype_pseudobulk_metadata.csv")
# call R script which runs pseudobulk differential expression using DESeq2

os.system("Rscript general_celltype_pseudobulk_differential_expression.R")
# if output is 256, it is an error, to see the error message, run the above command in a terminal window to see error output
# calculate pct expression for filtering 

# this takes a long time to run
# this command calculates the percent of cells expressing so we don't have to, but we are not using the actual result to rank_genes_groups (marker genes)
sc.tl.rank_genes_groups(adata_no_doublets, groupby="simplified celltypes",pts=True) 
cluster_order = ['LC2', 'LC1',  'CSN1S1 macrophages','GPMNB macrophages','neutrophils', 'B cells', 'T cells', 'dendritic cells',   'fibroblasts']

filtered_degenes = {}

for c in cluster_order:
    # read in pseudobulk
    pseudobulk = pd.read_csv(pseudobulk_markers_folder+"cluster_"+c+"_marker_genes_pseudobulk.csv",index_col=0)
    overlap_genes = list(set(pseudobulk.index).intersection(adata_no_doublets.uns["rank_genes_groups"]["pts"].index))
    pseudobulk.loc[overlap_genes, "pts"] = adata_no_doublets.uns["rank_genes_groups"]["pts"].loc[overlap_genes,c]
    pseudobulk.loc[overlap_genes, "pts_rest"] = adata_no_doublets.uns["rank_genes_groups"]["pts_rest"].loc[overlap_genes,c]
    pseudobulk["pts_difference"] = pseudobulk["pts"]-pseudobulk["pts_rest"]
    
    
    # filtering happens here:
    filtered_degenes[c] = pseudobulk[(pseudobulk["padj"]<0.05)&(pseudobulk["pts"] > .3)&(pseudobulk["log2FoldChange"] > .4)]
    
    # save results to file
    pseudobulk.to_csv(pseudobulk_markers_folder+"/marker_genes_with_pcts/"+c.replace(" ","_").replace("/","_")+"_pseudobulk_marker_genes.csv")
    filtered_degenes[c].sort_values("pts_difference",ascending=False).to_csv(pseudobulk_markers_folder+"marker_genes_with_pcts/filtered_"+c.replace(" ","_").replace("/","_")+"_pseudobulk_marker_genes.csv")<jupyter_output><empty_output><jupyter_text>Repeat for epithelial cell subclusters<jupyter_code># Prepare pseudobulk data
adata_epithelial = adata_no_doublets[adata_no_doublets.obs["simplified celltypes"].isin(["LC1","LC2"])]
raw_epithelial = adata_raw[adata_epithelial.obs_names]
processed_meta = adata_epithelial.obs[['sample','time_post_partum_days', 'time_post_partum_weeks', 'milk_stage', 'donor',
                      'infant_sick_YN', 'weaning_YN', 'mastisis_YN', 'breast_soreness_YN', 'directly_breastfeeding_YN', 'any_formula_YN',
                      'mother_medications_YN',    'hormonal_birthcontrol_YN', 'daycare_YN', 'vaccines_reported_YN', 'vaccines_list', 'solid_foods_YN', 'reported_infant_medical_events_YN', 'reported_infant_medical_events_description', 
                      'milk stage ordered',"Epithelial Cell Subclusters"]]
pseudo_counts = ph.save_pseudobulk_counts(raw_epithelial, processed_meta,celltype_col="Epithelial Cell Subclusters")
pseudo_counts.to_csv("../../Data/processed_data/cell_subset_expression/epithelial_celltype_pseudobulk_counts.csv")

ph.save_metadata_for_pseudobulk_deseq(processed_meta, celltype_col="Epithelial Cell Subclusters",cell_counts=True).to_csv("../../Data/processed_data/cell_subset_expression/epithelial_celltype_pseudobulk_metadata.csv")
# call R script which runs pseudobulk differential expression using DESeq2

os.system("R epithelial_celltype_pseudobulk_differential_expression.R")
epi_pseudobulk_markers_folder = "../../Results/tables/Epithelial_Celltype_Pseudobulk_Marker_Genes/"
adata_epithelial.obs["Epithelial Cell Subclusters"].unique()
sc.tl.rank_genes_groups(adata_epithelial, groupby="Epithelial Cell Subclusters",pts=True) # this command calculates the percent of cells expressing so we don't have to

filtered_degenes = {}

for c in adata_epithelial.obs["Epithelial Cell Subclusters"].unique():
    # read in pseudobulk
    pseudobulk = pd.read_csv(epi_pseudobulk_markers_folder+"cluster_"+c+"_marker_genes_pseudobulk.csv",index_col=0)
    overlap_genes = list(set(pseudobulk.index).intersection(adata_epithelial.uns["rank_genes_groups"]["pts"].index))
    pseudobulk.loc[overlap_genes, "pts"] = adata_epithelial.uns["rank_genes_groups"]["pts"].loc[overlap_genes,c]
    pseudobulk.loc[overlap_genes, "pts_rest"] = adata_epithelial.uns["rank_genes_groups"]["pts_rest"].loc[overlap_genes,c]
    pseudobulk["pts_difference"] = pseudobulk["pts"]-pseudobulk["pts_rest"]
    
    
    # filtering happens here:
    filtered_degenes[c] = pseudobulk[(pseudobulk["padj"]<0.05)&(pseudobulk["pts"] > .3)&(pseudobulk["log2FoldChange"] > .4)]
    
    # save results to file
    pseudobulk.to_csv(epi_pseudobulk_markers_folder+"marker_genes_with_pcts/"+c.replace(" ","_").replace("/","_")+"_pseudobulk_marker_genes.csv")
    filtered_degenes[c].sort_values("pts_difference",ascending=False).to_csv(epi_pseudobulk_markers_folder+"marker_genes_with_pcts/filtered_"+c.replace(" ","_").replace("/","_")+"_pseudobulk_marker_genes.csv")
    
  
epi_pseudobulk_time_varying_folder = "../../Results/tables/Epithelial_Celltype_Pseudobulk_Time_Varying_Genes/"
#Filter Secretory Lactocyte time varying genes

des_res_sec = pd.read_csv(epi_pseudobulk_time_varying_folder+"cluster_Secretory Lactocytes_time_varying_genes_pseudobulk.csv",index_col=0)
celltype = "Secretory Lactocytes"
des_res_reduced = des_res_sec.loc[des_res_sec["padj"]<.05]
des_res_reduced = des_res_reduced.loc[des_res_reduced["log2FoldChange"].abs()>.001]
overlap_genes = list(set(des_res_reduced.index).intersection(set(adata_epithelial.uns["rank_genes_groups"]["pts"].index)))

des_res_reduced = des_res_reduced.loc[overlap_genes]
des_res_reduced["pts"] = adata_epithelial.uns["rank_genes_groups"]["pts"].loc[des_res_reduced.index,celltype]
des_res_reduced = des_res_reduced.loc[des_res_reduced["pts"]>.2]
des_res_reduced.to_csv(epi_pseudobulk_time_varying_folder+"secretory_lactocytes_filtered_time_differential_genes.csv")

# Filter LC1 time varying genes

des_res = pd.read_csv(epi_pseudobulk_time_varying_folder+"cluster_LC1_time_varying_genes_pseudobulk.csv",index_col=0)
celltype = "LC1"
res_res_reduced_lum = des_res.loc[des_res["padj"]<.05]
res_res_reduced_lum = res_res_reduced_lum.loc[res_res_reduced_lum["log2FoldChange"].abs()>.001]
overlap_genes = list(set(res_res_reduced_lum.index).intersection(set(adata_epithelial.uns["rank_genes_groups"]["pts"].index)))
res_res_reduced_lum = res_res_reduced_lum.loc[overlap_genes]
res_res_reduced_lum["pts"] = adata_epithelial.uns["rank_genes_groups"]["pts"].loc[res_res_reduced_lum.index,celltype]
res_res_reduced_lum = res_res_reduced_lum.loc[res_res_reduced_lum["pts"]>.2]

res_res_reduced_lum.to_csv(epi_pseudobulk_time_varying_folder+"LC1_filtered_time_differential_genes.csv")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11753586122539753, 'fraction_numerical': 0.006894251084176582, 'mean_word_length': 7.656400384985563, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 17, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '350178', 'license': 'MIT', 'n_tokens_mistral': 3419, 'n_tokens_neox': 3035, 'n_words': 451, 'path': 'Analysis/comparative_analysis/Pseudobulk_Differential_Expression.ipynb', 'repo_name': 'ShalekLab/MIT_Milk_Study', 'size': '8993'}