<jupyter_start><jupyter_text>Table of Contents1&nbsp;&nbsp;miRNA1.1&nbsp;&nbsp;Map files to patients1.1.1&nbsp;&nbsp;Drop unused patients1.2&nbsp;&nbsp;Load all data1.3&nbsp;&nbsp;Scale values1.4&nbsp;&nbsp;Save to files2&nbsp;&nbsp;mRNA2.1&nbsp;&nbsp;Map files to patients2.1.1&nbsp;&nbsp;Drop unused patients2.2&nbsp;&nbsp;Pick gene subset2.2.1&nbsp;&nbsp;Compute variance for all genes2.2.2&nbsp;&nbsp;Define gene subset2.3&nbsp;&nbsp;Load all data2.4&nbsp;&nbsp;Scale values2.5&nbsp;&nbsp;Save to files3&nbsp;&nbsp;DNA methylation3.1&nbsp;&nbsp;Map files to patients3.1.1&nbsp;&nbsp;Drop unused patients3.2&nbsp;&nbsp;Pick probe subset3.2.1&nbsp;&nbsp;Save to file3.2.2&nbsp;&nbsp;Compute variance for all selected probes3.2.3&nbsp;&nbsp;Define probe subset3.3&nbsp;&nbsp;Load all data3.4&nbsp;&nbsp;Replace nan by median value3.5&nbsp;&nbsp;Save to files4&nbsp;&nbsp;CNV4.1&nbsp;&nbsp;Map aliquotes to patients4.2&nbsp;&nbsp;Load all data4.3&nbsp;&nbsp;Are there unchanged genes?4.4&nbsp;&nbsp;Replace aliquote by patient IDs4.4.1&nbsp;&nbsp;Check duplicate patient IDs4.5&nbsp;&nbsp;Encode categorical variables4.5.1&nbsp;&nbsp;Drop unused patients4.6&nbsp;&nbsp;Keep only high-variance genes4.7&nbsp;&nbsp;Save to files  Functions<jupyter_code>def request_file_info(data_type):
    fields = [
        'file_name',
        'cases.submitter_id',
        'cases.samples.sample_type',
        'cases.project.project_id',
        'cases.project.primary_site',
        ]

    fields = ','.join(fields)

    files_endpt = 'https://api.gdc.cancer.gov/files'

    filters = {
        'op': 'and',
        'content':[
            {
            'op': 'in',
            'content':{
                'field': 'files.experimental_strategy',
                'value': [data_type]
                }
            }
        ]
    }

    params = {
        'filters': filters,
        'fields': fields,
        'format': 'TSV',
        'size': '200000'
        }

    response = requests.post(files_endpt, headers = {'Content-Type': 'application/json'}, json = params)

    return pd.read_csv(StringIO(response.content.decode('utf-8')), sep='\t')
def make_patient_file_map(df, base_dir):
    d = {}
    
    for _, row in df.iterrows():
        patient = row['cases.0.submitter_id']
        if patient in d:
            if not isinstance(d[patient], tuple):
                d[patient] = (
                    d[patient],
                    os.path.join(base_dir, row.id, row.file_name))
            else:
                d[patient] += os.path.join(base_dir, row.id, row.file_name),
        else:
            d[patient] = os.path.join(base_dir, row.id, row.file_name)

    return d
def min_max_scale(data, features, groups):
    train = data.loc[data.index.isin(groups['train']), features]

    scaler = MinMaxScaler()
    columns = train.columns
    scaler = scaler.fit(train[columns])
    
    data.loc[data.index.isin(groups['train']), features] = scaler.transform(
        train)
    data.loc[data.index.isin(groups['val']), features] = scaler.transform(
        data.loc[data.index.isin(groups['val']), features])
    data.loc[data.index.isin(groups['test']), features] = scaler.transform(
        data.loc[data.index.isin(groups['test']), features])
    
    return data
def table_to_patient_files(table, dir_path, round_digits=4):
    n = len(table)
    
    i = 0

    for index, row in table.iterrows():
        print('\r' + f'Save data to files: {str(i + 1)}/{n}', end='')
        i+= 1

        target_file = os.path.join(dir_path, str(index) + '.tsv')
        
        with open(target_file, 'w') as f:
            if round_digits is not None:
                f.write('\n'.join(str(round(value, round_digits)) for value in row.values))
            else:
                f.write('\n'.join(str(value) for value in row.values))

    print()
    print()
def merge_all_tables(table_list):
    n = len(table_list)
    
    final_table = pd.DataFrame()

    for i, table in enumerate(table_list):
        print('\r' + f'Merge data tables: {str(i + 1)}/{n}', end='')

        if final_table.empty:
            final_table = table
        else:
            final_table = final_table.join(table)
    
    print()
    print('Final table dimensions:', final_table.shape)
    print()

    return final_table
labels = pd.read_csv('labels.tsv', sep='\t')
labels.head()
id_groups = {
    'train': list(labels.loc[labels['group'] == 'train', 'submitter_id']),
    'val': list(labels.loc[labels['group'] == 'val', 'submitter_id']),
    'test': list(labels.loc[labels['group'] == 'test', 'submitter_id'])}<jupyter_output><empty_output><jupyter_text>miRNA  Map files to patientsUse the [GDC API](https://docs.gdc.cancer.gov/API/Users_Guide/Python_Examples/using-python-to-query-the-gdc-api) to retrieve mapping between file names and patient IDs. Collect results as Pandas `DataFrame`.<jupyter_code>%%time

miRNA_files = request_file_info(data_type='miRNA-Seq')
miRNA_files.shape
miRNA_files.head()
miRNA_files = miRNA_files[miRNA_files['cases.0.project.project_id'].str.startswith('TCGA')]
miRNA_files = miRNA_files[miRNA_files['file_name'].str.endswith('mirbase21.mirnas.quantification.txt')]
miRNA_files = miRNA_files[miRNA_files['cases.0.samples.0.sample_type'] == 'Primary Tumor']
miRNA_files.shape<jupyter_output><empty_output><jupyter_text>Some cases have more than one `*mirbase21.mirnas.quantification.txt` file. I suppose these are replicate runs for the same patient. Accordingly, checking some examples showed that `reads_per_million_miRNA_mapped` values were highly correlated. Here I will simply keep the first one in the table and drop the other ones. This decision may be revisited later on if deemed useful.<jupyter_code>print('All rows:       ', miRNA_files.shape[0])
print('Unique patients:', miRNA_files['cases.0.submitter_id'].unique().shape[0])
miRNA_files = miRNA_files[~miRNA_files.duplicated(subset=['cases.0.submitter_id'], keep='first')]
miRNA_files.shape
file_map = make_patient_file_map(df=miRNA_files, base_dir='/mnt/dataA/TCGA/raw/miRNA-seq/')<jupyter_output><empty_output><jupyter_text>Drop unused patientsKeep only patients present in label data.<jupyter_code>len(file_map)
labels = pd.read_csv('../data/labels.tsv', sep='\t')
len(labels['submitter_id'])
file_map = {k: file_map[k] for k in file_map if k in list(labels['submitter_id'])}
len(file_map)<jupyter_output><empty_output><jupyter_text>Load all data<jupyter_code>def load_all_miRNA_data(patient_file_map):
    n = len(patient_file_map)
    
    dfs = []
    for i, patient in enumerate(patient_file_map):
        print('\r' + f'Load data tables: {str(i + 1)}/{n}', end='')
        df = pd.read_csv(patient_file_map[patient], sep='\t',
                         usecols=['miRNA_ID', 'reads_per_million_miRNA_mapped'],
                         index_col=0)
        df.columns = [patient]
        dfs.append(df)

    print()
    print()

    return dfs
%%time

data_tables = load_all_miRNA_data(patient_file_map=file_map)
%%time

miRNA = merge_all_tables(data_tables)<jupyter_output>Merge data tables: 9616/9616
Final table dimensions: (1881, 9616)

CPU times: user 14min 16s, sys: 2.21 s, total: 14min 18s
Wall time: 14min 16s<jupyter_text>Scale values<jupyter_code>miRNA.shape
miRNA.head()
miRNA = min_max_scale(data=miRNA.T, features=miRNA.T.columns, groups=id_groups)
miRNA.head()
miRNA.loc[miRNA.index.isin(['TCGA-A8-A09K', 'TCGA-EW-A1OX', 'TCGA-A2-A3KD'])]<jupyter_output><empty_output><jupyter_text>Save to files<jupyter_code>%%time

table_to_patient_files(
    miRNA, dir_path=os.path.join(DATA_LOCATION, 'miRNA-seq'), round_digits=4)<jupyter_output><empty_output><jupyter_text>mRNA  Map files to patientsUse the [GDC API](https://docs.gdc.cancer.gov/API/Users_Guide/Python_Examples/using-python-to-query-the-gdc-api) to retrieve mapping between file names and patient IDs. Collect results as Pandas `DataFrame`.<jupyter_code>%%time

mRNA_files = request_file_info(data_type='RNA-Seq')
mRNA_files.shape
mRNA_files.head()<jupyter_output><empty_output><jupyter_text>mRNA_files = mRNA_files[mRNA_files['cases.0.project.project_id'].str.startswith('TCGA')]mRNA_files = mRNA_files[mRNA_files['file_name'].str.endswith('htseq.counts.gz')]mRNA_files = mRNA_files[mRNA_files['cases.0.samples.0.sample_type'] == 'Primary Tumor']mRNA_files.shape<jupyter_code>mRNA_files = mRNA_files[mRNA_files['cases.0.project.project_id'].str.startswith('TCGA')]
mRNA_files = mRNA_files[mRNA_files['file_name'].str.endswith('FPKM-UQ.txt.gz')]
mRNA_files = mRNA_files[mRNA_files['cases.0.samples.0.sample_type'] == 'Primary Tumor']
mRNA_files.shape<jupyter_output><empty_output><jupyter_text>As for miRNA above, when there is more than one file for a single patient just keep the first one. Again, this decision may be revisited later on if deemed useful.<jupyter_code>print('All rows:       ', mRNA_files.shape[0])
print('Unique patients:', mRNA_files['cases.0.submitter_id'].unique().shape[0])
mRNA_files = mRNA_files[~mRNA_files.duplicated(subset=['cases.0.submitter_id'], keep='first')]
mRNA_files.shape
file_map = make_patient_file_map(mRNA_files, base_dir='/mnt/dataA/TCGA/raw/RNA-seq_FPKM-UQ/')<jupyter_output><empty_output><jupyter_text>Drop unused patientsKeep only patients present in label data.<jupyter_code>len(file_map)
labels = pd.read_csv('../data/labels.tsv', sep='\t')
len(labels['submitter_id'])
file_map = {k: file_map[k] for k in file_map if k in list(labels['submitter_id'])}
len(file_map)<jupyter_output><empty_output><jupyter_text>Pick gene subsetMost genes are probably not informative. Check variance of each gene across patients and remove low-variance genes to reduce data size.Calculate variance iteratively for gene subsets, to avoid memory issues.  Compute variance for all genesRun dedicated script (found here in `src/scripts` directory) as follows:```bash$ conda activate ig$ python ./compute_gene_variance.py \   -i /mnt/dataA/TCGA/raw/RNA-seq_FPKM-UQ/   -s 10000 \   -o /mnt/dataA/TCGA/raw/RNA-seq_FPKM-UQ_variance.tsv```  Define gene subsetKeep a sensible number of genes showing the highest variance across patients.Running the model on mRNA data after keeping increasing numbers of genes (500, 750, 1'000, or 2'000) showed slight increases in training data C-index (~2 percentage points). Picked the subset yielding highest validation C-index.<jupyter_code># Load computed variance table
variance_table = pd.read_csv('/mnt/dataA/TCGA/raw/RNA-seq_count_variance.tsv', sep='\t',
                             header=None, index_col=0, names=['Var'])
print(f'# of represented mRNAs: {len(variance_table)}')
variance_table.head()
variance_table.sort_values('Var', ascending=False).head(10)
variance_table.describe()
variance_table.boxplot()
e = 1e-7
np.log10(variance_table + e).boxplot()
len(variance_table[variance_table['Var'] > 1e9])
variance_table[variance_table['Var'] > 1e9].sort_values(by='Var', ascending=False).plot(use_index=False)
# Variance quantile to drop all genes but the top 2k
q = (variance_table.shape[0] - 2e3) / variance_table.shape[0]
print(f'Percentile = {round(q * 100, 1)}%')

variance_table[variance_table['Var'] > variance_table['Var'].quantile(q=q)].shape
# Variance quantile to drop all genes but the top 1k
q = (variance_table.shape[0] - 1e3) / variance_table.shape[0]
print(f'Percentile = {round(q * 100, 1)}%')

variance_table[variance_table['Var'] > variance_table['Var'].quantile(q=q)].shape
# Get IDs for all genes to keep
genes_to_keep = list(variance_table[variance_table['Var'] > variance_table['Var'].quantile(q=q)].index)
genes_to_keep[:5]
len(genes_to_keep)<jupyter_output><empty_output><jupyter_text>Load all data<jupyter_code>def load_all_mRNA_data(patient_file_map, genes):
    n = len(patient_file_map)
    
    dfs = []
    for i, patient in enumerate(patient_file_map):
        print('\r' + f'Load data tables: {str(i + 1)}/{n}', end='')
        df = pd.read_csv(patient_file_map[patient], sep='\t', header=None,
                         index_col=0)
        df = df.loc[genes]
        df.columns = [patient]
        dfs.append(df)

    print()
    print()

    return dfs
%%time

data_tables = load_all_mRNA_data(patient_file_map=file_map, genes=genes_to_keep)
data_tables[0].head()
%%time

mRNA = merge_all_tables(data_tables)<jupyter_output>Merge data tables: 9605/9605
Final table dimensions: (1000, 9605)

CPU times: user 14min 12s, sys: 1.9 s, total: 14min 14s
Wall time: 14min 11s<jupyter_text>Scale values<jupyter_code>mRNA.shape
mRNA.head()
mRNA = min_max_scale(data=mRNA.T, features=mRNA.T.columns, groups=id_groups)
mRNA.head()<jupyter_output><empty_output><jupyter_text>Save to files<jupyter_code>%%time

table_to_patient_files(
    mRNA, dir_path=os.path.join(DATA_LOCATION, 'RNA-seq'), round_digits=4)<jupyter_output><empty_output><jupyter_text>DNA methylationAs explained in the GDC [Methylation Liftover Pipeline page](https://docs.gdc.cancer.gov/Data/Bioinformatics_Pipelines/Methylation_LO_Pipeline/), data was generated using either Illumina Infinium Human Methylation 27 (HM27; 27'578 probes) or HumanMethylation450 (HM450; 485'577 probes).Use the intersection of the probes between the two (25'978 probes).  Map files to patientsUse the [GDC API](https://docs.gdc.cancer.gov/API/Users_Guide/Python_Examples/using-python-to-query-the-gdc-api) to retrieve mapping between file names and patient IDs. Collect results as Pandas `DataFrame`.<jupyter_code>%%time

DNAm_files = request_file_info(data_type='Methylation Array')
DNAm_files.shape
DNAm_files.head()
DNAm_files = DNAm_files[DNAm_files['cases.0.project.project_id'].str.startswith('TCGA')]
DNAm_files = DNAm_files[DNAm_files['file_name'].str.endswith('gdc_hg38.txt')]
DNAm_files = DNAm_files[DNAm_files['cases.0.samples.0.sample_type'] == 'Primary Tumor']
DNAm_files.shape<jupyter_output><empty_output><jupyter_text>As for miRNA above, when there is more than one file for a single patient just keep the first one. Again, this decision may be revisited later, if deemed useful.<jupyter_code>print('All rows:       ', DNAm_files.shape[0])
print('Unique patients:', DNAm_files['cases.0.submitter_id'].unique().shape[0])
DNAm_files = DNAm_files[~DNAm_files.duplicated(subset=['cases.0.submitter_id'], keep='first')]
DNAm_files.shape
file_map = make_patient_file_map(DNAm_files, base_dir='/mnt/dataA/TCGA/raw/Methylation/')<jupyter_output><empty_output><jupyter_text>Drop unused patientsKeep only patients present in label data.<jupyter_code>len(file_map)
labels = pd.read_csv('../data/labels.tsv', sep='\t')
len(labels['submitter_id'])
file_map = {k: file_map[k] for k in file_map if k in list(labels['submitter_id'])}
len(file_map)<jupyter_output><empty_output><jupyter_text>Pick probe subsetCheck a small subset of 100 files to get an idea of the types of data and their proportions in the dataset.<jupyter_code>%%time

# List number of probes per file
total_n_lines = []
i = 0

for patient in file_map:
    i += 1
    if i > 100:
        break
    if isinstance(file_map[patient], tuple):
        total_n_lines = [total_n_lines.append(len(list(pd.read_csv(
            f, sep='\t', index_col=0, usecols=['Composite Element REF', 'Beta_value']).index)))
                         for f in file_map[patient]]
    else:
        total_n_lines.append(len(list(pd.read_csv(
            file_map[patient], sep='\t', index_col=0,
            usecols=['Composite Element REF', 'Beta_value']).index)))
set(total_n_lines)
element = 27578
total_n_lines.count(element)
element = 485577
total_n_lines.count(element)
eg_file = file_map[list(file_map.keys())[100]]
total_n_lines = len(list(pd.read_csv(
    eg_file, sep='\t', index_col=0, usecols=['Composite Element REF', 'Beta_value']).index))
total_n_lines
DNAm = pd.read_csv(eg_file, sep='\t', index_col=0, usecols=['Composite Element REF', 'Beta_value'])
DNAm.head()
DNAm.shape
probes = list(DNAm.index)
probes[:5]
DNAm = pd.read_csv('/mnt/dataA/TCGA/raw/Methylation/0000c40e-9d45-4446-9dd9-a4676224d0ce/jhu-usc.edu_GBM.HumanMethylation450.7.lvl-3.TCGA-19-5955-01A-11D-1697-05.gdc_hg38.txt',
                   sep='\t', index_col=0, usecols=['Composite Element REF', 'Beta_value'])
DNAm.head()
DNAm.shape
DNAm[DNAm.index == 'cg00000292']
%%time

print('# Illumina Infinium Human Methylation 27 (HM27) probes found in HumanMethylation450 (HM450):')
print(sum([p in probes for p in list(DNAm.index)]))
print()<jupyter_output># Illumina Infinium Human Methylation 27 (HM27) probes found in HumanMethylation450 (HM450):
485577
CPU times: user 23min 30s, sys: 19.9 ms, total: 23min 30s
Wall time: 23min 30s<jupyter_text>Use the intersection of the probe sets in the two technologies (25'978 according to the two example files used here).<jupyter_code>probe_set = [p for p in list(DNAm.index) if p in probes]
len(probe_set)
probe_set[:10]<jupyter_output><empty_output><jupyter_text>Save to file<jupyter_code>target_file = '../data/DNAm_probes.tsv'<jupyter_output><empty_output><jupyter_text>with open(target_file, 'w') as f:    f.write('\n'.join(str(value) for value in probe_set))<jupyter_code>pd.read_csv(target_file, sep='\t', header=None, names=['Probes']).head()
pd.read_csv(target_file, sep='\t', header=None, names=['Probes']).shape<jupyter_output><empty_output><jupyter_text>Compute variance for all selected probesRun dedicated script (found here in `src/scripts` directory) as follows:```bash$ conda activate ig$ python ./compute_DNAm_variance.py \   -i /mnt/dataA/TCGA/raw/Methylation/ \   -s 5000 \   -l ../data/labels.tsv \   -p ../data/DNAm_probes.tsv \   -o /mnt/dataA/TCGA/raw/DNA_methylation_Beta-value_variance.tsv```  Define probe subsetKeep a sensible number of probes showing the highest variance across patients.<jupyter_code># Load computed variance table
variance_table = pd.read_csv('/mnt/dataA/TCGA/raw/DNA_methylation_Beta-value_variance.tsv', sep='\t',
                             header=None, index_col=0, names=['Var'])
variance_table.head()
variance_table.sort_values('Var', ascending=False).head(10)
variance_table.describe()
variance_table.boxplot()
e = 1e-7
np.log10(variance_table + e).boxplot()
variance_table[variance_table['Var'] > 1e-2].sort_values(by='Var', ascending=False).plot(use_index=False)
# Variance quantile to drop all genes but the chosen number of probes showing the top variance
# (since there are quite a lot of missing values, the result is lower than the selected number of probes)
q = (variance_table.shape[0] - 4e3) / variance_table.shape[0]
print(f'Percentile = {round(q * 100, 1)}%')

variance_table[variance_table['Var'] > variance_table['Var'].quantile(q=q)].shape
# Get IDs for all probes to keep
probes_to_keep = list(variance_table[variance_table['Var'] > variance_table['Var'].quantile(q=q)].index)
probes_to_keep[:5]
len(probes_to_keep)<jupyter_output><empty_output><jupyter_text>Load all data<jupyter_code>def load_all_DNAm_data(patient_file_map, probes):
    n = len(patient_file_map)
    
    dfs = []
    for i, patient in enumerate(patient_file_map):
        print('\r' + f'Load data tables: {str(i + 1)}/{n}', end='')
        df = pd.read_csv(patient_file_map[patient], sep='\t', index_col=0,
                         usecols=['Composite Element REF', 'Beta_value'])
        df = df.loc[probes]
        df.columns = [patient]
        dfs.append(df)

    print()
    print()

    return dfs
%%time

data_tables = load_all_DNAm_data(patient_file_map=file_map, probes=probes_to_keep)
data_tables[0].head()
%%time

DNAm = merge_all_tables(data_tables)<jupyter_output>Merge data tables: 10257/10257
Final table dimensions: (3601, 10257)

CPU times: user 24min 14s, sys: 3.37 s, total: 24min 17s
Wall time: 24min 43s<jupyter_text>__Note:__ There is no need to scale DNA methylation Beta Values here, since they are already in the interval [0, 1].<jupyter_code>DNAm.shape
DNAm.head()
DNAm=DNAm.T
DNAm.head()<jupyter_output><empty_output><jupyter_text>Replace `nan` by median value<jupyter_code>DNAm.isna().sum()
DNAm.median()
DNAm.fillna(DNAm.median()).isna().sum()
DNAm = DNAm.fillna(DNAm.median())<jupyter_output><empty_output><jupyter_text>Save to files<jupyter_code>%%time

table_to_patient_files(
    DNAm, dir_path=os.path.join(DATA_LOCATION, 'DNAm'), round_digits=4)<jupyter_output><empty_output><jupyter_text>CNV  Map aliquotes to patientsUse the [GDC API](https://docs.gdc.cancer.gov/API/Users_Guide/Python_Examples/using-python-to-query-the-gdc-api) to retrieve mapping of aliquote IDs to patients found in the 33 tables (one per cancer entity). Collect results as Pandas `DataFrame`.<jupyter_code>%%time

response = requests.get('https://api.gdc.cancer.gov/cases',
                        params={
                            'fields': 'cases.samples.portions.analytes.aliquots',
                            'format': 'TSV',
                            'size': '200000'
                        })
cnv_aliquotes = pd.read_csv(StringIO(response.content.decode("utf-8")), sep="\t")
cnv_aliquotes = cnv_aliquotes.set_index('submitter_id')
cnv_aliquotes.shape
cnv_aliquotes.head()
cnv_aliquotes = cnv_aliquotes[cnv_aliquotes.index.str.startswith('TCGA')]
print('Unique patients:', cnv_aliquotes.index.unique().shape[0])
cnv_aliquotes.index.duplicated
aliquote_map = {}

aliquote_cols = [col for col in cnv_aliquotes if col.startswith('aliquot_ids')]

for index, row in cnv_aliquotes.iterrows():
    if index in aliquote_map.keys():
        print(index, 'already in aliquote map!')
    
    aliquote_map[index] = []
    row = row.loc[~row.isnull()]
    cols = [col for col in aliquote_cols if col in row]
    
    for col in cols:
        aliquote_map[index].append(row[col])
len(aliquote_map)<jupyter_output><empty_output><jupyter_text>Load all data<jupyter_code>!ls /mnt/dataA/TCGA/raw/CNV/ | wc -l
def load_all_CNV_data():
    base_dir = '/mnt/dataA/TCGA/raw/CNV/'
    project_dirs = [(os.path.join(base_dir, d)) for d in os.listdir(base_dir)]
    files = [(os.path.join(pd, f)) for pd in project_dirs
             for f in os.listdir(pd) if f.endswith('focal_score_by_genes.txt')]    
    
    dfs = []
    n = len(files)

    for i, f in enumerate(files):
        print('\r' + f'Load data tables: {str(i + 1)}/{n}', end='')
        df = pd.read_csv(f, sep='\t', index_col=0)
        df = df.iloc[:, 2:]
        dfs.append(df)

    print()
    print()

    return dfs
data_tables = load_all_CNV_data()
data_tables[0].head()
cnv = merge_all_tables(data_tables)
cnv.shape<jupyter_output><empty_output><jupyter_text>Are there unchanged genes?Any genes with the same value for all patients will be completely uninformative.<jupyter_code># Any genes with same CNV value across all samples?
unchanged_genes = []

for index, row in cnv.iterrows():
    if len(pd.unique(row)) < 2:
        unchanged_genes.append(index)
print('# of unchanged genes:', len(unchanged_genes))<jupyter_output># of unchanged genes: 0<jupyter_text>Replace aliquote by patient IDs<jupyter_code>def lookup_submitter_id(df, aliquote_id_map):
    cols = df.columns
    n = len(cols)

    for i, col in enumerate(cols):
        print('\r' + f'Replace col names by patient IDs: {str(i + 1)}/{n}', end='')
        for patient in aliquote_id_map:
            if col in aliquote_id_map[patient]:
                patient_id = patient
            
        df.rename(columns={col: patient_id}, inplace=True)
   
    return df
%%time

cnv = lookup_submitter_id(cnv, aliquote_map)
print()
cnv.head()<jupyter_output><empty_output><jupyter_text>Check duplicate patient IDsAre there duplicate column names (more than one aliquote from the same patient)? If so, keep only one after checking that all values are equal.<jupyter_code># Any non-replaced column names
old_names = []
for col in cnv.columns:
    if not col.startswith('TCGA-'):
        old_names.append(col)

print('# columns with old name (not replaced by patient ID):', len(old_names))
# Non-unique column names
len(cnv.columns) - len(set(cnv.columns))
from collections import Counter

c = Counter(cnv.columns)

duplicated_cols = [col for col in c if c[col] > 1]
len(duplicated_cols)
set(c.values())
# Are repeats different?

not_equal_cols = []

for col in duplicated_cols:
    n_repeats = cnv.loc[:, col].shape[1]
    for i in range(n_repeats):
        if not cnv.loc[:, col].iloc[:, 0].equals(cnv.loc[:, col].iloc[:, 1]):
            not_equal_cols.append(col)
len(set(not_equal_cols))
not_equal_cols[:5]
# How different are they? Look at a couple of examples
eg = cnv.loc[:, not_equal_cols[60]]
eg.head()
print('# equal rows:', len(np.where(eg.iloc[:, 0] == eg.iloc[:, 1])[0]))
print('# different rows:', len(np.where(eg.iloc[:, 0] != eg.iloc[:, 1])[0]))<jupyter_output># equal rows: 17985
# different rows: 1744<jupyter_text>Variable number of differences. My non-systematic check returned values from as low as 5 different values to as high as 1'744. For now, just drop all but one column at random.<jupyter_code>cnv = cnv.T
cnv.head()
cnv.shape
# Non-unique column names
cnv[cnv.index.duplicated()].shape
len(cnv)
cnv = cnv.drop_duplicates(keep='first')
len(cnv)
# Non-unique column names
cnv[cnv.index.duplicated()].shape
cnv_copy = cnv.copy()
cnv_copy.shape
idxs = []
to_drop = []
n = cnv_copy.shape[0]

for i, idx in enumerate(cnv_copy.index):
    print('\r' + f'Check patient IDs for repeats: {str(i + 1)}/{n}', end='')
    if idx not in idxs:
        idxs.append(idx)
    else:
        to_drop.append(i)
len(to_drop)
cnv.drop(cnv.index[to_drop], inplace=True)
cnv_copy.shape
cnv.shape<jupyter_output><empty_output><jupyter_text>Encode categorical variablesEncode values (0, 1, or -1) to use with categorical embeddings.<jupyter_code>cnv.head()
label_encoders = {}
for gene in cnv.columns:
    label_encoders[gene] = LabelEncoder()
    label_encoders[gene].fit_transform(cnv.loc[:, gene])

for gene in cnv.columns:
    cnv.loc[:, gene] = label_encoders[gene].fit_transform(cnv.loc[:, gene])
pd.unique(cnv.iloc[:, 0])
pd.unique([int(cnv[col].nunique()) for col in cnv.columns])<jupyter_output><empty_output><jupyter_text>All genes have embedding dimension 3. According to the rule-of-thumb from the Fast.ai course the embedding size should be 2.<jupyter_code>cat_dim = 3
emb_size = min(50, (cat_dim + 1) // 2)
print(f'Embedding size of {emb_size}')
print(f'(categorical dimension {cat_dim} and a total of {cnv.shape[1]} features)')<jupyter_output>Embedding size of 2
(categorical dimension 3 and a total of 19729 features)<jupyter_text>Drop unused patientsKeep only patients present in label data.<jupyter_code>labels = pd.read_csv('../data/labels.tsv', sep='\t')
len(labels['submitter_id'])
cnv.shape
len(cnv.index)
%%time

for patient in cnv.index:
    if patient not in list(labels['submitter_id']):
        cnv.drop(patient, inplace=True)
cnv.shape<jupyter_output><empty_output><jupyter_text>Keep only high-variance genes<jupyter_code>var = cnv.var()
var[var > 1e-2].sort_values(ascending=False).plot(use_index=False)
# Variance quantile to drop all but the chosen number of genes showing the top variance
q = (len(var) - 5e3) / len(var)
print(f'Percentile = {round(q * 100, 1)}%')

var[var > var.quantile(q=q)].shape
# Variance quantile to drop all but the chosen number of genes showing the top variance
q = (len(var) - 2e3) / len(var)
print(f'Percentile = {round(q * 100, 1)}%')

var[var > var.quantile(q=q)].shape
# Get IDs for all genes to keep
genes_to_keep = list(var[var > var.quantile(q=q)].index)
genes_to_keep[:5]
len(genes_to_keep)
cnv = cnv[genes_to_keep]
cnv.shape<jupyter_output><empty_output><jupyter_text>Save to files<jupyter_code>%%time

table_to_patient_files(
    cnv, dir_path=os.path.join(DATA_LOCATION, 'CNV'), round_digits=4)<jupyter_output><empty_output><jupyter_text>Watermark<jupyter_code>%watermark --iversions
%watermark -v
print()
%watermark -u -n<jupyter_output>requests 2.22.0
numpy    1.18.1
pandas   1.0.1

CPython 3.6.7
IPython 7.11.1

last updated: Tue Jul 28 2020
{'fraction_non_alphanumeric': 0.13764367816091955, 'fraction_numerical': 0.020330459770114943, 'mean_word_length': 4.646116406408437, 'pattern_counts': {'":': 0, '<': 150, '<?xml version=': 0, '>': 163, 'https://': 7, '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': '39', 'dir': 'jupyter-structured-clean-dedup', 'id': '406766', 'license': 'MIT', 'n_tokens_mistral': 10606, 'n_tokens_neox': 9903, 'n_words': 2093, 'path': 'data/preprocess_omics.ipynb', 'repo_name': 'luisvalesilva/multisurv', 'size': '27840'}