basic

In [1]:

import pandas as pd
import numpy as np
import scanpy as sc
import scMethtools as scm
import matplotlib.pyplot as plt

import pandas as pd import numpy as np import scanpy as sc import scMethtools as scm import matplotlib.pyplot as plt

 ####   ####  #    # ##### #    # #####  ####   ####  #       ####  
 #      #    # ##  ##   #   #    #   #   #    # #    # #      #      
  ####  #      # ## #   #   ######   #   #    # #    # #       ####  
      # #      #    #   #   #    #   #   #    # #    # #           # 
 #    # #    # #    #   #   #    #   #   #    # #    # #      #    # 
  ####   ####  #    #   #   #    #   #    ####   ####  ######  ####                                                

Version: 1.0.0, Tutorials: https://ngdc.cncb.ac.cn/methbank/scm

In [4]:

import anndata as ad
adata = ad.read_h5ad('/xtdisk/methbank_baoym/zongwt/single/data/GSE56789/promoters_after.h5ad')

import anndata as ad adata = ad.read_h5ad('/xtdisk/methbank_baoym/zongwt/single/data/GSE56789/promoters_after.h5ad')

In [5]:

adata.var

adata.var

Out[5]:

	chromosome	start	end	covered_cell	var	sr_var	mean	pct_cell	feature_select	feature_select_var	Accession	Gene	Distance
index
chr1:3000001-3100000	chr1	3000001	3100000	48	0.053359	0.045232	0.535771	0.923077	False	False	ENSMUSG00000051951	Xkr4	176124
chr1:3100001-3200000	chr1	3100001	3200000	45	0.060264	0.054389	0.517756	0.865385	False	False	ENSMUSG00000051951	Xkr4	76124
chr1:3200001-3300000	chr1	3200001	3300000	46	0.066577	0.061711	0.613043	0.884615	False	False	ENSMUSG00000051951	Xkr4	0
chr1:3300001-3400000	chr1	3300001	3400000	45	0.053558	0.063810	0.638244	0.865385	False	False	ENSMUSG00000051951	Xkr4	0
chr1:3400001-3500000	chr1	3400001	3500000	45	0.061905	0.057469	0.688444	0.865385	False	False	ENSMUSG00000051951	Xkr4	0
...	...	...	...	...	...	...	...	...	...	...	...	...	...
chrY:90400001-90500000	chrY	90400001	90500000	38	0.133379	0.088469	0.541974	0.730769	False	True	ENSMUSG00000096178	Gm20837	0
chrY:90500001-90600000	chrY	90500001	90600000	40	0.087185	0.056909	0.457425	0.769231	False	False	ENSMUSG00000096178	Gm20837	55470
chrY:90600001-90700000	chrY	90600001	90700000	39	0.112958	0.074317	0.538282	0.750000	False	True	ENSMUSG00000096178	Gm20837	155470
chrY:90700001-90800000	chrY	90700001	90800000	49	0.049214	0.039676	0.401429	0.942308	False	False	ENSMUSG00000096178	Gm20837	255470
chrY:90800001-90900000	chrY	90800001	90900000	47	0.028525	0.029491	0.437000	0.903846	False	False	ENSMUSG00000096178	Gm20837	355470

26477 rows × 13 columns

In [7]:

#annotation
gtf_file='/xtdisk/methbank_baoym/zongwt/single/genome/mouse/mm10.ncbiRefSeq.gtf.gz'
genomic_regions = scm.dmr.DmrGenomicRegions(
    gtf_file
)

#annotation gtf_file='/xtdisk/methbank_baoym/zongwt/single/genome/mouse/mm10.ncbiRefSeq.gtf.gz' genomic_regions = scm.dmr.DmrGenomicRegions( gtf_file )

/asnas/baoym_group/zongwt/software/miniconda3/envs/py39/lib/python3.9/site-packages/scMethtools/dmr/annotation.py:42: FutureWarning: The behavior of Series.replace (and DataFrame.replace) with CategoricalDtype is deprecated. In a future version, replace will only be used for cases that preserve the categories. To change the categories, use ser.cat.rename_categories instead.
  df['Feature'] = df['Feature'].replace(feature_mapping)

In [10]:

dmr_list = scm.pp.dmr_df(adata)
dmr_gr = scm.dmr.df_to_pyranges(dmr_list)

dmr_list = scm.pp.dmr_df(adata) dmr_gr = scm.dmr.df_to_pyranges(dmr_list)

In [11]:

#注释时间怎么这么长
#dmr list太大
result = genomic_regions.annotate(dmr_gr)
result

#注释时间怎么这么长 #dmr list太大 result = genomic_regions.annotate(dmr_gr) result

Out[11]:

	group	names	scores	logfoldchanges	pvals	pvals_adj	Chromosome	Start	End	Feature	gene_id	gene_name
0	ESC	chr1:3000001-3100000	5.789519	1.207263	2.341092e-06	4.635439e-06	chr1	3000001	3100000	promoter(1-5kb)	Gm26206	Gm26206
1	ESC	chr1:3100001-3200000	6.902176	1.453446	3.010144e-08	7.540882e-08	chr1	3100001	3200000	promoter(1-5kb)	Gm26206	Gm26206
2	ESC	chr1:3200001-3300000	3.786908	0.779765	5.181727e-04	7.902119e-04	chr1	3200001	3300000	promoter(1-5kb)	Gm18956,Xkr4	Gm18956,Xkr4
3	ESC	chr1:3300001-3400000	0.043596	0.007627	9.654365e-01	9.686926e-01	chr1	3300001	3400000	promoter(1-5kb)	Xkr4	Xkr4
4	ESC	chr1:3400001-3500000	-3.877076	-0.584082	3.568692e-04	5.534692e-04	chr1	3400001	3500000	exon	Xkr4	Xkr4
...	...	...	...	...	...	...	...	...	...	...	...	...
52949	MII oocyte	chrY:90400001-90500000	-0.438300	-0.170020	6.655875e-01	6.923150e-01	chrY	90400001	90500000	promoter(1-5kb)	Gm20837	Gm20837
52950	MII oocyte	chrY:90500001-90600000	-1.786858	-0.763352	9.166696e-02	1.085693e-01	chrY	90500001	90600000	Distal intergenic
52951	MII oocyte	chrY:90600001-90700000	-3.709552	-1.626837	1.969971e-03	2.819706e-03	chrY	90600001	90700000	Distal intergenic
52952	MII oocyte	chrY:90700001-90800000	9.698689	1.155139	2.182844e-12	1.299351e-11	chrY	90700001	90800000	promoter(1-5kb)	Gm52481,Mid1-ps1,Erdr1,Gm33272,G530011O06Rik	Gm52481,Mid1-ps1,Erdr1,Gm33272,G530011O06Rik
52953	MII oocyte	chrY:90800001-90900000	2.675648	0.403259	1.051206e-02	1.396667e-02	chrY	90800001	90900000	promoter(1-5kb)	LOC108168645	LOC108168645

52954 rows × 12 columns

In [19]:

import matplotlib.pyplot as plt
import numpy as np


def plot_peak_annotation_wheel(annotation, col='Feature',out_file=None):
    '''
    Generates a pie chart of the genomic annotations of identified VMR

    Args:
        annnotation           Connection to the .loom file to use
        out_file              Name and location of the output file

    Remarks:

    '''
    simple_annotation = np.array([x.split(' ')[0].lower() for x in annotation[col]])

    counts = {}
    for x in np.unique(simple_annotation):
        counts[x] = np.sum(simple_annotation == x)

    labels = 'promoter(1-5kb)','promoter(<=1kb)','five_prime_utr','three_prime_utr'

    fig, ax = plt.subplots(figsize=(12, 10))
    plt.rcParams["font.size"] = "10"

    size = 0.4
    vals = np.array([counts[x] for x in labels])

    cmap = plt.get_cmap("tab20c")
    outer_colors = cmap(np.arange(2) * 4 + 1)[[1, 0], :]
    inner_colors = cmap(np.array([0, 1, 2, 3, 4]))[::-1, :]

    ax.pie([vals[0], sum(vals[1:])], radius=1, colors=outer_colors,
           wedgeprops=dict(width=size, edgecolor='w'), startangle=90)

    def func(pct, allvals):
        absolute = int(pct / 100. * np.sum(allvals))
        return "{:.1f}%\n({:d})".format(pct, absolute)

    wedges, texts, autotexts = ax.pie(vals.flatten(), radius=1 - size, colors=inner_colors,
                                      autopct=lambda pct: func(pct, vals), pctdistance=1.2,
                                      wedgeprops=dict(width=size, edgecolor='w'), startangle=90)

    ax.set(aspect="equal")
    ax.set_title('Intra- vs. Intergenic Peaks', fontsize=24)

    ax.legend(wedges, labels,
              title="Regions",
              loc="center left",
              bbox_to_anchor=(1, 0, 0.5, 1), frameon=False)

    # plt.savefig(out_file, dpi=300, bbox_inches='tight')

import matplotlib.pyplot as plt import numpy as np def plot_peak_annotation_wheel(annotation, col='Feature',out_file=None): ''' Generates a pie chart of the genomic annotations of identified VMR Args: annnotation Connection to the .loom file to use out_file Name and location of the output file Remarks: ''' simple_annotation = np.array([x.split(' ')[0].lower() for x in annotation[col]]) counts = {} for x in np.unique(simple_annotation): counts[x] = np.sum(simple_annotation == x) labels = 'promoter(1-5kb)','promoter(<=1kb)','five_prime_utr','three_prime_utr' fig, ax = plt.subplots(figsize=(12, 10)) plt.rcParams["font.size"] = "10" size = 0.4 vals = np.array([counts[x] for x in labels]) cmap = plt.get_cmap("tab20c") outer_colors = cmap(np.arange(2) * 4 + 1)[[1, 0], :] inner_colors = cmap(np.array([0, 1, 2, 3, 4]))[::-1, :] ax.pie([vals[0], sum(vals[1:])], radius=1, colors=outer_colors, wedgeprops=dict(width=size, edgecolor='w'), startangle=90) def func(pct, allvals): absolute = int(pct / 100. * np.sum(allvals)) return "{:.1f}%\n({:d})".format(pct, absolute) wedges, texts, autotexts = ax.pie(vals.flatten(), radius=1 - size, colors=inner_colors, autopct=lambda pct: func(pct, vals), pctdistance=1.2, wedgeprops=dict(width=size, edgecolor='w'), startangle=90) ax.set(aspect="equal") ax.set_title('Intra- vs. Intergenic Peaks', fontsize=24) ax.legend(wedges, labels, title="Regions", loc="center left", bbox_to_anchor=(1, 0, 0.5, 1), frameon=False) # plt.savefig(out_file, dpi=300, bbox_inches='tight')

In [21]:

import logging
logger = logging.getLogger()
logging.basicConfig(
    format='%(asctime)s %(levelname)-8s %(message)s',
    level=logging.INFO,
    datefmt='%H:%M:%S')

import logging logger = logging.getLogger() logging.basicConfig( format='%(asctime)s %(levelname)-8s %(message)s', level=logging.INFO, datefmt='%H:%M:%S')

In [25]:

import seaborn as sns
def grouped_boxplot(adata, color_by, value_column, colors=None):
    """
    绘制分组箱线图

    参数:
    adata (AnnData): 包含数据和分组信息的AnnData对象
    color_by (str): 用于分组的列名
    value_column (str): 要绘制箱线图的数值列名
    colors (list): 可选参数，用于指定颜色的列表
    """
    plt.figure(figsize=(10, 6))
    
    if colors is None:
        # 使用Matplotlib的调色盘为不同分组着色，将其转化为Seaborn颜色列表
        n_colors = len(adata.obs[color_by].unique())
        if n_colors > 10:
            colors = [plt.cm.get_cmap('tab20', n_colors)(i) for i in range(n_colors)]
        else:
            colors = [plt.cm.get_cmap('tab10', n_colors)(i) for i in range(n_colors)]
    
    # 绘制箱线图
    ax = sns.boxplot(x=color_by, y=value_column, data=adata.obs, palette=colors)
    
    # 绘制中线
    medians = adata.obs.groupby(color_by)[value_column].median()
    xtick_labels = list(adata.obs[color_by].unique())
    
    for xtick, median in zip(ax.get_xticks(), medians):
        ax.text(xtick, median, f'{median:.2f}', ha='center', va='bottom', color='k', fontsize=10)
    
    plt.xticks(np.arange(len(xtick_labels)), xtick_labels)
    
    plt.title(f'Grouped Boxplot of {value_column} by {color_by}')
    plt.xlabel(color_by)
    plt.ylabel(value_column)
    plt.show()

import seaborn as sns def grouped_boxplot(adata, color_by, value_column, colors=None): """ 绘制分组箱线图 参数: adata (AnnData): 包含数据和分组信息的AnnData对象 color_by (str): 用于分组的列名 value_column (str): 要绘制箱线图的数值列名 colors (list): 可选参数，用于指定颜色的列表 """ plt.figure(figsize=(10, 6)) if colors is None: # 使用Matplotlib的调色盘为不同分组着色，将其转化为Seaborn颜色列表 n_colors = len(adata.obs[color_by].unique()) if n_colors > 10: colors = [plt.cm.get_cmap('tab20', n_colors)(i) for i in range(n_colors)] else: colors = [plt.cm.get_cmap('tab10', n_colors)(i) for i in range(n_colors)] # 绘制箱线图 ax = sns.boxplot(x=color_by, y=value_column, data=adata.obs, palette=colors) # 绘制中线 medians = adata.obs.groupby(color_by)[value_column].median() xtick_labels = list(adata.obs[color_by].unique()) for xtick, median in zip(ax.get_xticks(), medians): ax.text(xtick, median, f'{median:.2f}', ha='center', va='bottom', color='k', fontsize=10) plt.xticks(np.arange(len(xtick_labels)), xtick_labels) plt.title(f'Grouped Boxplot of {value_column} by {color_by}') plt.xlabel(color_by) plt.ylabel(value_column) plt.show()

In [26]:

#画的是value数值型的分类图
grouped_boxplot(adata,color_by='Cell_type',value_column='')

#画的是value数值型的分类图 grouped_boxplot(adata,color_by='Cell_type',value_column='')

/tmp/ipykernel_48110/4082876098.py:20: MatplotlibDeprecationWarning: The get_cmap function was deprecated in Matplotlib 3.7 and will be removed two minor releases later. Use ``matplotlib.colormaps[name]`` or ``matplotlib.colormaps.get_cmap(obj)`` instead.
  colors = [plt.cm.get_cmap('tab10', n_colors)(i) for i in range(n_colors)]
/tmp/ipykernel_48110/4082876098.py:23: FutureWarning: 

Passing `palette` without assigning `hue` is deprecated and will be removed in v0.14.0. Assign the `x` variable to `hue` and set `legend=False` for the same effect.

  ax = sns.boxplot(x=color_by, y=value_column, data=adata.obs, palette=colors)
/tmp/ipykernel_48110/4082876098.py:26: FutureWarning: The default of observed=False is deprecated and will be changed to True in a future version of pandas. Pass observed=False to retain current behavior or observed=True to adopt the future default and silence this warning.
  medians = adata.obs.groupby(color_by)[value_column].median()

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[26], line 1
----> 1 grouped_boxplot(adata,color_by='Cell_type',value_column='Treatment')

Cell In[25], line 26, in grouped_boxplot(adata, color_by, value_column, colors)
     23 ax = sns.boxplot(x=color_by, y=value_column, data=adata.obs, palette=colors)
     25 # 绘制中线
---> 26 medians = adata.obs.groupby(color_by)[value_column].median()
     27 xtick_labels = list(adata.obs[color_by].unique())
     29 for xtick, median in zip(ax.get_xticks(), medians):

File /asnas/baoym_group/zongwt/software/miniconda3/envs/py39/lib/python3.9/site-packages/pandas/core/groupby/groupby.py:2532, in GroupBy.median(self, numeric_only)
   2459 @final
   2460 def median(self, numeric_only: bool = False) -> NDFrameT:
   2461     """
   2462     Compute median of groups, excluding missing values.
   2463 
   (...)
   2530     Freq: MS, dtype: float64
   2531     """
-> 2532     result = self._cython_agg_general(
   2533         "median",
   2534         alt=lambda x: Series(x, copy=False).median(numeric_only=numeric_only),
   2535         numeric_only=numeric_only,
   2536     )
   2537     return result.__finalize__(self.obj, method="groupby")

File /asnas/baoym_group/zongwt/software/miniconda3/envs/py39/lib/python3.9/site-packages/pandas/core/groupby/groupby.py:1998, in GroupBy._cython_agg_general(self, how, alt, numeric_only, min_count, **kwargs)
   1995     result = self._agg_py_fallback(how, values, ndim=data.ndim, alt=alt)
   1996     return result
-> 1998 new_mgr = data.grouped_reduce(array_func)
   1999 res = self._wrap_agged_manager(new_mgr)
   2000 if how in ["idxmin", "idxmax"]:

File /asnas/baoym_group/zongwt/software/miniconda3/envs/py39/lib/python3.9/site-packages/pandas/core/internals/base.py:367, in SingleDataManager.grouped_reduce(self, func)
    365 def grouped_reduce(self, func):
    366     arr = self.array
--> 367     res = func(arr)
    368     index = default_index(len(res))
    370     mgr = type(self).from_array(res, index)

File /asnas/baoym_group/zongwt/software/miniconda3/envs/py39/lib/python3.9/site-packages/pandas/core/groupby/groupby.py:1973, in GroupBy._cython_agg_general.<locals>.array_func(values)
   1971 def array_func(values: ArrayLike) -> ArrayLike:
   1972     try:
-> 1973         result = self._grouper._cython_operation(
   1974             "aggregate",
   1975             values,
   1976             how,
   1977             axis=data.ndim - 1,
   1978             min_count=min_count,
   1979             **kwargs,
   1980         )
   1981     except NotImplementedError:
   1982         # generally if we have numeric_only=False
   1983         # and non-applicable functions
   1984         # try to python agg
   1985         # TODO: shouldn't min_count matter?
   1986         # TODO: avoid special casing SparseArray here
   1987         if how in ["any", "all"] and isinstance(values, SparseArray):

File /asnas/baoym_group/zongwt/software/miniconda3/envs/py39/lib/python3.9/site-packages/pandas/core/groupby/ops.py:830, in BaseGrouper._cython_operation(self, kind, values, how, axis, min_count, **kwargs)
    828 ids, _, _ = self.group_info
    829 ngroups = self.ngroups
--> 830 return cy_op.cython_operation(
    831     values=values,
    832     axis=axis,
    833     min_count=min_count,
    834     comp_ids=ids,
    835     ngroups=ngroups,
    836     **kwargs,
    837 )

File /asnas/baoym_group/zongwt/software/miniconda3/envs/py39/lib/python3.9/site-packages/pandas/core/groupby/ops.py:540, in WrappedCythonOp.cython_operation(self, values, axis, min_count, comp_ids, ngroups, **kwargs)
    536 self._validate_axis(axis, values)
    538 if not isinstance(values, np.ndarray):
    539     # i.e. ExtensionArray
--> 540     return values._groupby_op(
    541         how=self.how,
    542         has_dropped_na=self.has_dropped_na,
    543         min_count=min_count,
    544         ngroups=ngroups,
    545         ids=comp_ids,
    546         **kwargs,
    547     )
    549 return self._cython_op_ndim_compat(
    550     values,
    551     min_count=min_count,
   (...)
    555     **kwargs,
    556 )

File /asnas/baoym_group/zongwt/software/miniconda3/envs/py39/lib/python3.9/site-packages/pandas/core/arrays/categorical.py:2738, in Categorical._groupby_op(self, how, has_dropped_na, min_count, ngroups, ids, **kwargs)
   2736     if kind == "transform":
   2737         raise TypeError(f"{dtype} type does not support {how} operations")
-> 2738     raise TypeError(f"{dtype} dtype does not support aggregation '{how}'")
   2740 result_mask = None
   2741 mask = self.isna()

TypeError: category dtype does not support aggregation 'median'

In [116]:

def stacked_bar(adata, groupby, colorby,orientation='vertical', ax=None,color=None,figsize=(10,6),fontsize=10, show=None,
    save=True):
    # 获取 obs 数据
    obs = adata.obs
    # 检查 uns 字典中是否有颜色信息
    if f'{colorby}_colors' in adata.uns:
        colors = adata.uns[f'{colorby}_colors']
        
    # 创建透视表，用于绘制堆积图
    pivot_table = obs.pivot_table(index=groupby, columns=colorby, aggfunc='size', fill_value=0, observed=False)
    print(pivot_table)
    if color is not None:
        colors = color[:len(pivot_table.columns)]
    else:
        colors = None

    if ax is None:
        fig, ax = plt.subplots(figsize=figsize,dpi=200)
    else:
        fig = ax.get_figure()
        
    # 绘制堆积条形图
    if orientation == 'horizontal':
        pivot_table.plot(kind='barh', stacked=True, color=colors, ax=ax)
        ax.set_xlabel('Counts')
        ax.set_ylabel(groupby)
    else:
        pivot_table.plot(kind='bar', stacked=True, color=colors, ax=ax)
        ax.set_xlabel(groupby)
        ax.set_ylabel('Counts')
    
    # 设置图形标题和标签
    ax.set_title(f'Stacked Bar Plot of {groupby} by {colorby}')
    ax.legend(title=colorby)
    
    # 设置左边和下边的坐标刻度为透明色
    ax.yaxis.tick_left()
    ax.xaxis.tick_bottom()
    ax.xaxis.set_tick_params(color='none')
    ax.yaxis.set_tick_params(color='none')

    plt.show()

def stacked_bar(adata, groupby, colorby,orientation='vertical', ax=None,color=None,figsize=(10,6),fontsize=10, show=None, save=True): # 获取 obs 数据 obs = adata.obs # 检查 uns 字典中是否有颜色信息 if f'{colorby}_colors' in adata.uns: colors = adata.uns[f'{colorby}_colors'] # 创建透视表，用于绘制堆积图 pivot_table = obs.pivot_table(index=groupby, columns=colorby, aggfunc='size', fill_value=0, observed=False) print(pivot_table) if color is not None: colors = color[:len(pivot_table.columns)] else: colors = None if ax is None: fig, ax = plt.subplots(figsize=figsize,dpi=200) else: fig = ax.get_figure() # 绘制堆积条形图 if orientation == 'horizontal': pivot_table.plot(kind='barh', stacked=True, color=colors, ax=ax) ax.set_xlabel('Counts') ax.set_ylabel(groupby) else: pivot_table.plot(kind='bar', stacked=True, color=colors, ax=ax) ax.set_xlabel(groupby) ax.set_ylabel('Counts') # 设置图形标题和标签 ax.set_title(f'Stacked Bar Plot of {groupby} by {colorby}') ax.legend(title=colorby) # 设置左边和下边的坐标刻度为透明色 ax.yaxis.tick_left() ax.xaxis.tick_bottom() ax.xaxis.set_tick_params(color='none') ax.yaxis.set_tick_params(color='none') plt.show()

In [7]:

scm.pl.stacked_plot(adata,groupby='Cell_type',orientation='horizontal',colorby='Treatment',color=scm.pl.zeileis_palette(),save=".png")

scm.pl.stacked_plot(adata,groupby='Cell_type',orientation='horizontal',colorby='Treatment',color=scm.pl.zeileis_palette(),save=".png")

Treatment   2i  serum/LIF  No Treatment
Cell_type                              
ESC         12         21             0
MII oocyte   0          0            12
saving figure to file ./figures/scvelo_stacked.png

In [ ]:

adata.uns

adata.uns

In [43]:

s = scm.pl.red_palette()

s = scm.pl.red_palette()

In [44]:

def venn(sets={}, out='./', palette='bgrc',
             ax=False, ext='png', dpi=300, fontsize=3.5):
    
    from ..utils import venny4py
    venny4py(sets=sets,out=out,ce=palette,asax=ax,ext=ext,
             dpi=dpi,size=fontsize)
    return ax

def venn(sets={}, out='./', palette='bgrc', ax=False, ext='png', dpi=300, fontsize=3.5): from ..utils import venny4py venny4py(sets=sets,out=out,ce=palette,asax=ax,ext=ext, dpi=dpi,size=fontsize) return ax

Out[44]:

['#F0C3C3',
 '#E07370',
 '#CB3E35',
 '#A22E2A',
 '#5A1713',
 '#D3396D',
 '#DBC3DC',
 '#85539B',
 '#5C2B80',
 '#5C4694']

In [ ]: