目录
- dataReproducibility: offline enrichment, seeded weighted median, prune stale tables13小时前
- manuscriptReframe as genetic prioritisation, add Bonferroni correction and ambient-RNA sensitivity25分钟前
- outReframe as genetic prioritisation, add Bonferroni correction and ambient-RNA sensitivity25分钟前
- scriptsReframe as genetic prioritisation, add Bonferroni correction and ambient-RNA sensitivity25分钟前
- .gitignoreTrack R result tables in out/ (DEGs, WGCNA, scoring, MR, coloc, enrichment, GTEx)15小时前
- LICENSEAdd MIT license19小时前
- README.mdReframe as genetic prioritisation, add Bonferroni correction and ambient-RNA sensitivity25分钟前
- REPRODUCIBILITY.mdCorrect the quantitative enrichment comparison in REPRODUCIBILITY.md13小时前
- manuscript_en_with_figures.mdReframe as genetic prioritisation, add Bonferroni correction and ambient-RNA sensitivity25分钟前
- manuscript_zh_with_figures.mdReframe as genetic prioritisation, add Bonferroni correction and ambient-RNA sensitivity25分钟前
- mr_candidate_genes.csvCode, figures and manuscript for the plaque signature cis-MR study22小时前
- run.shCode, figures and manuscript for the plaque signature cis-MR study22小时前
邀请码
版权所有：中国计算机学会技术支持：开源发展技术委员会
京ICP备13000930号-9
京公网安备 11010802047560号
Gene-set scoring in the human atherosclerotic plaque atlas (Traeuble et al. 2025)
Code, figures and manuscript for this project: https://gitlink.org.cn/yanhaidong/vein-graft-atherosclerosis-atlas
Scores an arbitrary gene list (e.g. a 104-gene signature) across every cell type of the integrated human plaque atlas, on all 259,116 cells, on a laptop with ~8 GB RAM.
Dataset
plaque atlas(CELLxGENE Discover)db70986c-7d91-49fe-a399-a4730be394acX=log1p(CP10k)-normalised counts, CSR (X/indptr,X/indices,X/data)cell_type_level1(13 types),cell_type_level2(23 types, author-level),cell_typecell_type_level1keeps Fibromyocyte separate from Smooth Muscle Cell, which matters for SMC-derived plaque phenotypes.Why this pipeline
SeuratDisk/hdf5rare not installable here, and converting a 3 GB h5ad would needTwo passes:
--subsamplecells total) needed for z-scores and for the AddModuleScore control bins;Runtime on this machine: ~2 min for all 259k cells (peak RSS ≈ 600 MB).
Scores (per cell)
score_meanscore_amSeurat::AddModuleScorescore_z(x - mean_g) / sd_g), library-size/level independentn_genes_detlib_sizescore_amis the primary ranking score;score_zis the most robust to genes with very different detection rates.Usage
1. Get the gene list out of WGCNA
R objects live only inside the session that made them, so export the vector once (
write.csv(data.frame(gene = candidate_genes_final), "candidate_genes_final.csv", row.names = FALSE)) or save the workspace.scripts/02_genes_from_wgcna.Rconverts either form to the one-symbol-per-line list the scorer wants, and pre-checks every symbol against the atlasvarbefore the 2-min scoring run:It prints
matched in atlas: N / M, writesdata/genes.txtwith the matched symbols only, and drops the misses intodata/genes_not_in_atlas.txt(aliases, lncRNAs or probe IDs that the atlas calls differently).2. Score the list
or directly:
Useful flags:
--chunk 5000(cells per block),--subsample 60000(cells for pass 1),--nbin 24(AddModuleScore bins),--no-plots,--max-cells N(smoke test on the first N cells),--seed 1.Outputs (per run,
out/<tag>/)cell_scores.tsv.gzlib_size,n_genes_detgene_mean_by_celltype.tsvgene_pct_by_celltype.tsvscore_summary_by_celltype.tsvscore_summary_console.txtgenes_used.txt,genes_not_in_atlas.txtvar/feature_name)figures/violin_score_by_celltype.{png,pdf}figures/dotplot_genes_by_celltype.pngfigures/heatmap_genes_by_celltype.pngfigures/umap_score.png,figures/umap_celltype.pngsessionInfo.txtSanity check with canonical markers
data/genes_demo_markers.txt(29 lineage markers) was used to validate the pipeline. It recovers exactly the expected biology (full atlas,cell_type_level1,score_am):Per-gene dotplot confirms the drivers (SMC:
ACTA2/MYH11/CNN1/TAGLN, fibroblast:COL1A1/COL3A1/COL5A1/DCN/LUM/POSTN, macrophage:CD68/C1QA/CD163/APOE, EC:PECAM1/VWF/CDH5/KDR/ACKR1), so the extraction and the CSR indexing are correct.Result: the 104-gene WGCNA module
Input:
data/candidate_genes_final.csv(104 symbols, NF-kB / cytokine / ER-stress / apoptosis / lipid-handling programme). 103 of 104 are present in the atlasvar;SOD2is genuinely absent (the atlas carriesSOD1,SOD2-OT1,SOD3), so it is listed indata/genes_not_in_atlas.txtand the run uses 103 genes (out/genes__cell_type_level1).cell_type_level1(13 author labels, all 259,116 cells)Kruskal-Wallis across cell types: chi2 = 65,230, df = 12, p < 2.2e-16 (
score_summary_console.txt).score_amandscore_zagree on the top 3; SMC and fibromyocyte are the lowest of all 13 types, i.e. the module is not a smooth-muscle / fibroblast programme.cell_type_level2(23 labels,out/genes__cell_type_level2)Within the monocyte-macrophage compartment the lipid-loaded / inflammatory
PLIN2+/TREM1+subset is the highest scorer and has the largestscore_zof any label, whileTREM2+/FoamyandHMOX1+macrophages are close to the macrophage average — the module tracks the inflammatory rather than the foamy/efferocytotic phenotype.Per-gene drivers (full atlas means,
gene_mean_by_celltype.tsv)CXCL84.71 (93% of cells),NFKBIA2.75 (86%),FOS2.22 (71%),IL1B1.91 (59%),NFE2L21.41 (56%), thenICAM1,HSP90AA1,MAP2K3.NFKBIA2.30,CXCL82.05,FOS2.01,IL1B1.90,CXCL21.48.FOS1.89,HSP90AA11.60,NFKBIA1.48,CALM11.40,HSP90AB11.27.CALM11.76,HSP90AB11.65,HSP90AA11.43), i.e. no module-specific activation.Controls (both in
out/)CXCL8/NFKBIA/FOS/IL1B/NFE2L2(_module_minus5)_null_random103)All 13 labels side by side:
out/validation_summary.tsv.cor(score_am, n_genes_detected) = 0.044andcor(score_am, log10 lib_size) = -0.018across all cells, and -0.02 / 0.01 / 0.05 within neutrophil / monocyte / macrophage — the control bins removed the depth effect.CXCL8/IL1B/NFKBIA/FOS/NFE2L2), so this is biology, not a single-gene outlier — but it does mean the neutrophil ranking is driven by the immediate-early / NF-kB arm of the module.Bottom line
The module is an inflammatory myeloid (NF-kB / cytokine / ER-stress) programme. It is highest in neutrophils (1.5% of the atlas, 97% of cells positive, 20.7x over-representation in the global top 1%), then in monocytes and macrophages, with the
PLIN2+/TREM1+inflammatory macrophage as the top-scoring myeloid subset. It is lowest in SMC / fibromyocyte / fibroblast, so it does not support an SMC- or fibroblast-derived origin for this WGCNA module.MR candidate genes, Figure 3, and pathway enrichment
1.
mr_candidate_genes.csv— 18 genes for MR / target follow-upscripts/04_mr_candidates.Rcombines three sources:out/genes__cell_type_level2(myeloid mean vs non-myeloid mean, log2FC, % of myeloid cells positive);data/dgidb_annotation.tsv; 101/103 genes annotated, 78/103 are in theDRUGGABLE GENOMEcategory);CXCL8,IL1B,NFKBIA,FOS,NFE2L2) + the five strongest myeloid-enriched drivers + Druggable Genome genes with myeloid log2FC >= 1, >= 10% positive cells and mean >= 0.05, until 18 genes.Expression columns are mean log1p(CP10k) in neutrophils / monocytes / the best macrophage subset (from
cell_type_level2).representative_agentsin the CSV lists approved or clinical-stage agents per gene (canakinumab for IL1B, dimethyl fumarate for NFE2L2, TZDs for PPARG, …).Notes:
SOD2is not in the atlas, so it cannot appear in any expression-filtered list (data/genes_not_in_atlas.txt).NFE2L2is not annotatedDRUGGABLE GENOMEby DGIdb (no direct small molecule; dimethyl fumarate/bardoxolone act on the NRF2 pathway), so itsdruggable_genomeflag isFALSEalthough itsrepresentative_agentsis populated.JAK2(0.39),CHUK(0.48),IKBKG,TBK1(0.34),IRAK4,MAPK13(1.51 but only 14% of cells),CASP3,BCL2,HSP90AA1.out/mr_candidates_all_genes.tsv, so the cut-offs can be re-tuned with--n-max,--min-pct,--min-mean,--min-lfc.These 18 genes are what the cis-MR below uses; the actual MR estimates are in Cis-MR: do the candidate genes causally affect CAD?.
2. Figure 3 — core genes at
cell_type_level2scripts/03_core_gene_figures.RstreamsXonce (same chunked CSR reader as01_score_genes.R) and keeps only the 18 candidate columns, so the whole figure costs ~40 MB of RAM and one pass over the 3 GB file.Outputs in
out/core_genes__cell_type_level2/figures/:fig3A_violin_core_genes.png/.pdffig3B_umap_core_genes.pngfig3C_umap_celltype_reference.pngfig3_core_genes_combined.pngData behind the panels:
core_gene_expression.tsv.gz(all 259,116 cells x 18 genes),core_gene_mean_by_celltype.tsv,core_gene_pct_by_celltype.tsv,core_gene_overall_summary.tsv(mean, % positive, top cell type, myeloid log2FC).3. GO / KEGG enrichment —
scripts/05_enrichment.pygseapy
enrichrwith the atlas’ own gene universe as background (data/atlas_background_genes.txt, 25,377 symbols), so the Fisher exact test is comparable toclusterProfiler::enricherrather than Enrichr’s generic background. Libraries: GO BP/MF/CC (2021), KEGG 2021, Reactome 2022, WikiPathways 2024, MSigDB Hallmark 2020. The seven GMT files are committed underdata/gmt/and read locally, so the test runs offline and gives identical P/FDR values on every machine.Key pathways (best term per family, KEGG 2021 unless stated):
Outputs:
out/enrichment/enrichment_all.tsv(every term, every library),enrichment_key_pathways.tsv,fig_enrichment_barplot.png/.pdf(A: top 5 terms per library, B: NF-κB / TNF / cytokine-receptor axis).Interpretation caveat. The 104-gene module was itself assembled from NF-κB, cytokine, ER-stress, apoptosis and lipid programmes, so these enrichments confirm what the module is — they are not independent evidence for the phenotype. The independent evidence is the cell-type scoring in the result section above. Two further points worth knowing: all 103 genes fall inside KEGG “Lipid and atherosclerosis” (a 194-gene term in this background), i.e. the module is essentially a subset of that pathway; and terms named after infectious diseases (Hepatitis B, measles, Kaposi) are driven by the same shared cytokine genes, so they are hidden from the figure (still present in
enrichment_all.tsv).Cis-MR: do the candidate genes causally affect CAD?
scripts/06_fetch_mr_inputs.sh(data) +scripts/07_run_mr.R(TwoSampleMR) run a two-sample cis-MR of the 18 candidate genes on coronary artery disease. Both summary-statistic sources are streamed and cut down to the rows we need, so the whole analysis runs locally from public files and nothing genome-wide has to be kept on disk.2019-12-11-cis-eQTLsFDR0.05-ProbeLevel-...txt.gz)GTEx_Analysis_v8_eQTL.tar, tissues Whole_Blood, Artery_Coronary, Artery_Tibial)cad.add.160614.website.txt)Two details that matter and are easy to get wrong:
beta = Z / sqrt(2p(1-p)(N+Z^2)),se = 1 / sqrt(2p(1-p)(N+Z^2)).IL8and ERO1A asERO1L.scripts/06_fetch_mr_inputs.shmaps the old names and writes the current HGNC symbol back into the table. Without this, two perfectly good genes look like they have no cis-eQTL.Instruments are cis SNPs with p < 5e-8 (GTEx: each gene’s own permutation threshold), distance-pruned at 250 kb — LD clumping is switched on only when a local PLINK binary and 1000G
--bfileare supplied (--plink,--bfile), so the default path is offline. All instruments have F > 16, i.e. no weak-instrument problem.Results, per 1-SD higher genetically predicted expression
CYBBis X-linked and GCST003116 contains autosomes only, so it has no outcome data. Diagnostics are inout/mr/mr_report.md: MR-Egger intercepts are all non-significant (no directional pleiotropy), Steiger directionality is correct for every gene, and heterogeneity is limited to CXCL1 (Q p = 0.046) and ICAM1 (p = 0.014), both with 5 and 2 SNPs.Reading of the two hits: higher NFKBIA (IκBα) expression is associated with lower CAD risk and higher CXCL8 (IL-8) expression with higher risk. Both are consistent with the atlas result (the module is NF-κB/cytokine driven) and survive the multiple-testing look (17 genes, Bonferroni 0.05/17 = 2.9e-3) only for NFKBIA — CXCL8 at p = 9.2e-3 does not. The weighted-median estimate for CXCL8 agrees (OR 1.06, p = 0.03), so the signal is not driven by a single SNP.
Variant-level follow-up: the Colocalisation section asks whether the eQTL and CAD signals at these two loci share a causal variant. Neither locus carries a genome-wide significant CAD association (min P = 2.7e-4 and 7.2e-4), so coloc is underpowered there (PP.H4 = 0.03 and 0.20) — MR can still be significant because it aggregates small per-SNP effects.
Outputs:
out/mr/mr_report.md(tables + embedded figures),mr_results_main.tsv,mr_results_all_methods.tsv,mr_instruments_summary.tsv,mr_harmonised_all.tsv,mr_leaveoneout.tsv,mr_singlesnp.tsv,mr_heterogeneity.tsv,mr_pleiotropy.tsv,mr_directionality.tsv(the last five only for genes withColocalisation: do the eQTL and CAD signals share a variant?
cis-MR says whether the transcript matters; it cannot say which variant carries the effect.
scripts/08_coloc.R(coloc 5.2.3) therefore asks, for the two MR hits, whether the cis-eQTL and the CAD association in the same window are consistent with one shared causal variant (PP.H4) or with two distinct causal variants (PP.H3).2019-12-11-cis-eQTLsFDR-ProbeLevel-...txt.gz, 3.9 GB)data/mr/_raw/)Why the MR input tables cannot be reused
coloc.abf()is only valid on complete regional summary statistics.data/mr/eqtlgen_cis.tsvkeeps FDR<0.05 pairs only anddata/mr/cad_subset.tsvkeeps the 12,454 instrument SNPs only; running coloc on those truncated tables inflates the eQTL likelihood and makes PP.H4 meaningless.scripts/06_fetch_mr_inputs.sh --coloctherefore streams the full eQTLGen file (cached underdata/coloc/_raw/, 3.9 GB, git-ignored, safe to delete) and extracts every SNP in the window from both sources.Both traits are harmonised on rsID: every eQTLGen Z-score is converted to beta/SE with
beta = Z / sqrt(2·MAF·(1-MAF)·(N + Z^2))(Vosa et al. 2021), both betas are put on the CAD effect allele (complementing strand-flipped pairs), palindromic SNPs with MAF > 0.42 are dropped, and each gene is tested withcoloc.abf(p12 = 1e-5)pluscoloc::sensitivity()over a grid of priors. The eQTL side is a quantitative trait (N = per-gene eQTLGen sample size,sdY = 1), the CAD side is case-control (s = 60,801/184,305 = 0.33).Results
Reading of the results
Outputs:
out/coloc/coloc_report.md(tables, per-gene numbers, embedded figures),coloc_results.tsv(PP.H0-H4, lead SNPs, effect sizes, interpretation),coloc_sensitivity.tsv(PP.H4 over the p12 grid),coloc_region_<gene>.tsv(harmonised regional data) andplots/coloc_locus_<gene>.png|.pdfplusplots/coloc_pp_summary.png|.pdf.Protein and tissue replication: is the blood eQTL the right proxy?
Both steps above rest on one exposure, a whole-blood cis-eQTL. Two questions follow: does the protein behave the same way, and does vascular tissue reproduce the RNA signal strongly enough to re-run either analysis?
1. Plasma pQTL (public GWAS Catalog)
UKB-PPP (Synapse account) and deCODE (download request form) cannot be fetched unattended, but the GWAS Catalog hosts harmonised public pQTL summary statistics for both proteins, with complete regional statistics and no registration:
These files are GRCh38 + rsIDs and already carry beta, SE and the effect-allele frequency, so
scripts/09_fetch_pqtl.shonly has to cut the ±500 kb windows: there is no Z → beta step and no external allele-frequency file. The CAD GWAS is hg19, but every join is on rsID, never on position.Cis-MR, per 1 SD of plasma protein. The CXCL8 pQTL is strong (rs79656070, β = −0.134, P = 1.7e-18, F = 77) and points the same way as the blood eQTL:
NFKBIA cannot be tested at the protein level with public data: its only GWAS Catalog pQTL has n = 2,721 and its best cis-SNP reaches P = 9.8e-3, nowhere near the 5e-8 needed for an instrument. That is a data-availability limit, not a negative result.
Coloc, pQTL vs CAD — same code, same windows, same priors as above:
The CXCL8 pQTL is enormous (P = 1.7e-18) and still gives PP.H4 = 0.06, because the CAD side of the same window contains nothing significant to colocalise with (min CAD P = 2.7e-4). Moving from mRNA to protein shifts PP.H4 from 0.033 to 0.057: the ceiling is set by the outcome, not by the exposure. Protein-level data therefore reinforce the direction of the MR finding and leave the coloc null where it was.
2. GTEx v8 artery: can vascular tissue carry the analysis?
scripts/10_gtex_artery_check.Rreads the GTEx v8 single-tissue eGenes and significant variant-gene pairs already indata/mr/_gtex/(Artery_Aorta, Artery_Coronary, Artery_Tibial, Whole_Blood), recovers each tissue’s sample size fromminor_allele_count / (2·MAF), and looks every lead SNP up in the CAD GWAS:No artery instrument exists for either gene (0 of 6 gene–tissue pairs).
So for these two genes GTEx artery can support neither MR nor coloc; the full table is
out/tissue/gtex_tissue_eqtl.tsv. Coloc is blocked a second time by the data structure: the v8 archive distributes only the significant variant-gene pairs, andcoloc.abf()needs complete regional statistics. Where the tissue does have power the same code works: at the GTEx whole-blood CXCL8 lead SNP the Wald ratio is OR = 1.14 (P = 0.07 at n = 670), the same direction as eQTLGen (1.07) and the plasma pQTL (1.29).3. Why a blood eQTL can under-report a vascular target
Tissue specificity is a genuine limitation of the coloc step above and belongs in the interpretation of the PP.H4 values:
The GTEx check above shows that public artery data cannot currently settle the question — at n = 213–584 these two genes are not eGenes in coronary artery or aorta, so the missing artery-based replication reflects power and data availability, not absence of a vessel-wall effect. What would settle it: a well-powered artery eQTL (e.g. STARNET or GTEx v10 arteries), a plaque single-cell/pseudobulk eQTL, or the pQTL route above, which measures the circulating protein that a drug would actually engage.
Figures, manuscript and end-to-end reproduction
The figures in the manuscript are assembled from the tables in
out/, not from screenshots, so every panel is a vector PDF that can be regenerated:11_bulk_gse310013.Rdownloads the GSE310013 raw counts once intodata/bulk/(git-ignored) and runs a paired DESeq2 model (~ donor + day). With the default filter (>= 10 counts in >= 4 libraries, padj < 0.05, no LFC cut-off) this gives 11,650 DEGs, of which 117 fall in KEGG hsa05417, capturing 101 of the 103 scorable signature genes. An earlier, unarchived run of the same comparison reported 10,848 DEGs / 111 overlapping genes; the difference is the expressed-gene filter and the LFC cut-off, andout/bulk/venn_counts.tsvis the authoritative table for the current run. The 215-gene KEGG set is cached indata/atherosclerosis_kegg.txtso the figure does not depend on network access to Enrichr.12_wgcna_modules.Rrebuilds the module structure on the 5,048 most variable DEGs (plus all expressed signature genes): signed hybrid network, soft power 16, modules renumbered by |r| with culture day so that ME1 is the module most anti-correlated with Day 7 (r = -0.99, P = 9.3e-10). 101 signature genes fall in ME1/ME2 and 76 have kME > 0.8.20_ambient_sensitivity.Ris the reviewer-facing control for the cell-type ranking. It estimates the ambient profile from the 10% lowest-complexity cells, unmixes every cell at contamination fractions phi = 0-0.30 (c’ = max(0, c - phiTa), renormalised to the atlas scale T2 recovered per cell from X = log1p(c/T2*1e4)) and rescores all 259,116 cells. At phi = 0 the recomputedscore_meanreproducesout/genes__cell_type_level2/cell_scores.tsv.gzexactly, which is the built-in check that the implementation is correct. Outputs (ambient_profile.tsv,scores_by_celltype_by_phi.tsv,summary.tsv, the phi-sweep figure andambient_sensitivity.md) go toout/ambient_sensitivity/.13_assemble_figures.Rwrites:Composite figures are 180 mm wide (single column 90 mm, double column 180 mm), so they drop into a manuscript template without rescaling; each panel is also saved on its own for journals that want panels uploaded separately.
figure3/fig3Asubsamples the UMAP to ~70,000 cells (keeping every PLIN2+/TREM1+ macrophage) so the PDF stays a few MB;fig3Bbins the embedding instead of drawing 259k points, and the violin uses all cells.The manuscript is
manuscript/manuscript.md(Abstract, Introduction, Results, Discussion, Limitations, Methods, availability statements, references, legends). Numbers quoted in the text are the ones written by the scripts above; if a script is re-run with different thresholds, the text has to be updated with the newvenn_counts.tsv/module_trait.tsv/mr_results_main.tsvvalues.Caveats
Xis used as stored (log-normalised). No re-normalisation, no scaling, no imputation — this mirrors whatSeurat::AddModuleScore/PercentageFeatureSetwould see afterNormalizeData.var/feature_name; anything unmatched is written todata/genes_not_in_atlas.txt(alias problems: e.g. use current HGNC symbols).score_amcontrol bins come from a 60k-cell subsample, so values are stable to ~1e-3 but not bit-identical to a Seurat run on the full object.cell_type_level1/2), i.e. author annotations, not re-clustered here.cell_type(level 0) merges SMC + fibromyocyte.obs/disease,obs/dataset,obs/donor_idare available for stratified follow-up runs.Files