Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion src/rdds/gicam/__main__.py
Original file line number Diff line number Diff line change
Expand Up @@ -79,7 +79,7 @@ def _explore(args):
help='Number of CPU cores to allocate for processing',
default=cpu_count() - 1)
subparser.add_argument('--replace_overwrite_vrs',
help='Write GICAM inference value to VrsModelPrediction field instead of separate GICAM' +
help='Write GICAM inference value to MivmirScore field instead of separate GICAM' +
'(not to be used in production)',
default=False)
def _infer_vcf(args):
Expand Down
2 changes: 1 addition & 1 deletion src/rdds/gicam/model/model.py
Original file line number Diff line number Diff line change
Expand Up @@ -468,7 +468,7 @@ def score_variants(self, variants: List[Variant]) -> np.ndarray:
score_mivmir: np.ndarray = np.zeros((n_samples, 1)) # [batch_dim, feature_dim]
score_genmod: np.ndarray = np.zeros((n_samples, 1))
for i, variant in enumerate(variants):
score_mivmir[i, 0] = variant.INFO['VrsModelPrediction']
score_mivmir[i, 0] = variant.INFO['MivmirScore']
rank_score_normalized_str: str = variant.INFO['RankScoreNormalized'] # format: str: case_name:rank_score
rank_score = float(rank_score_normalized_str.split(':')[1])
score_genmod[i, 0] = rank_score
Expand Down
10 changes: 5 additions & 5 deletions src/rdds/gicam/vcf_inference/infer_vcf.py
Original file line number Diff line number Diff line change
Expand Up @@ -26,15 +26,15 @@ def _infer_gicam_fn(vcf_file_path: str,

vcf_reader = VCFReader(vcf_file_path)
if not replace_overwrite_vrs_annotation:
vcf_reader.add_info_to_header({'ID': 'GICAM',
'Description': 'Rank score from GICAM model (joint MIVMIR and Genmod) (5 points precision)',
vcf_reader.add_info_to_header({'ID': 'GicamScore',
'Description': 'Rank score from GICAM model (joint MIVMIR and Genmod genetic models) (5 points precision)',
'Type': 'Float',
'Number': '1'})
# TODO: Add GICAM version to header
# Make a copy of the input VCF which is also the output file
subprocess_output_file_name = os.path.join(subprocess_work_dir, f'{variant_index_start}.vcf')
vcf_writer = VcfWriter(subprocess_output_file_name,
vcf_reader, # Reuse original file header, with VrsModelPrediction appended
vcf_reader, # Reuse original file header, with MivmirScore appended
mode='w')
# Load and preprocess variants
# Force load complete VCF into RAM as list of variants, drop out of scope variants
Expand All @@ -47,9 +47,9 @@ def _infer_gicam_fn(vcf_file_path: str,
scores = gicam.score_variants(variants=variants)
for i, (variant, score) in enumerate(zip(variants, scores)):
if replace_overwrite_vrs_annotation:
variant.INFO['VrsModelPrediction'] = f'{score:.5f}'
variant.INFO['MivmirScore'] = f'{score:.5f}'
else:
variant.INFO['GICAM'] = f'{score:.5f}'
variant.INFO['GicamScore'] = f'{score:.5f}'
vcf_writer.write_record(variant)
vcf_writer.close()
vcf_reader.close()
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -209,9 +209,9 @@
"def visualize_numerical_feature_vs_rank_score(df, annotation):\n",
" fig = plt.figure(figsize=FIGSIZE)\n",
" ax = fig.add_subplot()\n",
" sb.scatterplot(df[[annotation, 'VrsModelPrediction']],\n",
" sb.scatterplot(df[[annotation, 'MivmirScore']],\n",
" x=annotation,\n",
" y='VrsModelPrediction',\n",
" y='MivmirScore',\n",
" ax=ax)\n",
" fig.tight_layout()\n",
" fig.suptitle(f'n={len(df[annotation].dropna())} ({len(df)})')\n",
Expand Down Expand Up @@ -269,15 +269,15 @@
" data = df[[annotation]]\n",
" for index, row in data.iterrows():\n",
" sentence = _preprocess_text(row[annotation])\n",
" _preprocessed_df = pd.DataFrame(data={'VrsModelPrediction': df.loc[index].VrsModelPrediction,\n",
" _preprocessed_df = pd.DataFrame(data={'MivmirScore': df.loc[index].MivmirScore,\n",
" annotation: sentence,\n",
" 'case_name': df.loc[index].case_name},\n",
" index=[index])\n",
" if preprocessed_df is None:\n",
" preprocessed_df = _preprocessed_df\n",
" else: \n",
" preprocessed_df = pd.concat((preprocessed_df, _preprocessed_df), axis=0)\n",
" facet = sb.relplot(data=preprocessed_df, x='case_name', y='VrsModelPrediction', hue=annotation, style=annotation, kind='scatter', height=FIGSIZE[1],\n",
" facet = sb.relplot(data=preprocessed_df, x='case_name', y='MivmirScore', hue=annotation, style=annotation, kind='scatter', height=FIGSIZE[1],\n",
" aspect = FIGSIZE[0] / FIGSIZE[1])\n",
" facet.set(xticklabels=[])\n",
" facet.set(title=f'n={len(preprocessed_df.dropna())} ({len(df)})')\n",
Expand All @@ -292,13 +292,13 @@
"outputs": [],
"source": [
"def print_nan_entries(df, annotation):\n",
" nan_entries = df.loc[df[annotation].isna()][['VrsModelPrediction', annotation]]\n",
" nan_entries = df.loc[df[annotation].isna()][['MivmirScore', annotation]]\n",
" if len(nan_entries) > 0:\n",
" print(nan_entries, flush=True)\n",
" mean = nan_entries.VrsModelPrediction.mean()\n",
" std = nan_entries.VrsModelPrediction.std()\n",
" mean = nan_entries.MivmirScore.mean()\n",
" std = nan_entries.MivmirScore.std()\n",
" print(f'mean {mean:.4f} std {std:.4f}')\n",
" return nan_entries.VrsModelPrediction"
" return nan_entries.MivmirScore"
]
},
{
Expand All @@ -316,8 +316,8 @@
" 'Compounds_family_id',\n",
" 'CompoundsNormalized',\n",
" 'CSQ_ENSP',\n",
" 'VrsModelExplanation',\n",
" 'VrsModelPrediction']\n",
" 'MivmirExplanation',\n",
" 'MivmirScore']\n",
"\n",
"annotation_nan_statistics = {} # Keep track of pathogenicity scores for empty annotations\n",
"\n",
Expand Down Expand Up @@ -357,9 +357,9 @@
" if d < 0.2 : # if bias are less than\n",
" continue\n",
" facet = sb.displot(data=nan_scores, kind='hist', binwidth=0.1)\n",
" facet.set_xlabels('VrsModelPrediction Score')\n",
" facet.set_xlabels('MivmirScore Score')\n",
" facet.set(xlim=(-0.1, 1.1))\n",
" facet.set(title=f'VrsModelPrediction\\non {annotation} NaNs n={len(nan_scores)}')\n",
" facet.set(title=f'MivmirScore\\non {annotation} NaNs n={len(nan_scores)}')\n",
" facet.tight_layout()"
]
},
Expand Down Expand Up @@ -411,11 +411,11 @@
" df_annotation_info = _df_annotation_info\n",
" else:\n",
" df_annotation_info = pd.concat((df_annotation_info, _df_annotation_info), axis=1)\n",
"df_annotation_info['VrsModelPrediction'] = df.VrsModelPrediction\n",
"df_annotation_info['MivmirScore'] = df.MivmirScore\n",
"df_annotation_info['case_name'] = df.case_name\n",
"\n",
"for annotation in annotations:\n",
" facet = sb.relplot(data=df_annotation_info, x=f'{annotation}%', y='VrsModelPrediction', kind='scatter', height=FIGSIZE[1],\n",
" facet = sb.relplot(data=df_annotation_info, x=f'{annotation}%', y='MivmirScore', kind='scatter', height=FIGSIZE[1],\n",
" aspect = FIGSIZE[0] / FIGSIZE[1])\n",
" facet.set(title=f'Annotation magnitude {annotation}')\n",
" facet.tight_layout()"
Expand All @@ -429,14 +429,14 @@
"outputs": [],
"source": [
"annotation_magnitudes = df_annotation_info.iloc[:, :-5].sum(axis=1)\n",
"facet = sb.relplot(x=annotation_magnitudes, y=df_annotation_info.VrsModelPrediction, hue=df_annotation_info.case_name, style=df_annotation_info.case_name, kind='scatter', height=FIGSIZE[1],\n",
"facet = sb.relplot(x=annotation_magnitudes, y=df_annotation_info.MivmirScore, hue=df_annotation_info.case_name, style=df_annotation_info.case_name, kind='scatter', height=FIGSIZE[1],\n",
" aspect = FIGSIZE[0] / FIGSIZE[1])\n",
"for index, row in df_annotation_info.iterrows():\n",
" facet.ax.text(annotation_magnitudes.loc[index], row.VrsModelPrediction, row.case_name, fontsize=6, rotation=40)\n",
" facet.ax.text(annotation_magnitudes.loc[index], row.MivmirScore, row.case_name, fontsize=6, rotation=40)\n",
"facet.set(xlabel='Total Annotation Magnitude')\n",
"facet.set(title=f'Overall annotation magnitude')\n",
"facet.tight_layout()\n",
"pd.concat((annotation_magnitudes, df_annotation_info[['case_name', 'VrsModelPrediction']]), axis=1).sort_values(0)"
"pd.concat((annotation_magnitudes, df_annotation_info[['case_name', 'MivmirScore']]), axis=1).sort_values(0)"
]
},
{
Expand All @@ -450,13 +450,13 @@
"df_investigate_framshift_missense = pd.DataFrame()\n",
"df_investigate_framshift_missense['consequence'] = df.CSQ_Consequence\n",
"df_investigate_framshift_missense['annotation_magnitude'] = annotation_magnitudes\n",
"df_investigate_framshift_missense['VrsModelPrediction'] = df.VrsModelPrediction\n",
"df_investigate_framshift_missense['MivmirScore'] = df.MivmirScore\n",
"df_investigate_framshift_missense['case_name'] = df.case_name\n",
"df_investigate_framshift_missense.sort_values('consequence', inplace=True)\n",
"\n",
"facet = sb.relplot(y=df_investigate_framshift_missense.consequence,\n",
" x=df_investigate_framshift_missense.annotation_magnitude,\n",
" hue=df_investigate_framshift_missense.VrsModelPrediction,\n",
" hue=df_investigate_framshift_missense.MivmirScore,\n",
" height=FIGSIZE[1],\n",
" aspect=FIGSIZE[0] / FIGSIZE[1])\n",
"for index, row in df_investigate_framshift_missense.iterrows():\n",
Expand Down Expand Up @@ -626,7 +626,7 @@
"metadata": {},
"outputs": [],
"source": [
"points = df[['VrsModelPrediction','case_name']]\n",
"points = df[['MivmirScore','case_name']]\n",
"points"
]
},
Expand Down Expand Up @@ -671,7 +671,7 @@
"num_plot = cosmo(\n",
" points = points,\n",
" point_label_by = 'case_name',\n",
" point_color_by = 'VrsModelPrediction',\n",
" point_color_by = 'MivmirScore',\n",
" links = links_numerical,\n",
" link_source_by = 'sources',\n",
" link_target_by = 'targets',\n",
Expand Down Expand Up @@ -780,7 +780,7 @@
"txt_plot = cosmo(\n",
" points = points,\n",
" point_label_by = 'case_name',\n",
" point_color_by = 'VrsModelPrediction',\n",
" point_color_by = 'MivmirScore',\n",
" links = links_text,\n",
" link_source_by = 'sources',\n",
" link_target_by = 'targets',\n",
Expand Down Expand Up @@ -835,7 +835,7 @@
"num_text_plot = cosmo(\n",
" points = points,\n",
" point_label_by = 'case_name',\n",
" point_color_by = 'VrsModelPrediction',\n",
" point_color_by = 'MivmirScore',\n",
" links = links_numerical_text,\n",
" link_source_by = 'sources',\n",
" link_target_by = 'targets',\n",
Expand Down Expand Up @@ -871,14 +871,14 @@
"outputs": [],
"source": [
"facet = sb.relplot(x=df.case_name,\n",
" y=df.VrsModelPrediction,\n",
" y=df.MivmirScore,\n",
" hue=df.case_name,\n",
" style=df.case_name,\n",
" height=FIGSIZE[1],\n",
" aspect = FIGSIZE[0] / FIGSIZE[1],\n",
" palette='rocket')\n",
"for index, row in df.iterrows():\n",
" facet.ax.text(row.case_name, row.VrsModelPrediction, row.case_name, fontsize=6, rotation=40)\n",
" facet.ax.text(row.case_name, row.MivmirScore, row.case_name, fontsize=6, rotation=40)\n",
"facet.set(xticklabels=[])"
]
},
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -73,7 +73,7 @@ def view_vcf_rank_results(vcf_file_path: str,
vrs_rank_score = np.empty(n_variants)
frq = np.empty(n_variants)
vrs_model_explanation = np.empty(n_variants, dtype=object)
parse_only_fields = ['POS', 'CHROM', 'RankScore', 'RankScoreNormalized', 'VrsModelPrediction', 'Frq']
parse_only_fields = ['POS', 'CHROM', 'RankScore', 'RankScoreNormalized', 'MivmirScore', 'Frq']
pbar = ProgressBar(max_value=n_variants)
_LOGGER.info(f'{vcf_file_path}, n={n_variants} variants')
for i, variant in enumerate(variants):
Expand All @@ -97,12 +97,12 @@ def view_vcf_rank_results(vcf_file_path: str,
genmod_rank_score_normalized[i] = parsed_variant.RankScoreNormalized_value
else:
genmod_rank_score_normalized[i] = np.nan
if 'VrsModelPrediction' in parsed_variant.parsed_fields:
vrs_rank_score[i] = parsed_variant.VrsModelPrediction
if 'MivmirScore' in parsed_variant.parsed_fields:
vrs_rank_score[i] = parsed_variant.MivmirScore
else:
vrs_rank_score[i] = np.nan
try:
vrs_model_explanation[i] = variant.INFO['VrsModelExplanation']
vrs_model_explanation[i] = variant.INFO['MivmirExplanation']
except KeyError:
pass
pbar.update(i)
Expand Down
14 changes: 7 additions & 7 deletions src/rdds/variant_rank_score/vcf_inference/predict_on_vcf.py
Original file line number Diff line number Diff line change
Expand Up @@ -20,19 +20,19 @@ def _subprocess_predict_on_vcf_part(vcf_file_path: str,
variant_index_stop: int):
from ..model import VariantRankScoreModel
vcf_reader = VCFReader(vcf_file_path)
vcf_reader.add_info_to_header({'ID': 'VrsModelPrediction',
'Description': 'Rank score from VRS model (5 points precision)',
vcf_reader.add_info_to_header({'ID': 'MivmirScore',
'Description': 'Rank score from MIVMIR model (5 points precision)',
'Type': 'Float',
'Number': '1'})
vcf_reader.add_info_to_header({'ID': 'VrsModelExplanation',
'Description': 'List of annotation impact scores on VrsModelPrediction (2 points precision)',
vcf_reader.add_info_to_header({'ID': 'MivmirExplanation',
'Description': 'List of annotation impact scores on MivmirScore (2 points precision)',
'Type': 'String',
'Number': '.'})

# Make a copy of the input VCF which is also the output file
subprocess_output_file_name = os.path.join(subprocess_work_dir, f'{variant_index_start}.vcf')
vcf_writer = VcfWriter(subprocess_output_file_name,
vcf_reader, # Reuse original file header, with VrsModelPrediction appended
vcf_reader, # Reuse original file header, with MivmirScore appended
mode='w')

# Load and preprocess variants
Expand All @@ -55,7 +55,7 @@ def _subprocess_predict_on_vcf_part(vcf_file_path: str,
df: pd.DataFrame = vrs_model.score_variant(parsed_variants)
for i, variant in enumerate(variants):
df_i = df.iloc[i]
variant.INFO['VrsModelPrediction'] = f'{df_i.pathogenicity_score:.5F}'
variant.INFO['MivmirScore'] = f'{df_i.pathogenicity_score:.5F}'
# Sort the explanations in decreasing importance (positive = more contributing to higher scoring result)
explanations_sorted_in_decreasing_importance = df_i.sort_values(ascending=False)
vrs_model_explanations = '['
Expand All @@ -66,7 +66,7 @@ def _subprocess_predict_on_vcf_part(vcf_file_path: str,
continue
vrs_model_explanations += f'{key}={contribution_score:.2F},'
vrs_model_explanations += ']'
variant.INFO['VrsModelExplanation'] = vrs_model_explanations
variant.INFO['MivmirExplanation'] = vrs_model_explanations
vcf_writer.write_record(variant)
vcf_writer.close()
vcf_reader.close()
Expand Down
10 changes: 5 additions & 5 deletions src/tests/gicam/inference_test_data.vcf
Original file line number Diff line number Diff line change
Expand Up @@ -84,8 +84,8 @@
##INFO=<ID=SomaticMutationsHeteroplasmy,Number=A,Type=String,Description="The variant was found as heteroplasmic in Mitomap Somatic Mutations datasets (Y=yes,N=no)">
##INFO=<ID=SomaticMutationsHomoplasmy,Number=A,Type=String,Description="The variant was found as homoplasmic in Mitomap Somatic Mutations datasets (Y=yes,N=no)">
##INFO=<ID=dbSNP,Number=A,Type=String,Description="dbSNP ID of the variant">
##INFO=<ID=VrsModelPrediction,Number=1,Type=Float,Description="Rank score from VRS model (5 points precision)">
##INFO=<ID=VrsModelExplanation,Number=.,Type=String,Description="List of annotation impact scores on VrsModelPrediction (2 points precision)">
##INFO=<ID=MivmirScore,Number=1,Type=Float,Description="Rank score from Mivmir model (5 points precision)">
##INFO=<ID=MivmirExplanation,Number=.,Type=String,Description="List of annotation impact scores on MivmirScore (2 points precision)">
##INFO=<ID=Annotation,Number=.,Type=String,Description="Annotates what feature(s) this variant belongs to.">
##INFO=<ID=GeneticModels,Number=.,Type=String,Description="':'-separated list of genetic models for this variant.">
##INFO=<ID=ModelScore,Number=.,Type=String,Description="PHRED score for genotype models.">
Expand Down Expand Up @@ -214,6 +214,6 @@
##rs_dbSNP150=rs_dbSNP150 from dbNSFP file
##LoFtool=LoFtool score for gene
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT SAMPLENAME
1 2 1_C_A C A 68 . VrsModelPrediction=0.98754;RankScore=casename:8.0;RankScoreNormalized=casename:1.0;RankScoreMinMax=casename:-29.0:8.0;RankResult=4|1|3 GT:DP:AD:GQ:PL:RNC 0/1:32:13,18:67:68,0,76:..
1 3 1_3_AG_A AG A 40 . VrsModelPrediction=0.00123;RankScore=casename:6.0;RankScoreNormalized=casename:0.9459459459459459;RankScoreMinMax=casename:-29.0:8.0;RankResult=4|1|1 GT:DP:AD:GQ:PL:RNC 1/1:16:3,13:18:40,17,0:..
1 4 1_4_C_G C G 49 . VrsModelPrediction=1.00000;RankScore=casename:8.0;RankScoreNormalized=casename:1.0;RankScoreMinMax=casename:-29.0:8.0;RankResult=4|1|3 GT:DP:AD:GQ:PL:RNC 0/1:12:5,7:49:49,0,60:..
1 2 1_C_A C A 68 . MivmirScore=0.98754;RankScore=casename:8.0;RankScoreNormalized=casename:1.0;RankScoreMinMax=casename:-29.0:8.0;RankResult=4|1|3 GT:DP:AD:GQ:PL:RNC 0/1:32:13,18:67:68,0,76:..
1 3 1_3_AG_A AG A 40 . MivmirScore=0.00123;RankScore=casename:6.0;RankScoreNormalized=casename:0.9459459459459459;RankScoreMinMax=casename:-29.0:8.0;RankResult=4|1|1 GT:DP:AD:GQ:PL:RNC 1/1:16:3,13:18:40,17,0:..
1 4 1_4_C_G C G 49 . MivmirScore=1.00000;RankScore=casename:8.0;RankScoreNormalized=casename:1.0;RankScoreMinMax=casename:-29.0:8.0;RankResult=4|1|3 GT:DP:AD:GQ:PL:RNC 0/1:12:5,7:49:49,0,60:..
10 changes: 5 additions & 5 deletions src/tests/gicam/test_inference.py
Original file line number Diff line number Diff line change
Expand Up @@ -19,9 +19,9 @@ def test_vcf_inference(work_dir, overwrite_vrs_annotation):

infer_vcf(vcf_file_path=test_data_path, cpu_cores=1, replace_overwrite_vrs_annotation=overwrite_vrs_annotation)
reader = VCFReader(output_file, 'r')
target_annotation = 'GICAM'
target_annotation = 'GicamScore'
if overwrite_vrs_annotation:
target_annotation = 'VrsModelPrediction'
target_annotation = 'MivmirScore'
assert target_annotation in reader.info_fields
for variant in list(reader):
assert isinstance(variant.INFO[target_annotation], float)
Expand All @@ -42,10 +42,10 @@ def test_vcf_inference_cli(work_dir, cpu_cores):
shell=True, stderr=sp.STDOUT)

reader = VCFReader(output_file, 'r')
assert 'GICAM' in reader.info_fields
assert 'GicamScore' in reader.info_fields
for variant in list(reader):
assert isinstance(variant.INFO['GICAM'], float)
assert 0 <= variant.INFO['GICAM'] <= 1
assert isinstance(variant.INFO['GicamScore'], float)
assert 0 <= variant.INFO['GicamScore'] <= 1
reader.close()

def test_score_variant():
Expand Down
Loading