<jupyter_start><jupyter_text>StrategyWe will create different [loom](http://linnarssonlab.org/loompy/index.html) file for each ORF ID.The problem however is, we have 2397475 ~ 2.5M ORF ids, so we'll have to have so may files.However, we'll try to minimize the organization burden by having separating everything by transcript ID.Still need 2.5M individual transcript ID files.  Potential issues.1. If the user comes to us with their own query, there's no way for us to not force them to useRiboCop, or at least somehow get this through querying the bam. But what easy way is out there anyway?2. This ORF list is conservative. Ideally we have more than 10M ORFs. But probably not feasible tosupport them all.So we should make peace with these two limitations.<jupyter_code>! rm outfile.loom
transcript_id = 'ENST00000626826.1'
orfs_for_txid = annotation[annotation.transcript_id == transcript_id]

"""
# Stra
#for orf in orfs
loompy.create(filename, matrix, row_attrs, col_attrs)
with loompy.new("outfile.loom") as dsout:
    for sample in samples:
        with loompy.connect(sample) as dsin:
            logging.info(f"Appending {sample}.")
            dsout.add_columns(ds.layers, col_attrs=dsin.col_attrs, row_attrs=dsin.row_attrs)
"""
out_root_dir = '/staging/as/skchoudh/ribopod_hg38_looms/'
for line_number, row in  tqdm(orfs_for_txid.iterrows(), total=orfs_for_txid.shape[0]):
    annotation_orf_id = row['ORF_ID']
    out_dir = os.path.join(out_root_dir, annotation_orf_id)
    mkdir_p(out_dir)
    loom_file_path = os.path.join(out_dir, 'loomfile.loom')
    orf_index = create_index_for_annotation_row(row)
    chrom = row['chrom']
    sample = sample_list[0]      
    filepath = sample[2]
    line = getline(filepath, line_number+2)
    line_splitted = line.split('\t')
    file_orf_id = line_splitted[0]
    assert file_orf_id == annotation_orf_id
    profile = eval(line_splitted[-1])

    row_attrs = ['{}_{}'.format(chrom, pos) for pos in orf_index]

    row_attrs = { "chr_pos": np.array(row_attrs) }
    col_attrs = { "sample":  np.array([sample[0]])}
    profile = np.array([profile]).reshape(len(row_attrs['chr_pos']), 1)
    loompy.create(loom_file_path, profile,
                  row_attrs, col_attrs)

    with loompy.connect(loom_file_path) as dsout:
        for sample in tqdm(sample_list[1:]):
            filepath = sample[2]
            # +2 because the indexing is 1-based and the first 
            # line is a header
            line = getline(filepath, line_number+2)
            """
            ['ORF_ID',
            'ORF_type',
            'status',
            'phase_score',
            'read_count',
            'length',
            'valid_codons',
            'transcript_id',
            'transcript_type',
            'gene_id',
            'gene_name',
            'gene_type',
            'chrom',
            'strand',
            'start_codon',
            'profile']
            """
            line_splitted = line.split('\t')
            file_orf_id = line_splitted[0]
            assert file_orf_id == annotation_orf_id
            profile = eval(line_splitted[-1])
            col_attrs = { "sample":  np.array([sample[0]])}
            profile = np.array([profile]).reshape(len(row_attrs['chr_pos']), 1)

            dsout.add_columns(profile, col_attrs=col_attrs)# row_attrs=row_attrs)
    #break

transcript_id = 'ENST00000626826.1'
orfs_for_txid = annotation[annotation.transcript_id == transcript_id]

"""
# Stra
#for orf in orfs
loompy.create(filename, matrix, row_attrs, col_attrs)
with loompy.new("outfile.loom") as dsout:
    for sample in samples:
        with loompy.connect(sample) as dsin:
            logging.info(f"Appending {sample}.")
            dsout.add_columns(ds.layers, col_attrs=dsin.col_attrs, row_attrs=dsin.row_attrs)
"""
out_root_dir = '/staging/as/skchoudh/ribopod_hg38_looms/'
for sample_index, sample in enumerate(tqdm(sample_list)):
    df = pd.read_table(sample[2], usecols=['ORF_ID', 'profile'])
    for orf_id
    df.apply()
    
    annotation_orf_id = row['ORF_ID']
    out_dir = os.path.join(out_root_dir, annotation_orf_id)
    mkdir_p(out_dir)
    loom_file_path = os.path.join(out_dir, 'loomfile.loom')
    orf_index = create_index_for_annotation_row(row)
    chrom = row['chrom']
    sample = sample_list[0]      
    filepath = sample[2]
    line = getline(filepath, line_number+2)
    line_splitted = line.split('\t')
    file_orf_id = line_splitted[0]
    assert file_orf_id == annotation_orf_id
    profile = eval(line_splitted[-1])

    row_attrs = ['{}_{}'.format(chrom, pos) for pos in orf_index]

    row_attrs = { "chr_pos": np.array(row_attrs) }
    col_attrs = { "sample":  np.array([sample[0]])}
    profile = np.array([profile]).reshape(len(row_attrs['chr_pos']), 1)
    loompy.create(loom_file_path, profile,
                  row_attrs, col_attrs)

    with loompy.connect(loom_file_path) as dsout:
        for sample in tqdm(sample_list[1:]):
            filepath = sample[2]
            # +2 because the indexing is 1-based and the first 
            # line is a header
            line = getline(filepath, line_number+2)
            """
            ['ORF_ID',
            'ORF_type',
            'status',
            'phase_score',
            'read_count',
            'length',
            'valid_codons',
            'transcript_id',
            'transcript_type',
            'gene_id',
            'gene_name',
            'gene_type',
            'chrom',
            'strand',
            'start_codon',
            'profile']
            """
            line_splitted = line.split('\t')
            file_orf_id = line_splitted[0]
            assert file_orf_id == annotation_orf_id
            profile = eval(line_splitted[-1])
            col_attrs = { "sample":  np.array([sample[0]])}
            profile = np.array([profile]).reshape(len(row_attrs['chr_pos']), 1)

            dsout.add_columns(profile, col_attrs=col_attrs)# row_attrs=row_attrs)
    #break

with loompy.connect('/staging/as/skchoudh/ribopod_hg38_looms/ENST00000626826.1_102197843_102197914_72/loomfile.loom') as ds:
    print(list(ds.col_attrs['SRP']))
def batch(iterable, n=50):
    l = len(iterable)
    for ndx in range(0, l, n):
        yield iterable[ndx:min(ndx + n, l)]

def read_batch_tsv(sample_list):
    dfs = [pd.read_table(tsv,
                         usecols=['ORF_ID', 'profile'], 
                         dtype={'ORF_ID':'str', 'profile': 'str'}, 
                         memory_map=True) for srp, srx, tsv in sample_list]
    srps = [srp for srp, srx, tsv in sample_list]
    srxs = [srx for srp, srx, tsv in sample_list]
    col_attrs = {"study": np.array(srps), "experiment": np.array(srxs) }
    return dfs, col_attrs

def write_loom_file(loom_file_path, matrix, col_attrs, row_attrs=None):
    """Create/Update loom file.
    
    Parameters
    ----------
    filepath: string
              Path to write loomfile
    matrix: array_like
            Data matrix
    col_attrs: dict
               A dict of lists with same length as the columns in matrix
    row_attrs: dict
               A dict of lists with same length as the rows in matrix
    
    """
    # Check if loomfile exists
    if os.path.isfile(loom_file_path):
        # Read its columns attributes
        # to see if the new column_attrs are same as previous
        with loompy.connect(loom_file_path) as ds:
            assert sorted(list(col_attrs.keys())) == sorted(list(ds.col_attrs.keys()))
            # Check if nothing is being added again
            for key in col_attrs.keys():
                assert len(set(col_attrs[key]).intersection(set(ds.col_attrs[key]))) == 0
            
        # Add columns
        with loompy.connect(loom_file_path) as dsout:
            dsout.add_columns(matrix, col_attrs=col_attrs)# row_attrs=row_attrs)
    else:
        loompy.create(loom_file_path,
                      matrix,
                      row_attrs=row_attrs,
                      col_attrs=col_attrs)

def get_row_attrs_from_orf(annotation, orf_id):
    orf_row = annotation.loc[annotation.ORF_ID==orf_id].iloc[0]
    orf_index = create_index_for_annotation_row(orf_row)
    chrom = orf_row['chrom']
    row_attrs = ['{}_{}'.format(chrom, pos) for pos in orf_index]
    row_attrs = { "chr_pos": np.array(row_attrs) }
    return row_attrs


def write_loom_for_dfs(loom_file_path, list_of_dfs, col_attrs, row_attrs, orf_id):
    # Given a list of dfs,
    dfs_subset = [df[df.ORF_ID==orf_id].iloc[0] for df in list_of_dfs]
    profile_stacked = [literal_eval(df.profile) for df in dfs_subset]
    matrix = np.array(profile_stacked).T
    write_loom_file(loom_file_path, matrix, col_attrs, row_attrs)

    
    
def write_loom_batches(sample_list, annotation_filepath, batch_size=50):
    print('Reading annotation ... ')
    annotation = pd.read_table(annotation_filepath)
    print('Done!')
    out_root_dir = '/staging/as/skchoudh/ribopod_hg38_looms/'

    ORF_IDS = annotation.ORF_ID.tolist()
    TX_IDS = annotation.transcript_id.tolist()
    for batch_sample in batch(sample_list, 5):
        # Read all the tsvs in this batch
        print('Reading batch tsv ... ')
        list_of_dfs, col_attrs = read_batch_tsv(batch_sample)        
        print('Done! ')
        # For each ORF in all tsvs
        # write a loom file
        # organized by `transcript_id/ORF_ID.loom`
        for tx_id, orf_id in zip(TX_IDS, ORF_IDS):
            row_attrs = get_row_attrs_from_orf(annotation, orf_id)            
            out_dir = os.path.join(out_root_dir, tx_id)
            mkdir_p(out_dir)
            loom_file_path = os.path.join(out_dir, '{}.loom'.format(orf_id))
            write_loom_for_dfs(loom_file_path, list_of_dfs, col_attrs, row_attrs, orf_id)
        del list_of_dfs
        del col_attrs
write_loom_batches(sample_list, annotation_filepath, batch_size=30)
%%timeit
df = pd.read_table('/staging/as/skchoudh/re-ribo-analysis/hg38/SRP010679/ribocop_results_Feb2019_longest/SRX118286_translating_ORFs.tsv',
                            usecols=['ORF_ID', 'profile'])
%%timeit

df = pd.read_table('/staging/as/skchoudh/re-ribo-analysis/hg38/SRP010679/ribocop_results_Feb2019_longest/SRX118286_translating_ORFs.tsv',
                   usecols=['ORF_ID', 'profile'], 
                   dtype={'ORF_ID':'str', 'profile': 'str'}, 
                   memory_map=True)

%timeit profiles = df.profile.apply(lambda x: eval(x))
loom_file_path = '/staging/as/skchoudh/ribopod_hg38_looms/ENST00000487959.5/ENST00000487959.5_147844554_147844688_135.loom'
with loompy.connect(loom_file_path) as ds:
    print(ds.values)
filepath = '/staging/as/skchoudh/re-ribo-analysis/hg38/SRP010679/ribocop_results_Feb2019_longest/SRX118286_translating_ORFs.tsv'
%timeit
def test_getline(filepath, list_of_lines):
    for line_number in tqdm(list_of_lines):
        line = getline(filepath, line_number+2)
        line_splitted = line.split('\t')
        file_orf_id = line_splitted[0]
        profile = eval(line_splitted[-1])
test_getline(filepath, range(df.shape[0]))<jupyter_output>3%|▎         | 67875/2397475 [00:52<26:25, 1469.25it/s]
{'fraction_non_alphanumeric': 0.11378752471687938, 'fraction_numerical': 0.02300916771526155, 'mean_word_length': 2.764208389715832, 'pattern_counts': {'":': 9, '<': 5, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '148999', 'license': 'BSD-2-Clause', 'n_tokens_mistral': 4064, 'n_tokens_neox': 3820, 'n_words': 758, 'path': 'notebooks/loompy-converter.ipynb', 'repo_name': 'saketkc/re-ribo-smk', 'size': '11126'}