s6_cnv.py

# reference: https://infercnvpy.readthedocs.io/en/latest/notebooks/tutorial_3k.html

import argparse
# deal with input
parser = argparse.ArgumentParser(description='infer cnv')
parser.add_argument('-i', '--input', help='input file name (h5ad)', type=str, required=True)
args = parser.parse_args()

import os
import sys
import traceback
import warnings
from pathlib import Path

import scanpy as sc
import seaborn as sns
import infercnvpy as cnv
import matplotlib.pyplot as plt

warnings.filterwarnings("ignore")
sc.settings.set_figure_params(figsize=(5, 5))
sc.settings.verbosity = 3

input_file = args.input
GSE_code = input_file.split('_')[1] # d6_GSExxxx_xxx.h5ad
output_file = f'd7_{GSE_code}_cnv.h5ad'
outputdir = Path('cnv')


print(f'''
-------------------------
Input file:
{input_file}

Do:
infer cnv

Output file:
{output_file} (h5ad)
{str(outputdir.absolute())} (dir)
-------------------------
      ''')

def my_chromosome_heatmap(adata, groupby, vmax, **kwargs):
    # Save the original function
    original_heatmap = sc.pl.heatmap

    # Create a new function that calls the original function without a norm parameter
    def new_heatmap(*args, **kwargs):
        kwargs.pop('norm', None)  # Remove the norm parameter if it exists
        return original_heatmap(*args, vmax=vmax, **kwargs)

    # Replace the original function with our new function
    sc.pl.heatmap = new_heatmap

    # Call chromosome_heatmap, which will now use our new function
    result = cnv.pl.chromosome_heatmap(adata, groupby=groupby, **kwargs)

    # Restore the original function
    sc.pl.heatmap = original_heatmap
    return result

if not outputdir.exists():
    outputdir.mkdir(parents=True)

adata = sc.read_h5ad(input_file)
cnv.tl.infercnv(
    adata,
    reference_key="manual_ann",
    reference_cat=[
        "Endothelial",
    ],
    window_size=250,
)

adata.write_h5ad(output_file) # cache
cnv.pl.chromosome_heatmap(adata, groupby="manual_ann", dendrogram=True)
plt.savefig(outputdir / 'f1_celltype_cnv_heatmap.png')

my_chromosome_heatmap(adata, groupby="manual_ann", vmax=1, dendrogram=True)
plt.savefig(outputdir / 'f2_celltype_cnv_heatmap_revised.png')


cnv.tl.pca(adata)
cnv.pp.neighbors(adata)
cnv.tl.leiden(adata)

cnv.pl.chromosome_heatmap(adata, groupby="cnv_leiden", dendrogram=True)
plt.savefig(outputdir / 'f3_leiden_cnv_heatmap.png')

my_chromosome_heatmap(adata, groupby="cnv_leiden", vmax=1, dendrogram=True)
plt.savefig(outputdir / 'f4_leiden_cnv_heatmap_revised.png')
adata.write_h5ad(output_file) # cache

cnv.tl.umap(adata)
cnv.tl.cnv_score(adata)

adata.write_h5ad(output_file) # cache

fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(11, 11))
ax4.axis("off")
cnv.pl.umap(
    adata,
    color="cnv_leiden",
    legend_loc="on data",
    legend_fontoutline=2,
    ax=ax1,
    show=False,
)
cnv.pl.umap(adata, color="cnv_score", ax=ax2, show=False)
cnv.pl.umap(adata, color="manual_ann", ax=ax3)
plt.savefig(outputdir / 'f5_cnv_umap.png')

fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(
    2, 2, figsize=(12, 11), gridspec_kw=dict(wspace=0.5)
)
ax4.axis("off")
sc.pl.umap(adata, color="cnv_leiden", ax=ax1, show=False)
sc.pl.umap(adata, color="cnv_score", ax=ax2, show=False)
sc.pl.umap(adata, color="manual_ann", ax=ax3)
plt.savefig(outputdir / 'f6_cnv_umap2.png')

adata.obs["cnv_status"] = "normal"
adata.obs.loc[
    adata.obs['cnv_leiden'].isin(["13", "17", "18", "19", "22", "14", "25", "26", "11", "20"]), "cnv_status"
] = "tumor"
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5), gridspec_kw=dict(wspace=0.5))
cnv.pl.umap(adata, color="cnv_status", ax=ax1, show=False)
sc.pl.umap(adata, color="cnv_status", ax=ax2)
plt.savefig(outputdir / 'f7_status_assigned.png')

cnv.pl.chromosome_heatmap(adata[adata.obs["cnv_status"] == "tumor", :])
plt.savefig(outputdir / 'f8_tumor_cnv_heatmap.png')

cnv.pl.chromosome_heatmap(adata[adata.obs["cnv_status"] == "normal", :])
plt.savefig(outputdir / 'f9_normal_cnv_heatmap.png')