## 4 Materials and Methods
### Statistics  
Statistical analyses were performed using R-4.1.3 ([[http://www.r-project.org](http://www.r-project.org)]([http://www.r-project.org](http://www.r-project.org))) [103] with the additional packages beeswarm_0.6.0 [104], tidyverse_1.3.1 [105], qqman_0.1.4 [106] and SNPassoc_2.0-11 [107]. Associations between demographic and clinical factors and _SLC10A1_ mRNA expression and protein levels were investigated using Wilcoxon-Mann-Whitney tests or [PERSON]'s rank correlation coefficients as appropriate. Univariate and multivariate association analyses between the genotyped variants and _SLC10A1_ mRNA expression or protein levels were performed using the generalized linear model capabilities of the SNPassoc package. For this purpose, _SLC10A1_ mRNA and NTCP protein measurements were first power-transformed (_SLC10A1_: \(\lambda\) = 0.33, NTCP: \(\lambda\) = 0.7) to fulfill Gaussian distribution assumption. The transformation was determined via the R-package car_3.0-2 [108]. Distributions of transformed measurements were confirmed using a quantile-quantile plot and Shapiro-Wilk tests. In multivariate analysis, we used seven demographic and clinical factors (gender, age, CRP, cholestasis, medication, alcohol, and nicotine consumption) as covariates. In addition, associations between the _SLC10A1_ genetic variants and the bile acid concentrations were analyzed using Kruskal-Wallis tests. Association study of genome-wide imputed genotypes was performed with SNPTEST v2.5 using the frequentist association test option with expected genotype counts. Imputation of the GSA-SNP array data was performed using Minimac3 software on the Michigan Imputation Server [109]. File preparation was performed using PLINK software and HRC preparation checking tool version 4.3.0 ([[https://www.well.ox.ac.uk/~wrayner/tools/](https://www.well.ox.ac.uk/~wrayner/tools/)]([https://www.well.ox.ac.uk/~wrayner/tools/](https://www.well.ox.ac.uk/~wrayner/tools/))). Vcf-files of autosomal SNPs were subject to imputation. HRC-Tool, eagle_v2.3 for phasing and Minimac3 were run using default parameters and hrcr.1.2016 haplotypes as reference panel. The association between each of the genetic variants and the transformed measurements was analyzed in the additive genetic model. For the analyses of linkage disequilibrium and haplotypes the Haploview version 4.2 software ([[http://www.broadinstitute.org/haploview](http://www.broadinstitute.org/haploview)]([http://www.broadinstitute.org/haploview](http://www.broadinstitute.org/haploview)), Cambridge, MA, USA) was used. Hardy-Weinberg equilibrium calculations were performed to compare observed and expected allele and genotype frequencies. Regulatory motifs affected by genetic variants were identified using HaploReg version 4.1 ([[http://www.broadinstitute.org/mammals/haploreg/haploreg.php](http://www.broadinstitute.org/mammals/haploreg/haploreg.php)]([http://www.broadinstitute.org/mammals/haploreg/haploreg.php](http://www.broadinstitute.org/mammals/haploreg/haploreg.php))). Correlation coefficients between mRNA and protein levels and endogenous metabolites were determined using Spearman's rho (r\({}_{\text{S}}\)). Associations between SNPs (MAF > 1%) and metabolites were investigated using Kruskal-Wallis test for univariate analysis as well as ANOVA for median regression for multivariate analysis. To be more precise, for each metabolite phenotype and each chosen SNP, we applied function anova.rq in R-package quantreg-4.91 with rank test-statistic and Wilcoxon score function to compare two median regression fits with only the seven non-genetic demographic and clinical factors as covariates versus a model of the seven non-genetic factors plus the considered SNP. All statistical tests were two-sided and statistical significance was defined as \(p<0.05\). Where indicated, adjustment for multiple testing was carried out using the Benjamini-Hochberg procedure [110].