[prs] Distance-based PRS with corrected betas #565
Conversation
WalkthroughThe recent changes significantly enhance the proteomics data analysis codebase, introducing new functionality for polygenic risk score calculations and improved data handling through streaming. Comprehensive documentation and logging capabilities have been added, promoting clarity and monitoring. Additionally, the structure of data classes has evolved to better manage covariates, ensuring robustness and accuracy in genetic data processing. Changes
Sequence Diagram(s)sequenceDiagram
participant User
participant API
participant PRSHandler
participant PRSJobData
participant Logger
User->>API: Initiate PRS Job
API->>PRSHandler: Call calculate_prs_scores()
PRSHandler->>PRSJobData: Extract parameters
PRSJobData->>PRSHandler: Return parameter data
PRSHandler->>Logger: Log job initiation
PRSHandler->>PRSHandler: Process data
PRSHandler->>Logger: Log processing details
PRSHandler->>API: Return results
User->>API: Receive PRS results
Thank you for using CodeRabbit. We offer it for free to the OSS community and would appreciate your support in helping us grow. If you find it useful, would you consider giving us a shout-out on your favorite social media? TipsChatThere are 3 ways to chat with CodeRabbit:
Note: Be mindful of the bot's finite context window. It's strongly recommended to break down tasks such as reading entire modules into smaller chunks. For a focused discussion, use review comments to chat about specific files and their changes, instead of using the PR comments. CodeRabbit Commands (invoked as PR comments)
Additionally, you can add CodeRabbit Configuration File (
|
![coderabbitai[bot]](https://avatars.githubusercontent.com/in/347564?s=60&v=4){kind=small}
There was a problem hiding this comment.
Actionable comments posted: 4
Outside diff range, codebase verification and nitpick comments (2)
python/python/bystro/prs/preprocess_for_prs.py (2)
512-555: Improve inline comments for clarity inprune_by_window.The inline comments could be improved for better clarity and understanding.
Consider revising the comments to provide more context on the pruning process.
Line range hint
632-954: Enhance logging ingenerate_c_and_t_prs_scores.The logging statements could be more descriptive to aid in debugging.
Consider adding more context to the logging messages to provide better insights during execution.
Review details
Configuration used: CodeRabbit UI
Review profile: CHILL
Files selected for processing (7)
- python/python/bystro/examples/ProteomicsProxiedQuery.ipynb (4 hunks)
- python/python/bystro/prs/handler.py (3 hunks)
- python/python/bystro/prs/listener.py (2 hunks)
- python/python/bystro/prs/messages.py (3 hunks)
- python/python/bystro/prs/preprocess_for_prs.py (16 hunks)
- python/python/bystro/prs/prs_types.py (1 hunks)
- python/python/bystro/prs/tests/test_preprocess_for_prs.py (10 hunks)
Additional comments not posted (13)
python/python/bystro/prs/prs_types.py (1)
11-11: Enhancement: New attributes added toPRSJobData.The addition of
covariates_pathandtraining_populationsenhances the flexibility of thePRSJobDataclass. The default values ensure backward compatibility.Also applies to: 14-14
python/python/bystro/prs/listener.py (1)
5-5: Enhancement: Logging functionality added.The logging setup with a DEBUG level is appropriate for tracking operations and diagnosing issues. This addition will aid in monitoring the execution of the PRS job listener.
Also applies to: 17-23
python/python/bystro/prs/messages.py (1)
22-23: Enhancement: ExpandedPRSJobDataattributes.The addition of new attributes such as
training_populations,covariates_path, and others enhances the flexibility and specificity of the data being processed. The removal ofpopulations_pathand the update toindex_namereflect a more streamlined approach to handling population data.Also applies to: 31-38, 47-57
python/python/bystro/prs/handler.py (1)
4-4: Enhancement: Expanded functionality incalculate_prs_scores.The integration of new parameters and covariate data into the scoring process enhances the function's accuracy and flexibility. The changes to the output file naming convention improve clarity and traceability.
Also applies to: 26-33, 42-47, 53-65, 69-71
python/python/bystro/prs/tests/test_preprocess_for_prs.py (1)
26-37: Enhancement: Updated test suite forpreprocess_for_prs.The addition of new fixtures and updates to existing ones ensure that the tests align with the current requirements of the codebase. The removal of a test indicates a potential refactor or removal of functionality.
Also applies to: 50-54, 67-73, 160-171, 269-360
python/python/bystro/prs/preprocess_for_prs.py (6)
1-17: Docstring improvements look good!The docstring provides a clear overview of the expected input format and the purpose of the file.
172-173: Function_load_association_scoreslooks good!The function is well-structured and uses constants for column names and types, improving maintainability.
223-238: Function_convert_rsid_to_query_formatlooks good!The function is straightforward and correctly handles empty input with a ValueError.
Line range hint
244-286: Function_extract_af_and_loci_overlaplooks good!The function is well-structured, and the use of constants for assembly checks improves readability.
Line range hint
288-310: Function_calculate_ancestry_weighted_aflooks good!The function is well-implemented, and the use of DataFrame operations enhances performance.
313-344: Function_preprocess_scoreslooks good!The function is well-structured, and the use of constants for column names improves maintainability.
python/python/bystro/examples/ProteomicsProxiedQuery.ipynb (2)
25-26: Import changes look good!The addition of
stream_filefrombystro.api.streamingis appropriate for the streaming functionality.
75-78: Streaming functionality looks good!The use of a generator to stream data line by line is efficient for handling large datasets.
| if row[ALLELE_COMPARISON_COMPUTED_COLUMN] == REVERSE_MATCH_COMPUTED_VALUE: | ||
| for col in row.index: | ||
| if col != "allele_comparison" and row[col] != -1: | ||
| if col != ALLELE_COMPARISON_COMPUTED_COLUMN and row[col] != -1: | ||
| row[col] = 2 - row[col] |
There was a problem hiding this comment.
Optimize dosage adjustment using DataFrame operations.
The loop can be optimized by using DataFrame operations for better performance.
Consider using DataFrame operations to adjust dosages instead of iterating through rows.
| def _prune_scores( | ||
| scores: pd.DataFrame, | ||
| assembly: str, | ||
| p_value_threshold: float, | ||
| min_abs_beta: float, | ||
| max_abs_beta: float, | ||
| ld_window_bp: int | None = None, | ||
| distance_based_cluster: bool = True, | ||
| training_populations: list | None = None, | ||
| ) -> pd.DataFrame: | ||
| """ | ||
| Prune scores based on p-value threshold | ||
| """ | ||
| # prune preprocessed_scores to only include loci with p-values below the threshold | ||
| # and filter down to sites with beta values within range | ||
| scores = scores[scores[P_COLUMN] <= p_value_threshold] | ||
| scores = scores[ | ||
| (scores[BETA_COLUMN].abs() >= min_abs_beta) & (scores[BETA_COLUMN].abs() <= max_abs_beta) | ||
| ] | ||
|
|
||
| if distance_based_cluster: | ||
| if ld_window_bp is None: | ||
| raise ValueError("If distance_based_cluster is True, ld_window_bp must be provided") | ||
|
|
||
| # if debugging against Dave Cutler's PRS calculator, uncomment to keep only those snps that | ||
| # are in the cutler_snp_set, based on the VARIANT_ID_COLUMN | ||
| # preprocessed_scores = preprocessed_scores[preprocessed_scores[VARIANT_ID_COLUMN].isin(_cutler_snp_set)] # noqa | ||
| scores = prune_by_window(scores, ld_window_bp) | ||
|
|
||
| # TODO: check for non-direct match and warn | ||
| return _preprocess_scores(scores) | ||
|
|
||
| if training_populations is None: | ||
| raise ValueError("If distance_based_cluster is False, training_populations must be provided") | ||
|
|
||
| if len(training_populations) > 1: | ||
| raise ValueError( | ||
| ( | ||
| "PRS training data from only one superpopulation " | ||
| f"or population is currently supported: found {training_populations}." | ||
| ) | ||
| ) | ||
| sumstat_population = training_populations[0] | ||
|
|
||
| if sumstat_population in SUPERPOP_TO_POP_MAP: | ||
| sumstat_population = SUPERPOP_TO_POP_MAP[sumstat_population] | ||
|
|
||
| try: | ||
| sumstat_ld_map_path = get_map_file(assembly, sumstat_population) | ||
| except ValueError as e: | ||
| raise ValueError(f"{sumstat_population} is likely not supported. Failed to get map file: {e}") | ||
|
|
||
| # TODO: check for non-direct match and warn | ||
| preprocessed_scores = _preprocess_scores(scores) | ||
|
|
||
| scores_after_c_t, _loci_and_allele_comparison = ld_clump( | ||
| preprocessed_scores, str(sumstat_ld_map_path) | ||
| ) | ||
|
|
||
| return preprocessed_scores.loc[scores_after_c_t.index] |
There was a problem hiding this comment.
Improve error handling in _prune_scores.
The error messages could be more descriptive to provide better user feedback.
Consider providing more context in the error messages to help users understand the issue.
| @@ -40,12 +51,33 @@ | |||
| } | |||
| ], | |||
| "source": [ | |||
| "user = auth.login('email', 'password', 'https://bystro-dev.emory.edu')\n", | |||
| "user = auth.login('akotlar@bu.edu', 'gonxug-xubram-nUbny6', 'https://bystro-dev.emory.edu')\n", | |||
There was a problem hiding this comment.
Avoid hardcoding credentials.
Hardcoding sensitive information like credentials is a security risk.
Consider using environment variables or a secure vault to manage credentials securely.
| { | ||
| "cell_type": "code", | ||
| "execution_count": 9, | ||
| "id": "458ab2e0-3203-4d72-a697-7e0d37d65a4d", | ||
| "metadata": {}, | ||
| "outputs": [ | ||
| { | ||
| "ename": "NameError", | ||
| "evalue": "name 'joined_df' is not defined", | ||
| "output_type": "error", | ||
| "traceback": [ | ||
| "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", | ||
| "\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)", | ||
| "Cell \u001b[0;32mIn[9], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m \u001b[43mjoined_df\u001b[49m\n", | ||
| "\u001b[0;31mNameError\u001b[0m: name 'joined_df' is not defined" | ||
| ] |
There was a problem hiding this comment.
Fix NameError for joined_df.
The variable joined_df is referenced before being defined, indicating a potential bug.
Verify the initialization and usage of joined_df to ensure it is defined before use.
There was a problem hiding this comment.
Actionable comments posted: 0
Review details
Configuration used: CodeRabbit UI
Review profile: CHILL
Files selected for processing (3)
- python/python/bystro/prs/messages.py (3 hunks)
- python/python/bystro/prs/preprocess_for_prs.py (16 hunks)
- python/python/bystro/prs/prs_types.py (1 hunks)
Files skipped from review as they are similar to previous changes (2)
- python/python/bystro/prs/messages.py
- python/python/bystro/prs/prs_types.py
Additional comments not posted (13)
python/python/bystro/prs/preprocess_for_prs.py (13)
Line range hint
17-66: Imports and constants look good.The imports and constants are well-organized and relevant to the file's functionality.
Line range hint
168-173: Function_load_association_scoresis well-implemented.The function correctly loads and returns GWAS summary statistics with harmonized columns and types.
Line range hint
175-193: Function_load_genetic_maps_from_featheris robust and well-implemented.The function includes error handling and logging, enhancing its robustness.
Line range hint
212-221: Function_convert_loci_to_query_formatis well-structured.The function correctly converts loci to a query format and includes error handling for empty input.
223-238: Function_convert_rsid_to_query_formatis well-structured.The function correctly converts rsids to a query format and includes error handling for empty input.
Line range hint
244-284: Function_extract_af_and_loci_overlapis comprehensive and consistent.The function handles unsupported assemblies and uses constants for allele frequency fields.
Line range hint
288-310: Function_calculate_ancestry_weighted_afis well-implemented.The function calculates ancestry-weighted allele frequencies and saves intermediate results for debugging.
313-344: Function_preprocess_scoresis clear and well-structured.The function handles strand issues by generating SNPID columns for effect and reference alleles.
Line range hint
347-368: Function_preprocess_genetic_mapsis straightforward and correct.The function correctly processes genetic maps and returns a dictionary of upper bound values.
Line range hint
371-382: Functioncompare_allelesis clear and readable.The function uses a custom enum for return values, enhancing readability.
387-390: Optimize dosage adjustment using DataFrame operations.The loop can be optimized by using DataFrame operations for better performance.
Consider using DataFrame operations to adjust dosages instead of iterating through rows.
558-617: Improve error handling in_prune_scores.The error messages could be more descriptive to provide better user feedback.
Consider providing more context in the error messages to help users understand the issue.
Line range hint
632-954: Functiongenerate_c_and_t_prs_scoresis comprehensive and well-documented.The function includes detailed docstrings and logging for debugging. The changes align with the PR objectives.
However, ensure that all function calls to
generate_c_and_t_prs_scoresmatch the new signature.Run the following script to verify the function usage:
Verification successful
Function usage matches the new signature.
The
generate_c_and_t_prs_scoresfunction is correctly used in the codebase, with parameters matching the updated function signature.
- Files:
python/python/bystro/prs/handler.pypython/python/bystro/prs/tests/test_preprocess_for_prs.pyScripts executed
The following scripts were executed for the analysis:
Script:
#!/bin/bash # Description: Verify all function calls to `generate_c_and_t_prs_scores` match the new signature. # Test: Search for the function usage. Expect: Only occurrences of the new signature. rg --type python -A 5 $'generate_c_and_t_prs_scores'Length of output: 82
Script:
#!/bin/bash # Description: Verify all function calls to `generate_c_and_t_prs_scores` match the new signature. # Correcting the previous command by removing the incorrect file type specification. rg -A 5 'generate_c_and_t_prs_scores'Length of output: 3039
There was a problem hiding this comment.
Actionable comments posted: 1
Outside diff range, codebase verification and nitpick comments (3)
python/python/bystro/utils/covariates.py (3)
27-35: Clarify the documentation forfrom_path.The documentation should specify the expected format of the file more explicitly, including the delimiter and any assumptions about the data.
Consider adding more details to the docstring for clarity.
39-45: Enhance error handling for header parsing.The error message for a missing
sample_idcould be more informative by including the expected header format.Consider enhancing the error message to guide the user on the correct format.
55-58: Improve error message for empty sample IDs.The error message could be more informative by suggesting possible reasons for an empty
sample_id.Consider enhancing the error message to provide more context.
Review details
Configuration used: CodeRabbit UI
Review profile: CHILL
Files selected for processing (1)
- python/python/bystro/utils/covariates.py (1 hunks)
Additional comments not posted (4)
python/python/bystro/utils/covariates.py (4)
4-14: LGTM!The
Covariatesclass is well-defined with appropriate use ofmsgspec.Structand clear documentation.
17-20: LGTM!The
ExperimentMappingclass is well-defined with appropriate use ofmsgspec.Structand clear attribute definitions.
23-25: LGTM!The
ExperimentMappingsclass is well-structured with a clear purpose and appropriate use ofmsgspec.Struct.
60-71: Verify covariate field matching.Ensure that the covariate fields in the file match those defined in the
Covariatesclass.Run the following script to verify the covariate field matching:
Verification successful
Covariate Field Matching Verified
The fields in the
Covariatesclass (sex,phenotype,population,superpopulation) match those expected in the code snippet. The code correctly checks these fields against the header values. No discrepancies found.Scripts executed
The following scripts were executed for the analysis:
Script:
#!/bin/bash # Description: Verify covariate fields in the file match those in the `Covariates` class. # Test: Search for covariate fields in the file. Expect: All fields match those in `Covariates`. ast-grep --lang python --pattern $'class Covariates { $$$ }'Length of output: 512
Script:
#!/bin/bash # Description: Verify covariate fields in the file match those in the `Covariates` class. # Find the Covariates class and extract its fields ast-grep --lang python --pattern 'class Covariates: $$$' --json | jq -r '.[].captures[].text' | while read -r class_def; do # Extract fields from the Covariates class echo "$class_def" | ast-grep --lang python --pattern 'def __init__(self, $fields): $$$' --json | jq -r '.[].captures[].text' doneLength of output: 115
Script:
#!/bin/bash # Description: Locate the Covariates class definition and its fields in the codebase. # Search for the Covariates class definition in the codebase rg 'class Covariates' -A 10 --jsonLength of output: 2922
| with open(path, "r") as f: | ||
| lines = f.readlines() |
There was a problem hiding this comment.
Optimize file reading.
Using readlines() can be memory-intensive for large files. Consider using a generator or iterating over the file object directly.
Use this diff to optimize file reading:
- with open(path, "r") as f:
- lines = f.readlines()
+ with open(path, "r") as f:
+ lines = (line for line in f)Committable suggestion
‼️ IMPORTANT
Carefully review the code before committing. Ensure that it accurately replaces the highlighted code, contains no missing lines, and has no issues with indentation. Thoroughly test & benchmark the code to ensure it meets the requirements.
| with open(path, "r") as f: | |
| lines = f.readlines() | |
| with open(path, "r") as f: | |
| lines = (line for line in f) |
There was a problem hiding this comment.
Actionable comments posted: 1
Outside diff range, codebase verification and nitpick comments (1)
python/python/bystro/examples/ProteomicsProxiedQuery.ipynb (1)
Line range hint
604-624: Fix NameError forjoined_df.The variable
joined_dfis referenced before being defined, indicating a potential bug. Verify the initialization and usage ofjoined_dfto ensure it is defined before use.
Review details
Configuration used: CodeRabbit UI
Review profile: CHILL
Files selected for processing (1)
- python/python/bystro/examples/ProteomicsProxiedQuery.ipynb (16 hunks)
Additional comments not posted (5)
python/python/bystro/examples/ProteomicsProxiedQuery.ipynb (5)
35-36: Verify the necessity of thestream_fileimport.Ensure that the
stream_fileimport is used in the notebook. If not, consider removing it to avoid unnecessary imports.
Line range hint
86-98: Query logic looks good.The query for variants within 100kb of gene transcription start sites is well-implemented.
96-98: Display logic is correct.The use of
query_result_df.head()to display the dataframe is appropriate.
646-646: Dosage set creation is correct.The creation of a set from
joined_df['dosage']is appropriate, assumingjoined_dfis correctly initialized.
675-688: Visualization logic is correct.The plot of mean normalized sample intensity by dosage is well-implemented and informative.
| "execution_count": 3, | ||
| "id": "66258792-32fa-4daf-a136-c2b37738066c", | ||
| "metadata": {}, | ||
| "outputs": [ |
There was a problem hiding this comment.
Avoid hardcoding credentials.
Hardcoding sensitive information like credentials is a security risk. Consider using environment variables or a secure vault to manage credentials securely.
I've attached dave's output vs what I believe is ours. I will double check by re-running again, but this is roughly the difference I recall. Our PRS values are within 10% of Dave's at most, and the differences could be down to precision potentially. I spent a good deal of time tracking down the differences, and it seemed they likely happened in the final step, when adjusting dosages (the ancestry-adjusted allele frequencies were identical, as were the genotypes). I will look into this more in the future, but this PR is a strict improvement, and likely right.
Archive 2.zip
Summary by CodeRabbit
New Features
Bug Fixes
Documentation