Goal¶
- implement interval overlap
Tasks¶
- [ ]
Required files¶
-
In [1]:
# Imports
from basepair.imports import *
hv.extension('bokeh')
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
In [2]:
# Common paths
model_dir = Path(f"{ddir}/processed/chipnexus/exp/models/oct-sox-nanog-klf/models/n_dil_layers=9/")
modisco_dir = model_dir / f"modisco/all/profile/"
output_dir = Path("/srv/www/kundaje/avsec/chipnexus/oct-sox-nanog-klf/models/n_dil_layers=9/modisco/all/profile")
In [3]:
# create_tf_session(0)
In [4]:
data = f"{ddir}/processed/chipnexus/exp/models/genomewide/oct-sox-nanog-klf"
In [3]:
ls {data}
In [7]:
from pybedtools import BedTool
In [10]:
btg.head()
In [13]:
btg = BedTool(f"{ddir}/raw/annotation/mm10/mm10.genome.stride1000.w1000.no-blacklist.bed.gz")
dfg = btg.to_dataframe()
dfg['name'] = dfg.index
In [15]:
btg = BedTool.from_dataframe(dfg)
In [18]:
ls {data}
In [19]:
b = BedTool(f"{data}/Sox2-summits.200bp.bed.gz")
In [20]:
b.head()
In [25]:
feature = 'feature1'
In [53]:
import pybedtools
In [ ]:
pybedtools.
In [64]:
from basepair.config import get_data_dir
from basepair.preproc import label_bed
ddir = get_data_dir()
data = f"{ddir}/processed/chipnexus/exp/models/genomewide/oct-sox-nanog-klf"
dfo = label_bed(f"{ddir}/raw/annotation/mm10/mm10.genome.stride1000.w1000.no-blacklist.bed.gz",
{t: f"{data}/{t}-summits.200bp.bed.gz" for t in ['Oct4', 'Sox2', 'Nanog', 'Klf4']})
dfo.to_csv(f"{data}/1kb.osnk.tsv.gz", compression='gzip', sep='\t', index=False)
In [52]:
for t in ['Oct4', 'Sox2', 'Nanog', 'Klf4']:
print(t)
!zcat {data}/{t}-summits.200bp.bed.gz | wc -l
In [49]:
dfo.set_index(['chrom', 'start', 'end'], inplace=True)
dfo.sum(axis=0)
Out[49]:
In [50]:
dfo.mean(axis=0)
Out[50]:
In [28]:
intersected = btg.intersect(b, wa=True, u=True).to_dataframe()['name']
dfg[feature] = 0
dfg.loc[intersected, feature] = 1
In [ ]:
In [23]:
dfi.head()
Out[23]:
In [16]:
dfg.head()
Out[16]:
In [2]:
from gin_train.samplers import StratifiedRandomBatchSampler,iterable_cycle
In [3]:
classes = np.concatenate([np.ones((int(1e7),)), np.zeros((int(1e7),))])
In [4]:
import ipython_memory_usage.ipython_memory_usage as imu
In [5]:
imu.start_watching_memory()
In [25]:
sampler = StratifiedRandomBatchSampler(classes, [0.5, 0.5], 128)
In [38]:
it = iter(sampler)
In [45]:
for i in sampler:
pass
In [36]:
len(sampler)
Out[36]:
In [ ]:
iterable_cycle()
In [5]:
intervals_file = '/users/avsec/workspace/basepair/data/processed/chipnexus/exp/models/genomewide/oct-sox-nanog-klf/1kb.osnk.tsv.gz'
dataspec = '/users/avsec/workspace/basepair/data/processed/chipnexus/exp/models/oct-sox-nanog-klf/dataspec.yml'
In [8]:
from basepair.datasets import get_gw_StrandedProfile_datasets
train, valid = get_gw_StrandedProfile_datasets(dataspec = '/users/avsec/workspace/basepair/data/processed/chipnexus/exp/models/oct-sox-nanog-klf/dataspec.yml',
intervals_file = '/users/avsec/workspace/basepair/data/processed/chipnexus/exp/models/genomewide/oct-sox-nanog-klf/1kb.osnk.tsv.gz',
peak_width=1000,
seq_width=3088,
exclude_chr = ['chrX', 'chrY']
) # use the default train and valid chromosomes
In [9]:
ds = valid[0][1]
In [10]:
valid[0][0]
Out[10]:
In [11]:
it = ds.batch_train_iter(64, num_workers=6)
In [119]:
%tqdm_restart
In [12]:
for i in tqdm(range(100)):
# ds[i]
batch = next(it)
In [13]:
it = ds.batch_train_iter(cycle=False, batch_size=64, num_workers=6)
In [19]:
for i in tqdm(range(100)):
# ds[i]
batch = next(it)
In [96]:
from tqdm import tqdm
In [68]:
%tqdm_restart
In [8]:
from kipoi.data_utils import iterable_cycle
TODO¶
- there is a memory leak in multiprocessing, you should prevent it by storing the arrays as numpy dictionaries
In [5]:
train.tsv.df[0]
Out[5]:
In [1]:
from basepair.data import Dataset
In [2]:
from basepair.datasets import *
In [10]:
import torch
In [63]:
class StrandedProfile(Dataset):
def __init__(self, ds,
peak_width=200,
seq_width=None,
incl_chromosomes=None,
excl_chromosomes=None,
intervals_file=None,
shuffle=True, target_transformer=None):
"""Dataset for loading the bigwigs and fastas
Args:
ds (basepair.src.schemas.DataSpec): data specification containing the
fasta file, bed files and bigWig file paths
chromosomes (list of str): a list of chor
peak_width: resize the bed file to a certain width
bcolz: If True, the bigwig/fasta files are in the genomelake bcolz format
in_memory: If True, load the whole bcolz into memory. Only applicable when bcolz=True
shuffle: True
preprocessor: trained preprocessor object containing the .transform methods
"""
if isinstance(ds, str):
self.ds = DataSpec.load(ds)
else:
self.ds = ds
self.peak_width = peak_width
if seq_width is None:
self.seq_width = peak_width
else:
self.seq_width = seq_width
self.shuffle = shuffle
self.intervals_file = intervals_file
self.incl_chromosomes = incl_chromosomes
self.excl_chromosomes = excl_chromosomes
self.target_transformer = target_transformer
# not specified yet
self.fasta_extractor = None
self.bw_extractors = None
# Load chromosome lengths
fa = FastaFile(self.ds.fasta_file)
self.chrom_lens = {name: l for name, l in zip(fa.references, fa.lengths)}
del fa
self.tsv = TsvReader(self.intervals_file,
num_chr=False,
label_dtype=int,
mask_ambigous=-1,
incl_chromosomes=incl_chromosomes,
excl_chromosomes=excl_chromosomes,
)
self.dfm = self.tsv.df # use the data-frame from tsv
# self.dfmo = {"chrom": self.dfm}
if self.shuffle:
self.dfm = self.dfm.sample(frac=1)
def __len__(self):
return len(self.dfm)
def get_targets(self):
"""
'targets'
"""
assert self.intervals_file is not None
return self.tsv.get_targets()
def __getitem__(self, idx):
if self.fasta_extractor is None:
# first call
self.bw_extractors = {task: [BigwigExtractor(task_spec.pos_counts),
BigwigExtractor(task_spec.neg_counts)]
for task, task_spec in self.ds.task_specs.items()}
self.fasta_extractor = FastaExtractor(self.ds.fasta_file)
# Load the bias model if available
interval = Interval(self.dfm.iat[idx, 0], # chrom
self.dfm.iat[idx, 1], # start
self.dfm.iat[idx, 2]) # end
target_interval = resize_interval(deepcopy(interval), self.peak_width)
seq_interval = resize_interval(deepcopy(interval), self.seq_width)
# task = self.dfm.iat[idx, 3] # task
# TODO - add data augmentation
sequence = self.fasta_extractor([seq_interval])[0]
# cuts = {f"profile/{task}": run_extractors(self.bw_extractors[task],
# [target_interval],
# ignore_strand=spec.ignore_strand)[0]
# for task, spec in self.ds.task_specs.items()}
cuts = {}
task = ''
return {"inputs": sequence,
"targets": cuts,
"metadata": {"range": GenomicRanges(target_interval.chrom,
target_interval.start,
target_interval.stop,
idx),
"interval_from_task": task}}
In [64]:
class Ds(Dataset):
def __init__(self, ds, seq_width=1000, intervals_file=None):
"""Dataset for loading the bigwigs and fastas
Args:
ds (basepair.src.schemas.DataSpec): data specification containing the
fasta file, bed files and bigWig file paths
chromosomes (list of str): a list of chor
peak_width: resize the bed file to a certain width
bcolz: If True, the bigwig/fasta files are in the genomelake bcolz format
in_memory: If True, load the whole bcolz into memory. Only applicable when bcolz=True
shuffle: True
preprocessor: trained preprocessor object containing the .transform methods
"""
if isinstance(ds, str):
self.ds = DataSpec.load(ds)
else:
self.ds = ds
self.tsv = TsvReader(intervals_file,
num_chr=False,
label_dtype=int,
mask_ambigous=-1
)
self.seq_width = seq_width
self.dfm = self.tsv.df # use the data-frame from tsv
self.dfm[0] = self.dfm[0].astype(str)
self.fasta_extractor = None
def __len__(self):
return len(self.dfm)
def get_targets(self):
"""
'targets'
"""
assert self.intervals_file is not None
return self.tsv.get_targets()
def __getitem__(self, idx):
return torch.ones((10000, 10))
# return np.ones((10000, 10))
# if self.fasta_extractor is None:
# self.fasta_extractor = FastaExtractor(self.ds.fasta_file)
# interval = Interval(self.dfm.iat[idx, 0], # chrom
# self.dfm.iat[idx, 1], # start
# self.dfm.iat[idx, 2]) # end
# seq_interval = resize_interval(deepcopy(interval), self.seq_width)
# sequence = self.fasta_extractor([seq_interval])[0]
# # sequence = ''
# return {"input": sequence,
# "s2": sequence,
# "metadata": {"range": GenomicRanges(interval.chrom,
# interval.start,
# interval.stop,
# idx)}}
In [65]:
from torch.utils.data import DataLoader
In [66]:
%tqdm_restart
In [67]:
import torch
In [68]:
intervals_file = '/users/avsec/workspace/basepair/data/processed/chipnexus/exp/models/genomewide/oct-sox-nanog-klf/1kb.osnk.tsv.gz'
dataspec = '/users/avsec/workspace/basepair/data/processed/chipnexus/exp/models/oct-sox-nanog-klf/dataspec.yml'
In [76]:
train = StrandedProfile(dataspec, 10000, intervals_file=intervals_file)
it = train.batch_iter(batch_size=32, num_workers=12)
In [70]:
from kipoi.data_utils import numpy_collate
In [71]:
dl = DataLoader(train, num_workers=12, batch_size=32)
it = iter(dl)
In [75]:
next(it)
Out[75]:
In [28]:
np.array("asd")
Out[28]:
In [73]:
from basepair.data import to_numpy
In [78]:
for i in tqdm(range(10000)):
# if i % 1000 == 0:
# gc.collect()
batch = next(it)
# o = to_numpybatch)
# a= batch['inputs'].numpy()
# del batch
In [94]:
del train
In [97]:
it = train.batch_train_iter(batch_size=32, num_workers=6)
In [90]:
del it
In [98]:
for i in tqdm(range(1000)):
next(it)
In [84]:
import torch
In [ ]: