The new pecotmr interface

[danielnachun]

I wanted to start a discussion here in a place that is more visible than Slack to clarify what remains to be done bring xqtl-protocol into alignment with pecotmr. There have been many changes that needed to be made for pecotmr to be compliant with Bioconductor as well as general bioinformatics engineering best practices, and that will require some changes to both xqtl-protocol and other custom users of pecotmr. I will first start by laying out the basic structure of pecotmr now, and then discuss places where changes are still needed.

Classes

A key requirement of Bioconductor is usage of S4 classes. These classes are powerful abstractions for complex data structures like the ones we encounter with pecotmr. I will not describe all the classes in pecotmr, just the ones that are relevant to us right now:

• GenotypeHandle: this is a file handle that points to genotype data on disk that is any of plink1 (.bed/.bim/.fam), plink2 (.pgen/.pvar/.psam), GDS or VCF. When genotypes are requested for downstream analysis, they are lazily loaded from those files. This handle can point to either a single file, or a set of files if the genotypes are split by chromosome.
• QtlDataset: this is the main container for individual-level data for a single study. Each context is stored as a SummarizedExperiment object with the actual molecular trait measurements as the main assay, the molecular PCs and other covariates as colData, and the gene positions as rowData. There also slots for a GenotypeHandle from which genotypes can be requested, and and genotypeCovariates slot to hold the genotype PCs. The most important functions related this class are:
  • getResidualizedPhenotypes(context, traitId/region): this returns a list of phenotype matrices adjusted for (by default) genotype and phenotype covariates. The list is named by context, and multiple traits can be requested simultaneously (to support multi-trait analysis with mr.mash, mvsusie and fsusie. The multi-trait selection can be with either a list of traits, or by instead providing a GRanges range.
  • getResidualizedGenotypes(context, traitId/region, cisWindows): this is the same idea as getResidualizedPhenotypes but for genotypes. The genomic range to extract is either directly taken from a region if given, or by getting the positions of traitIds from the rowData and then adding the cis-window.
• QtlSumStats/GwasSumStats: these classes hold Granges objects with summary statistics for QTL and GWAS, respectively. For QTL, we keep the summary statistics separate by study, context, and trait, while for GWAS, we only need to track study. The only other large structure in these classes is a GenotypeHandle for the LD sketch.
• MultiStudyQtlDataset: this classes has slots for a QtlSumStats object and a list of QtlDataset objects, one per study. The requirement here is that the instantiation of this class requires either at least two QtlDataset objects, or at least one QtlDataset object and one QtlSumStats object. You can of course have both a list of QtlDataset objects and a QtlSumStats object too.
• QtlFineMappingResult/GwasFineMappingResult: These classes hold the results of a fine-mapping analysis with fineMappingPipeline. This is includes the SuSiE fits as well as the topLoci slot which can be returned as a data frame or GRanges object.
• TwasWeights: This class holds the results of a twasWeightsPipeline analysis, which are mainly just regression fits from which TWAS weights can be extracted.

Pipeline functions

• sumStatsQc: This is the sole home of all of our summary statistics QC now. It takes a QtlSumStats or GwasSumStats object and returns another with QC applied to it.
• mashPipeline: This pipeline takes either a QtlSumStats object or a QtlFineMappingResult and uses mashr to estimate the mash prior for input to mr.mash and mvsusie. It returns an object of type MashPrior.
• fineMappingPipeline: This pipeline takes in any of a QtlDataset, QtlSumStats, GwasSumStats or MultiStudyQtlDataset object, and returns a FineMappingResult after running any of susie, mvsusie or fsusie. When a MultiStudyQtlDataset object is given, the fine-mapping is run separately per QtlDataset object (since there is no "multi-study" individual-level fine-mapping right now). Multi-trait and multi-context fine-mapping are handled automatically when mvsusie and fsusie are selected as methods. Multi-study fine-mapping can be done with mvSuSIE RSS on a QtlSumStats/GwasSumStats) object (including if the the QtlSumStats object is inside of a MultiStudyQtlDataset object).
• twasWeightsPipeline: This pipeline takes in any of a QtlDataset, QtlSumStats, or MultiStudyQtlDataset, and optionally a FineMappingResult object and learns TWAS weights with regularized regression. Multi-trait and multi-context regression is handled automatically when mr.mash is requested. MultiStudyQtlDataset objects are handled the same way as in fineMappingPipeline, with multi-study regularized regression being possible only with mr.mash RSS. The ensemble method is run automatically if more than regression model is requested. If a FineMappingResult is provided, the weights can be extracted from this object and added to the ensemble model.
• qtlEnrichmentPipeline: This pipeline takes a QtlFineMappingResult and a GwasFineMappingResult and computes heritability enrichment of the GWAS for QTLs as originally implemented by fastenloc.
• colocboostPipeline: This pipeline takes in any of a QtlDataset, QtlSumStats, GwasSumStats or MultiStudyQtlDataset and runs colocboost. Since colocboost is able to directly handle all of these types, the main thing this pipeline does is just extract all the right inputs and passes them to colocboost. We may still need to discuss options for making the return value of this nicer.
• causalInferencePipeline: This pipeline takes at least one of a TwasWeights and QtlFineMappingResult object and a GwasSumStats object and computes TWAS Z-scores, and if a QtlFineMappingResult object was provided, MR Wald ratios. The return value is just a data frame, though we can iterate on this more if needed.
• ctwasPipeline: This pipeline takes a TWAS Z-scores and a GwasSumStats object and runs CTWAS on them. If multiple TWAS Z-scores exist for the same context and trait (this means multiple methods were used), it can iterate over these. And the user can request either marginal CTWAS per-context or multi-context CTWAS.
• colocPipeline: This pipeline takes at least one of a TwasWeights and QtlFineMappingResult object and either a GwasSumStats object or a GwasFineMappingResult object. When a GwasSumStats object is provided, window-level fine-mapping of the GWAS is run, while when a GwasFineMappingResult is provided, window-level fine-mapping is skipped. If the enrichment information from qtlEnrichmentPipeline is provided, then the pipeline is running SuSiE-enloc while when it is absent, it is just running the regular SuSiE-coloc.

Pipeline ordering - shared

This is the key logic everyone needs to understand moving forward. The above class and pipeline definitions are to make this section interpretable:

1. Construct the QtlDataset, QtlSumStats, MultiStudyDataset and GwasSumStats objects.
   • QtlDataset construction should not be parallelized or split by gene/region, because our classes now know how to extract specific columns from the phenotype data.
   • QtlSumStats and GwasSumStats objects do need to be constructed per gene/region because loading a full GWAS summary statistic would use a lot of memory.
   • MultiStudyQtlDataset objects that only contain QtlDataset objects should also be constructed once without splitting over gene/region. If they need to also contain QtlSumStats, those should be made per gene and added afterward.
2. Apply sumStatsQc to QtlSumStats (if used) and GwasSumStats objects. For now, if QtlSumStats need to be added to a MultiStudyQtlDataset, the QC should be done before adding the object to the MultiStudyQtlDataset.

Pipeline ordering - TWAS/MR

1. (mr.mash/mvsusie only and optional) Run the univariate fineMappingPipeline per-trait/per-context on the QtlDataset/QtlSumStats/MultiStudyQtlDataset object to get a FineMappingResult to pass to the mashPipeline. This can also double as a screening step to decide where muti-trait/multi-context fine-mapping should be performed
2. (mr.mash/mvsusie only) Run the mashPipeline to get a MashPrior. Step 3 is needed if the user wants to use a FineMappingResult to make it easier to select strong signals.
3. (mvsusie only) Run the twasWeightsPipeline for mr.mash with the QtlDataset/QtlSumStats/MultiStudyQtlDataset to get a TwasWeights object with a mr.mash fit with mixture weights and residual variance to pass to mvsusie.
4. Run the main fineMappingPipeline call on the QtlDataset/QtlSumStats/MultiStudyQtlDataset to get a FineMappingResult. If the analysis is univariate only, then this is the only time fineMappingPipeline needs to be called. Otherwise if step 1 was already run, the FineMappingResult from there can be provided here and only new methods (multivariate or additional univariate) are run. mvsusie requires the output of steps 2 and 3 as additional inputs.
5. Run the main twasWeightsPipeline call on the QtlDataset/QtlSumStats/MultiStudyQtlDataset to get a TwasWeights object. If a FineMappingResult is provided, its weights are extracted and used in the ensemble model without rerunning fine-mapping. For the multivariate case, mr.mash and mvsusie will have already been run, so the only thing left to do there is the ensemble model between those two methods.
6. Run the causalInferencePipeline on the FineMappingResult from step 4 (can get TWAS Z-scores and MR) and/or the TwasWeights output from step 5 (TWAS Z-scores only).

Pipeline ordering - genome-wide GWAS fine-mapping

1. Run the fineMappingPipeline on the GwasSumStats object from step 2 of the shared pipeline per LD block.

Pipeline ordering - univariate colocalization

1. (SuSiE-enloc only) - run the qtlEnrichmentPipeline on the QtlFineMappingResult from step 2/step 4 of the TWAS/MR pipeline and the GwasFineMappingResult of step 1 of the genome-wide GWAS fine-mapping pipeline.
2. Run the colocPipeline, optionally with the enrichment from step 1, on the QtlFineMappingResult from step 2/step 4 of the TWAS/MR pipeline and either the GwasSumStats object from step 2 of the shared pipeline or the GwasFineMappingResult from step 1 of the genome-wide GWAS fine-mapping pipeline.

Pipeline ordering - colocboost

1. Run colocBoostPipeline on one of: the QtlDataset object from step 1 of the shared pipeline, the QtlSumStats object from step 2 of the shared pipeline, or the MultiStudyQtlDataset object from step 1 or 2 of the shared pipeline and the GwasSumStats object from step 2 of the shared pipeline.

Key improvements from before

There are several critical improvements that have made or will make the protocol more intuitive and make unit-testing and development much easier:

• Front loading data management: The idea here is that main source of friction is getting the data assembled into the right shape. By handling all of this during the construction of QtlDataset/QtlSumStats/MultiStudyQtlDataset and forcing users to strictly organize their data is a sensible way, we greatly reduce failures in later steps due to data management issues
• Modularized I/O: Related to my first point, only the first step of constructing the QtlDataset/QtlSumStats/MultiStudyQtlDataset objects requires dealing with messy I/O. Once that is done, all subsequent pipeline steps are only reading RDS objects or raw genotypes.

Future improvements

There is room for improvement here. My current ideas:

• loadQtlDatasetFromManifest/loadQtlSumStatsFromManifest/loadGwasSumStatsFromManifest/loadMultiStudyQtlDatasetFromManifest functions: These functions would take in a manifest file with the right columns and call the constructors for the user. A couple of restrictions:
  • For individual level data, the user must either provide a tensorQTL-style BED file with molecular trait positions and the quantifications, or a matrix with only quantifications and a separate GTF or tabular file with molecular trait positions.
  • For sumstats, we can support loading from a GRange only if the underlying files are tabix indexed or are GWAS VCFs.

This is a huge document so please ask as many questions as needed to clarify things!