In [1]:
import pandas as pd
import numpy as np
import os
import sys
import time
import scanpy as sc
import matplotlib.pyplot as plt
1. Load scAND Module¶
In [2]:
sys.path.append('scAND-code/')
import scAND
In [3]:
np.random.seed(2019)
2. Load Data¶
The input of scAND contains three parts:
- A pandas.Series() represents the barcode of cells
- A pandas.Series() represents the region of peaks
- A pandas.DataFrame() contains 3 columns: Peaks, Cells, Counts (The number of Peaks and Cells start as 1). Each row of the DataFrame represents an element in scATAC-seq data.
In [4]:
Count_df = pd.read_csv('Example/Count_df.txt', sep='\t')
cells = pd.read_csv('Example/Count_Cells.txt', sep='\t', header=None)
cells = cells[0]
peaks = pd.read_csv('Example/Count_Peaks.txt', sep='\t', header=None)
peaks = peaks[0]
In [5]:
cells.head()
Out[5]:
In [6]:
peaks.head()
Out[6]:
In [7]:
Count_df.head()
Out[7]:
2.1 Construct Inputs from scATAC-seq Matrix¶
The get_scAND_inputs() function can be used to construct inputs from scATAC-seq matrix. Note that Each row and column of the used matrix should contain at least one element.
In [8]:
MM = pd.read_csv('Example/GSM1647122_GM12878vsHEK.dhsmatrix.txt.gz', sep='\t')
In [9]:
MM.iloc[:3,:10]
Out[9]:
In [10]:
scATAC = MM.iloc[:,4:]
scATAC.index = MM.iloc[:,:4].apply(lambda x: x['chr']+':'+str(x['start'])+'-'+str(x['end']), axis=1)
In [11]:
scATAC.iloc[:3, :5]
Out[11]:
In [12]:
Count_df, cells, peaks = scAND.get_scAND_inputs(scATAC)
In [13]:
cells.head()
Out[13]:
In [14]:
peaks.head()
Out[14]:
In [15]:
Count_df.head()
Out[15]:
3. Run scAND¶
Run_scAND(Count_df, d, weights, cells, peaks, random_seed=2019, L2_norm=True, Binary=True, Graph_norm=True, return_peaks=False, verbose=True)
The parameters of Run_scAND() function:
- Count_df: A pandas.DataFrame() contains 3 columns: Peaks, Cells, Counts (The number of Peaks and Cells start as 1). Each row of the DataFrame represents an element in scATAC-seq data.
- d: The dimensions of the low-dimensional representation of scAND.
- weights: The parameter beta in scAND model. The input should be a list() or a np.array(). scAND can calculate results of different parameters simultaneously.
- cells: A pandas.Series() represents the barcode of cells
- peaks: A pandas.Series() represents the region of peaks
- random_seed: The random seed.
- L2_norm: Logical, should L2 regularization be applied to the scAND representation? Because L2 regularization is related to the dimensions used, we recommend not to do L2 regularization here, but to do L2 regularization in Get_Result() function.
- Binary: Logical, should binarization be applied to the scATAC-seq matrix?
- Graph_norm: Logical, should graph normalization be applied to the adjacency matrix of network?
- return_peaks: Logical, should the scAND representation of peaks be returned?
- verbose: Logical, should the function run silently?
In [16]:
# parameters
beta_list = [0, 0.01, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95]
start_time = time.time()
Rep_cells = scAND.Run_scAND(Count_df=Count_df, d=50, weights=beta_list, cells=cells, peaks=peaks, random_seed=2019,
L2_norm=False, Binary=True, Graph_norm=True, return_peaks=False,verbose=True)
end_time = time.time()
print('Time Costing: %f s' %(end_time-start_time))
In [17]:
used_dim = 10
used_beta = 0.8
used_rep = scAND.Get_Result(Rep_cells, beta=used_beta, dim=used_dim, L2_norm=True)
In [18]:
used_rep.iloc[:5, :5]
Out[18]:
4. T-SNE Visualization or UMAP Visualization Using the scanpy Package¶
In [19]:
metadata = pd.read_csv('Example/metadata.txt', sep='\t', header=None, index_col=0)
metadata.columns = ['label']
metadata.head()
Out[19]:
In [20]:
adata = sc.AnnData(used_rep)
adata.var_names_make_unique()
adata.obs = metadata.copy()
In [21]:
sc.pp.neighbors(adata, random_state=2019)
sc.tl.tsne(adata, random_state=2019)
sc.tl.umap(adata, random_state=2019)
fig, ax = plt.subplots(figsize=(3, 3))
sc.pl.tsne(adata, color='label', title='scAND', s=15, ax=ax)
fig, ax = plt.subplots(figsize=(3, 3))
sc.pl.umap(adata, color='label', title='scAND', s=15, ax=ax)
