<jupyter_start><jupyter_text>B) Concatenated omics logistic LASSO<jupyter_code>CV_REPEATS_DEMONSTRATIONS = 20
CV_REPEATS_FINAL_MODEL = 100
concatenated = pd.concat([
    raw_rna_matrix.T.add_suffix('.RNA').T,
    raw_protein_matrix.T.add_suffix('.Protein').T
], axis=0)
concatenated.T.head().T<jupyter_output><empty_output><jupyter_text>1. Using the intersection - shared normalization of blocks Subsetting is required before normalization in this variant as otherwise we have a lot of NaNs:<jupyter_code>concat_intersect = concatenated.loc[:, intersect_cm_and_tmdr]
conditions_intersect = concat_intersect.columns.map(clinical.Meningitis)
conditions_intersect.value_counts()
concatenated_or_zero = concatenated.fillna(0)
%%R -i concatenated_or_zero -i conditions_intersect
colnames(concatenated_or_zero) = remove_leading_X(colnames(concatenated_or_zero))
OMICS_MATRICES[['concatenated']] = concatenated_or_zero
CONDITION_VECTORS[['concatenated']] = conditions_intersect
lasso_pipeline_concat_intersect_shared = create_lasso_pipeline(
    'concatenated',
    smallest_group_n=7,
    subset=intersect_cm_and_tmdr,
    outliers=RNA_OUTLIERS,
    verbose=F,
    omics={'concatenated': concatenated_or_zero}
)<jupyter_output>[1] "Retaining: 55.77%"<jupyter_text>For increased performance, copy to the R environment once only:<jupyter_code>concat_intersect_lasso_shared = cross_validate(
    lasso_pipeline_concat_intersect_shared,
    add_supervision_block(
        {'x': subset(concatenated_or_zero, intersect_cm_and_tmdr, outliers=RNA_OUTLIERS)},
        intersect_cm_and_tmdr
    ),
    test_data=add_supervision_block(
        {'x': subset(concatenated_or_zero, cm_and_tmdr_validation, outliers=RNA_OUTLIERS)},
        cm_and_tmdr_validation
    ),
    n=CV_REPEATS_DEMONSTRATIONS,
    verbose=F,
    stratify=T,
    early_normalization=False
)
roc_auc = concat_intersect_lasso_shared.cross_validation_results.roc_auc
%%R -i roc_auc
(
    plot_roc_auc(roc_auc, mean_auc=T)
    + labs(
        title='Average train set cross-validation ROC: sparse logistic regression',
        subtitle='Shared normalization of RNA and Protein blocks'
    )
)
concat_intersect_lasso_shared.cross_validation.contributions.x.add_weighted_auc(
    concat_intersect_lasso_shared.sub_sampling_test_results.cv_auc
)
concat_intersect_contribs_shared = concat_intersect_lasso_shared.cross_validation.contributions.x.data
%%R -w 700 -i concat_intersect_contribs_shared

coefficients_volcano_plot(
    concat_intersect_contribs_shared,
    n=15, p_value='weighted_auc', neg=F, point_color='grey40',
    highlight_frequent=0.9
) + ylab('Mean contribution value')<jupyter_output><empty_output><jupyter_text>2. Using the intersection - separate normalizations<jupyter_code>from machine_learning.multi_block_pipeline import MultiBlockPipeline
from machine_learning.data_classes import BlocksWrapper
from machine_learning.combine import BlocksConcatenation
from machine_learning.adapter import SklearnAdapter
from machine_learning.utilities import Suffixer
def create_block_lasso_pipeline(smallest_group_n, outliers, subset=None, verbose=False, **kwargs):
    return MultiBlockPipeline(
        block_pipelines={
            'Protein': make_pipeline(
                preprocess('protein', smallest_group_n, outliers, verbose=verbose, **kwargs),
                Suffixer('.Protein')
            ),
            'RNA': make_pipeline(
                preprocess('rna', smallest_group_n, outliers, verbose=verbose, **kwargs),
                Suffixer('.RNA')
            ),
            'Meningitis': make_pipeline(
                OutliersFilter(outlier_patients=outliers, verbose=verbose),
                ConditionsVectorizer(mapping=clinical.Meningitis, case_class='Tuberculosis')
            )
        },
        combine=BlocksConcatenation(blocks=['Protein', 'RNA']),
        adapter=SklearnAdapter(Combined='X', Meningitis='y'),
        model=LogisticRegressionCV(
            penalty='l1', solver='liblinear',
            class_weight='balanced',
            max_iter=150
        ),
        predict=predict_proba
    )
lasso_pipeline_concat_intersect = create_block_lasso_pipeline(7, outliers=RNA_OUTLIERS)
train_data = (
    BlocksWrapper({
        'Protein': subset(raw_protein_matrix, intersect_cm_and_tmdr),
        'RNA': subset(raw_rna_matrix, intersect_cm_and_tmdr, outliers=RNA_OUTLIERS)
    })
    .align_indices()
    .add_supervision_block(intersect_cm_and_tmdr, name='Meningitis')
    .blocks
)
train_data['Protein'].index.map(clinical.Meningitis).value_counts()
train_data['RNA'].index.map(clinical.Meningitis).value_counts()<jupyter_output><empty_output><jupyter_text>The test data is only available for patients with single omics, not for patients with both omics, as this was the criteria for initial train-test split:<jupyter_code>assert not set(cm_and_tmdr_validation).intersection(set(patients_with_both))<jupyter_output><empty_output><jupyter_text>If I want to have a test set using these data, it would be degenerate with half of the data (for one or the other omic) missing. The tests could be performed twice, once with RNA missing, once with protein missing as such it would be the best to perform the tests outside of this function which allows for one test only:<jupyter_code>concat_intersect_lasso = cross_validate(
    lasso_pipeline_concat_intersect,
    train_data,
    test_data=None,
    n=CV_REPEATS_FINAL_MODEL,
    coefficients={'combined': 'coef_'},
    verbose=T,
    stratify=T,
    early_normalization=False
)<jupyter_output>Fitting cross-validation models...<jupyter_text>Test case 1: Use patients with Protein data only (and set missing RNA observations to zero):<jupyter_code>test_protein = (
    BlocksWrapper({
        'Protein': subset(raw_protein_matrix, cm_and_tmdr_validation, outliers=RNA_OUTLIERS),
        'RNA': subset(raw_rna_matrix, cm_and_tmdr_validation, outliers=RNA_OUTLIERS)
    })
    .nullify_exclusively('RNA')
    .transform_single_block('RNA', Suffixer('.RNA'), mark=True)
    .add_supervision_block(cm_and_tmdr_validation, name='Meningitis')
)

test_protein_dataset = test_protein.to_dataset(
    case_class='Tuberculosis',
    response=clinical.Meningitis
)
protein_test_result = concat_intersect_lasso.validate(test_protein_dataset)<jupyter_output><empty_output><jupyter_text>Test case 2: Use patients with RNA data only (and set missing Protein observations to zero):<jupyter_code>test_rna = (
    BlocksWrapper({
        'Protein': subset(raw_protein_matrix, cm_and_tmdr_validation, outliers=RNA_OUTLIERS),
        'RNA': subset(raw_rna_matrix, cm_and_tmdr_validation, outliers=RNA_OUTLIERS)
    })
    .nullify_exclusively('Protein')
    .transform_single_block('Protein', Suffixer('.Protein'), mark=True)
    .add_supervision_block(cm_and_tmdr_validation, name='Meningitis')
)

test_rna_dataset = test_rna.to_dataset(
    case_class='Tuberculosis',
    response=clinical.Meningitis
)
rna_test_result = concat_intersect_lasso.validate(test_rna_dataset)
rna_test_result_sub = concat_intersect_lasso.cross_validation.validate(test_rna_dataset)
protein_test_result_sub = concat_intersect_lasso.cross_validation.validate(test_protein_dataset)
roc_auc_concat = pd.concat([
    concat_intersect_lasso.cross_validation_results.roc_auc.assign(group='CV train: TMD+TMR vs CM', models='Cross-Validation', set='Train: TMD+TMR vs CM'),
    protein_test_result_sub.roc_auc.assign(group='Protein CV test: TMD vs CM', models='Cross-Validation', set='Train: Protein, TMD+TMR vs CM'),
    rna_test_result_sub.roc_auc.assign(group='RNA CV test: TMD vs CM', models='Cross-Validation', set='Train: RNA, TMD+TMR vs CM'),
    protein_test_result.roc_auc.assign(group='Protein full test: TMD vs CM', models='Full model', set='Test: Protein, TMD+TMR vs CM'),
    rna_test_result.roc_auc.assign(group='RNA full test: TMD vs CM', models='Full model', set='Test: RNA, TMD+TMR vs CM')
])
%%R -i roc_auc_concat -w 550 -h 550

colors = c(
    'Train: TMD+TMR vs CM'='#999999',
    'Test: Protein, TMD+TMR vs CM'='red',
    'Test: RNA, TMD+TMR vs CM'='purple',
    'Train: Protein, TMD+TMR vs CM'='red',
    'Train: RNA, TMD+TMR vs CM'='purple'
)

roc_auc_tmdr_cm(roc_auc_concat, annotation=0.15, colors=colors)
len(test_rna.blocks['RNA']), len(test_protein.blocks['Protein'])<jupyter_output><empty_output><jupyter_text>Features driving models with high AUC 100 models were trained with two omics at once but tested on data from patients having a single omic only (either RNA or Protein). The performance of each model in both scenarios (testing with RNA or with Protein only) was averaged:<jupyter_code>average_auc = (Series(protein_test_result_sub.cv_auc) + Series(rna_test_result_sub.cv_auc)) / 2
average_auc.head(3)<jupyter_output><empty_output><jupyter_text>and the features included in the successful models (having high AUC) were visualized, weighting by the contribution strength:<jupyter_code>concat_intersect_lasso.cross_validation.contributions.combined.add_weighted_auc(average_auc)
concat_intersect_lasso_contribs = concat_intersect_lasso.cross_validation.contributions.combined.data
%%R -w 700 -i concat_intersect_lasso_contribs
coefficients_volcano_plot(
    concat_intersect_lasso_contribs,
    n=15, p_value='weighted_auc', neg=F, point_color='grey40',
    highlight_frequent=0.75
) + ylab('Mean contribution value')<jupyter_output><empty_output><jupyter_text>On the bottom there are gene products frequently selected by the "bad" models. On the top there are gene products frequently selected in the "good" models and the gene products being the best contributors of the "good models" even if selected rarely.This visualization is difficult to interpret given that the validation was performed on one omic at a time, thus we cannot read to much out of that. There might be interactions between VAPB, NCBP3 and some proteins which are only improving the classification results (and therefore AUC) if both omics are provided.On the other hand, in the near future the clinical biomarker panels are more likely to be composed of one omic only, thus the positive results are still interesting.The top results - proteins: MMP10, TNFSF15, CDH1 and transcript of GPNMB are consistent with previous analyses.AKT1 and MBD3L2 seem to be novel findings.PSMC1P3 - HIV1 related?The gene products which contribute highly to good models but are rarely selected are likely not significant, and while some look nice (IL17RA) we should not cherry-pick. Instead those could be tested for significance with a permutation test.  Using internal CV? One idea would be to look at the AUC of internal CV. This may highlight some Protein-RNA selection preferences:<jupyter_code>concat_intersect_lasso.cross_validation.contributions.combined.add_weighted_auc(
    concat_intersect_lasso.cross_validation_results.cv_auc
)
concat_intersect_lasso_cv = concat_intersect_lasso.cross_validation.contributions.combined.data
%%R -w 700 -i concat_intersect_lasso_cv
coefficients_volcano_plot(
    concat_intersect_lasso_cv,
    n=15, p_value='weighted_auc', neg=F, point_color='grey40',
    highlight_frequent=0.75
) + ylab('Mean contribution value')<jupyter_output><empty_output><jupyter_text>TODO: - ROC/AUC - forward-regression for the top 4. - PLOS computional biology (look up the recent papers) < 8k definite This is kind of reassuring to see that the gene products prioritetized by the average model performance on single omics are comparable to those performing well in the internal CV.<jupyter_code>auc_comparison = DataFrame(dict(
    average_of_single_omics=average_auc,
    internal_cv=concat_intersect_lasso.cross_validation_results.cv_auc
))
auc_comparison.plot.scatter(x='average_of_single_omics', y='internal_cv');
from scipy.stats import spearmanr
spearmanr(auc_comparison.average_of_single_omics, auc_comparison.internal_cv)<jupyter_output><empty_output><jupyter_text>Maybe we should correct for the internal CV performance? or maybe not..  3. With a new test-set set aside ("small model") So this was fun, but I got nothing to evaluate on except for half empty test dataset (with either RNA of Protein data missing). Instead, I can set some CM aside and use TMS patients with data for both omics; this is not exactly what was done for O2-PLS, however is still interesting to see if we get the same trend (TMS-CM performing better).<jupyter_code>train_data['Protein'].index.map(clinical.Meningitis).value_counts()<jupyter_output><empty_output><jupyter_text>First off, no need to worry about outliers:<jupyter_code>patients_with_both[patients_with_both.isin(RNA_OUTLIERS)]<jupyter_output><empty_output><jupyter_text>I take all 4 TMS patients for whom there are both omics, which gives 11:4 balance for TBM:<jupyter_code>tms_with_both = patients_with_both[patients_with_both.str.endswith('TMS')]
len(tms_with_both)<jupyter_output><empty_output><jupyter_text>I take random 5 CM from what previously was the test set; this way, I have 10:4 balance for CM:<jupyter_code>cm_with_both_test = patients_with_both[patients_with_both.str.endswith('CM')].sample(4, random_state=0)
len(cm_with_both_test)
test_for_small_model = pd.concat([tms_with_both, cm_with_both_test])
test_for_small_model.map(clinical.Meningitis).value_counts()<jupyter_output><empty_output><jupyter_text>And I construct train set without those excluded for the test set:<jupyter_code>train_patients_small = Series(list(
    set(patients_with_both[patients_with_both.str.contains('TM|CM')])
    .difference(test_for_small_model)
))
train_patients_small.map(clinical.Meningitis).value_counts()
small_lasso_concat = cross_validate(
    lasso_pipeline_concat_intersect,
    train_data=(
        BlocksWrapper({
            'Protein': subset(raw_protein_matrix, train_patients_small, outliers=RNA_OUTLIERS),
            'RNA': subset(raw_rna_matrix, train_patients_small, outliers=RNA_OUTLIERS)
        })
        .align_indices()
        .add_supervision_block(train_patients_small, name='Meningitis')
        .blocks
    ),
    test_data=(
        BlocksWrapper({
            'Protein': subset(raw_protein_matrix, test_for_small_model, outliers=RNA_OUTLIERS),
            'RNA': subset(raw_rna_matrix, test_for_small_model, outliers=RNA_OUTLIERS)
        })
        .align_indices()
        .add_supervision_block(test_for_small_model, name='Meningitis')
        .blocks
    ),
    n=CV_REPEATS_DEMONSTRATIONS,
    coefficients={'combined': 'coef_'},
    verbose=False,
    stratify=T,
    early_normalization=False
)
roc_auc_concat_small = pd.concat([
    small_lasso_concat.cross_validation_results.roc_auc.assign(group='CV train: TMD+TMR vs CM', models='Cross-Validation', set='Train: TMD+TMR vs CM'),
    small_lasso_concat.sub_sampling_test_results.roc_auc.assign(group='CV test: TMS vs CM', models='Cross-Validation', set='Test: TMS vs CM'),
    small_lasso_concat.test_result.roc_auc.assign(group='Full test: TMS vs CM', models='Full model', set='Test: TMS vs CM')
])
%%R -i roc_auc_concat_small -w 550 -h 550

colors = c(
    'Train: TMD+TMR vs CM'='#999999',
    'Test: TMS vs CM'='navy'
)

roc_auc_tmdr_cm(roc_auc_concat_small, annotation=0.15, colors=colors)<jupyter_output><empty_output><jupyter_text>This time, there is no need to average anything as we got proper test results avilable for all the models:<jupyter_code>small_lasso_concat.cross_validation.contributions.combined.add_weighted_auc(
    small_lasso_concat.sub_sampling_test_results.cv_auc
)
small_lasso_contribs = small_lasso_concat.cross_validation.contributions.combined.data
%%R -w 700 -i small_lasso_contribs
coefficients_volcano_plot(
    small_lasso_contribs,
    n=15, p_value='weighted_auc', neg=F, point_color='grey40',
    highlight_frequent=0.75
) + ylab('Mean contribution value')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09108785283856645, 'fraction_numerical': 0.007865524896923565, 'mean_word_length': 5.7869134739560915, 'pattern_counts': {'":': 0, '<': 72, '<?xml version=': 0, '>': 71, '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': '18', 'dir': 'jupyter-structured-clean-dedup', 'id': '49695', 'license': 'MIT', 'n_tokens_mistral': 5415, 'n_tokens_neox': 5024, 'n_words': 1137, 'path': 'analyses/integration/Logistic_LASSO-B-concatenated.ipynb', 'repo_name': 'krassowski/meningitis-integration', 'size': '15765'}