主题
案例：糖尿病并发症 cis-pQTL MR 研究
这是导师课题的完整方法学（遗传学导向的糖尿病并发症靶点优先级分析），也是 方法引擎 各步骤在真实研究中的落地范例。分析用 R 实现，代码分布在
1_r9_mr_diabetes.Rmd … 8_sensitivity_AMD.Rmd八个脚本。以下 Methods 为原文（英文），按分析步骤分节。
Study overview
We performed a genetics-guided target prioritization analysis for diabetic comorbidity using a multi-step workflow implemented in R. The analytical pipeline consisted of: (i) selection of cis protein quantitative trait locus (pQTL) instruments, (ii) two-sample Mendelian randomization (MR) against diabetic complication GWAS datasets, (iii) validation analyses for type 1 and type 2 diabetes using external GWAS resources, (iv) overlap analysis of significant MR findings across diabetic traits, (v) HyPrColoc analysis to evaluate shared causal signals across proteins and disease traits, (vi) phenome-wide association study (PheWAS) profiling of prioritized variants/proteins, (vii) construction of a protein-disease network, and (viii) sensitivity analyses involving age-related macular degeneration (AMD), including both MR and LD score regression (LDSC) analyses. The codebase for these analyses was distributed across the scripts [1_r9_mr_diabetes.Rmd, 2_diabete_validation.Rmd, 3_r9_upsetR.Rmd, 4_r9_hyprcoloc_newdiabete.Rmd, 5_r9_data_clean.Rmd, 6_r9_phewas.R, 7_network.Rmd, 8_sensitivity_AMD.Rmd].
Protein instrument selection
Genetic instruments for circulating proteins were obtained from the local file data/protein_info.csv. From the resulting dataset, the following columns were retained for downstream analysis: exposure name, cis/trans annotation, exposure identifier, SNP identifier, beta, standard error, effect allele, other allele, effect allele frequency, annotated gene, gene consequence, biotype, and pvallog. Only cis instruments were included. Exposure P values were reconstructed as 10^(-pvallog). The exposure sample size was set to 54,219 for all protein instruments. Duplicated SNPs were checked, although repeated SNP identifiers were retained where they corresponded to different protein exposures.
Outcome GWAS data
Outcome GWAS summary statistics for diabetic complications were primarily obtained from local FinnGen R9 summary statistic files. The following files were used: retinopathy (summary_stats_finngen_R9_DM_RETINOPATHY_EXMORE.gz), maculopathy (summary_stats_finngen_R9_DM_MACULOPATHY_EXMORE.gz), neovascular glaucoma (summary_stats_finngen_R9_DM_NEOVASCULAR_GLAUCOMA.gz), nephropathy (summary_stats_finngen_R9_DM_NEPHROPATHY_EXMORE.gz), and neuropathy (summary_stats_finngen_R9_DM_NEUROPATHY.gz). For each file, only SNPs present in the cis-protein instrument set were retained. Outcome columns were renamed to conform to the TwoSampleMR format as follows: ref to other_allele.outcome, alt to effect_allele.outcome, beta to beta.outcome, sebeta to se.outcome, af_alt to eaf.outcome, pval to pval.outcome, rsids to SNP, and nearest_genes to id.outcome. The script-assigned outcome sample sizes were 10,413 for retinopathy, 3,572 for maculopathy, 1,100 for neovascular glaucoma, 4,111 for nephropathy, and 2,843 for neuropathy.
Validation GWAS data for type 1 and type 2 diabetes
Type 2 diabetes validation was performed using Mahajan.NatGenet2018b.T2D-noUKBB.European.txt, corresponding to a European T2D GWAS from Mahajan et al. A harmonized variant identifier (clean_variant) was created in both the exposure and outcome datasets. In the exposure dataset, clean_variant was defined by truncating the original variant identifier before the :imp:v1 suffix. In the outcome dataset, clean_variant was defined as Chr:Pos:NEA:EA. Exposure and outcome datasets were then matched by this field. The outcome sample size for this T2D validation dataset was set to 57,698.
Type 1 diabetes validation was performed using GCST90475661.tsv.gz. Because the file reported odds ratios rather than beta coefficients, beta values were calculated as log(odds_ratio). Standard errors were reconstructed from the reported confidence interval as (log(ci_upper) - log(ci_lower)) / (2 * 1.96). Effect allele frequency was taken from effect_allele_frequency. SNP matching was performed using rsID, and the outcome sample size was derived from the num_cases field after harmonization.
Two-sample Mendelian randomization
All exposure-outcome pairs were harmonized using harmonise_data() from the TwoSampleMR package with action = 2. Only harmonized variants with mr_keep == TRUE were retained for MR estimation. Instrument strength metrics were calculated for each harmonized variant. The proportion of variance explained (R2) was calculated as:
and the F statistic as F = R2 × (samplesize.exposure − 2) / (1 − R2).
MR estimation was implemented as a single-variant Wald ratio analysis for each harmonized SNP. Specifically, the causal effect estimate was calculated as beta.outcome / beta.exposure, and its standard error was approximated as se.outcome / beta.exposure. The 95% confidence interval was defined as walds_ratio ± 1.96 × se_mr, and a two-sided P value was calculated from the z statistic walds_ratio / se_mr. For each outcome, results were corrected for multiple testing using the Benjamini-Hochberg false discovery rate (FDR) procedure (p.adjust(..., method = "fdr")). Significant MR findings were defined as FDR_Adjusted_P < 0.05. Separate significant result files were generated for each trait and subsequently used for overlap and downstream analyses.
Definition and visualization of shared MR-positive signals
To identify protein signals shared across diabetic traits, the significant result files for retinopathy, neuropathy, nephropathy, maculopathy, type 2 diabetes, and type 1 diabetes were converted into binary tables by assigning each dataset a trait-specific indicator column equal to 1. These datasets were then merged by SNP using outer joins, and missing values were replaced with 0. The resulting matrix (significant_sum) was used to summarize overlap patterns across traits. UpSet plots were generated using both UpSetR and ComplexUpset, and the merged overlap table was written for downstream annotation and HyPrColoc selection.
HyPrColoc analysis
HyPrColoc analysis was performed for proteins shared across one or more diabetic traits, using the workflow defined in [4_r9_hyprcoloc_newdiabete.Rmd]. Protein-level GWAS files were stored as compressed .gz files under 1/r9/hyprcoloc/new_gwas_Data. For each candidate protein, a local genomic window spanning ±500 kb around the target gene position (GENPOS) was extracted from the protein GWAS file using the custom subset_gene_region() function.
To identify SNPs shared between the protein and disease GWAS datasets, a unique variant key was created in the protein data as CHROM_GENPOS_ALLELE0_ALLELE1 and in the disease data as chrom_pos_ref_alt. Only SNPs present in the protein dataset and in all disease datasets included in the current combination were retained. The order of SNPs across all matched datasets was checked explicitly, and combinations failing this order check were discarded. For each retained protein-disease combination, beta and standard error matrices were constructed with rows representing overlapping SNPs and columns representing the protein trait plus one or more disease traits. HyPrColoc was then run as hyprcoloc(pro_beta, pro_se, trait.names = traits, snp.id = rsid) using default settings, because no alternative priors or posterior thresholds were passed explicitly in the script.
The HyPrColoc output table (hyprcoloc_results$results) was augmented with protein, disease_combination, and snp_count, and saved as both CSV and RDS files. The fields retained in downstream cleaned summary tables included iteration, traits, posterior_prob, regional_prob, candidate_snp, posterior_explained_by_snp, and dropped_trait. Because the final binary summary of HyPrColoc-positive results was assembled after cleaning and restructuring these outputs in [5_r9_data_clean.Rmd], the operational definition of a HyPrColoc-positive result in the present workflow was based on the cleaned posterior summary outputs rather than a separately scripted hard threshold. HyPrColoc-positive results were defined as colocalization signals with posterior_prob > 0.7.
Phenome-wide association study (PheWAS)
PheWAS analyses were conducted using per-SNP phenotype scan files stored under 1/result/r9/phewas/, with one CSV file per SNP. The master MR-significant dataset significant_mr_all.csv was used as the reference table. For each SNP-specific PheWAS file, alleles were compared against the MR outcome effect allele using the custom compare_alleles() function. If the PheWAS effect allele (EA) matched the MR effect_allele.outcome, the direction of association (MatchType) was assigned directly from the sign of beta.outcome. If the PheWAS non-effect allele (NEA) matched the MR effect allele, the direction was inverted. Results from all SNP files were concatenated into a single table (phewasall_new0821.csv).
For each association, logP = -log10(P.value) was calculated. A Bonferroni threshold was defined as 0.05 / nrow(combined_results), and associations meeting or exceeding this threshold on the -log10(P) scale were annotated as "Strong." Gene symbols were extracted from the exposure field by retaining the substring before the first colon and stored in the final PheWAS result file phewas_sig_gene.csv.
Overall PheWAS visualization
Overall PheWAS plots were generated from phewas_sig_gene.csv. The plotted fields included phenotype domain (Domain), trait label (Trait), gene, P value, and direction (MatchType). Duplicate protein-trait entries were removed. For each phenotype domain and label combination, either the maximum absolute signed signal or the maximum unsigned logP value was retained, depending on the plot version. Jitter was added to x-axis phenotype positions to reduce overplotting. Three main visualizations were produced: a signed PheWAS plot using signed_logP, an unsigned PheWAS plot with all labels, and a reduced-label version restricted to top signals. Direction was encoded by shape (positive = 24, negative = 25) and phenotype domain by color.
Protein-disease network analysis
Protein-disease networks were constructed from phewas_sig_gene.csv. Edges were aggregated by gene and disease. In the main implementation, the mean walds_ratio across rows within each gene-disease pair was used to define edge direction: positive mean values were assigned "+", negative mean values "-", and missing values NA. The number of rows contributing to each gene-disease pair was used as an evidence count. In an alternative implementation, edge direction was assigned by majority vote of MatchType, with ties recorded as missing. Nodes were classified as either Protein or Disease. During downstream interpretation, TRIM40 was removed because it showed inconsistent directions across diseases.
Drug-target PheWAS and directional interpretation
To evaluate therapeutic potential and possible safety liabilities, proteins were separated into positive_targets and negative_targets based on the sign of walds_ratio in the protein-disease summary table. Positive targets were defined as proteins with walds_ratio > 0, and negative targets as proteins with walds_ratio < 0. Within each group, Bonferroni-adjusted and FDR-adjusted P values were calculated using p.adjust(..., method = "bonferroni") and p.adjust(..., method = "BH"), respectively.
To remove proteins with directionally conflicting PheWAS signals, two intermediate binary tables were used: positive_binary_csv.csv and negative_binary_csv.csv. For positive targets, genes were removed if they showed MatchType == "negative" with p_bonf < 0.05. For negative targets, genes were removed if they showed MatchType == "positive" with p_bonf < 0.05. The retained datasets were named positive_targets2 and negative_targets2. Separate PheWAS plots were then generated for these two groups using -log10(p_bonf) as the y-axis, phenotype domain as color, and association direction as point shape. This framework was used to distinguish potentially beneficial pleiotropic effects from potential adverse effects according to the intended direction of target modulation inferred from MR.
AMD sensitivity analysis
Sensitivity analyses involving age-related macular degeneration (AMD) were performed in two stages. First, the same cis-protein instrument selection and single-SNP Wald ratio MR workflow described above was applied to AMD using FinnGen R9 summary statistics from summary_stats_finngen_R9_H7_AMD.gz. Outcome fields were harmonized identically to the main diabetic complication analyses, and the outcome sample size was set to 8,913. Significant AMD-associated proteins were defined as FDR_Adjusted_P < 0.05.
Second, the genetic correlation between AMD and diabetic maculopathy/DME was evaluated using LD score regression implemented via GenomicSEM. AMD summary statistics were read from summary_stats_finngen_R9_H7_AMD.gz, and DME summary statistics were read from summary_stats_finngen_R9_DM_MACULOPATHY_EXMORE.gz. Summary files were reformatted to retain the columns SNP, A1, A2, BETA, SE, P, and N. Variants with missing values, missing rsIDs, or non-ACGT alleles were removed. Reformatted summary statistics were written to text files and munged with munge() using HapMap3 SNPs from w_hm3.snplist. LD score regression was then run with ldsc() using European LD reference files located in eur_w_ld_chr/. Two analyses were performed: one without prevalence correction and one including sample prevalence correction. For the prevalence-adjusted analysis, sample prevalence was calculated as 8913 / (8913 + 348936) for AMD and 3572 / (3572 + 308547) for DME, with population.prev = NA for both traits. The resulting observed-scale heritability matrix, intercept matrix, and genetic correlation estimates were interpreted as the main outputs.
Statistical significance
For MR analyses, significance was defined as FDR_Adjusted_P < 0.05 within each outcome-specific result table. HyPrColoc-positive results were defined as colocalization signals with posterior_prob > 0.7. For overall PheWAS visualization, a Bonferroni threshold was defined as 0.05 / number_of_tests based on the size of the combined PheWAS result table. For drug-target filtering, conflicting secondary phenotype associations were identified using p_bonf < 0.05.
代码脚本对照：
1_r9_mr_diabetes（主 MR）·2_diabete_validation（T1D/T2D 验证）·3_r9_upsetR（跨性状重叠）·4_r9_hyprcoloc_newdiabete（共定位）·5_r9_data_clean（清洗汇总）·6_r9_phewas（PheWAS）·7_network（蛋白-疾病网络）·8_sensitivity_AMD（AMD 敏感性 + LDSC）。方法出处见 参考文献。