Stem Cell

Effects of psoriasis and psoralen exposure on the somatic mutation landscape of the skin – Nature Genetics

Time Stamp: October 25, 2023 8:00 PM
Source Node: 214595

This study complies with all relevant ethical regulations. All donors gave informed consent for genetic research of the material and the study was approved by the research ethics committee of Christian-Albrechts University in Kiel (A100/12), the National Health Service (NHS) Research Ethics Committee (Yorkshire and The Humber—South Yorkshire Research Ethics Committee, REC ID 20/YH/0244, IRAS ID 286843) and by the Wellcome Trust Sanger Institute Human Materials and Data Management Committee (approval number 20/0085).

Statistics and reproducibility

This is a descriptive study with no interventions. No randomization of patients was applied. All dependent variables were generated computationally and no blinding was applied. Sample size was not predetermined. We excluded from the whole study six microbiopsies suspected of being sample swaps and/or contaminated with external DNA. A total of 101 microbiopsies were excluded only from the analysis of structural variants either due to the inability of the algorithm to find an optimal solution or if the goodness-of-fit for the optimal solution was <90%.

Human tissue attainment and processing

Tissue donation and fixation

Round punch biopsies, 4 mm in diameter, from lesional and adjacent nonlesional skin were donated by psoriasis patients presenting to the Department of Dermatology, University Hospital Schleswig-Holstein, Kiel between 2017 and 2019. Biopsies were fixed in RNAlater (ThermoFisher, catalog no. AM7021) upon collection following the manufacturer’s instructions and frozen. Patients were not compensated for participation in the study.

LCM of epidermis

Skin biopsies were processed using a Tissue Tek VIP 6 AI tissue processor (Sakura Finetek) and embedded in paraffin. The blocks were sectioned into 10-μm-thick sections using an Accu-Cut SRM 200 microtome (Sakura Finetek). Sections were fixed to 4 μm PEN membrane slides (Leica, catalog no. 11600288). Slides were stained with hematoxylin and eosin and imaged using a NanoZoomer 2.0-HT slide scanner (Hamamatsu Photonics). Samples of epidermis were dissected from this material using LCM microscopy (Leica, catalog no. LMD7000). Cells were lysed using Arcturus PicoPure DNA extraction kit (Applied Biosystems, catalog no. KIT0103) according to the manufacturer’s instructions.

DNA sequencing

Whole-exome sequencing

A total of 1,182 microbiopsies (Supplementary Table 2) from 111 individuals were whole-exome sequenced on Illumina NovaSeq 6000 machines using 150 bp paired-end reads and the Agilent SureSelect Human All Exon V.5 bait set (catalog no. S04380110). Paired-end reads were aligned to the human reference genome (build hg38) using BWA-MEM⁵⁸. PCR duplicates were marked using biobambam⁵⁹ and duplicate statistics were calculated using Picard (v.1.131) (http://broadinstitute.github.io/picard/). Sample contamination estimates were calculated using VerifyBamID (v.1.1.3)⁶⁰ and the on-target coverage was calculated using Samtools (v.1.11) depth command, considering only reads with base quality and mapping quality >30. The median on-target coverage for all samples was 56× (range, 18–190) for the dataset as a whole but 58× and 46× for microbiopsies from lesional and nonlesional skin, respectively. The difference in coverage is due to the greater volume of the microbiopsies from lesional skin (Extended Data Fig. 1), which results in greater library complexity and translates to lower PCR duplicate rate and higher coverage. We cut columns of cells and the thicker epidermis of lesional skin results in a greater sample volume for the same surface area covered.

Whole-genome sequencing

WGS was performed on 16 microbiopsies from three donors (patients 18, 21 and 34) who showed the signature of PUVA exposure in the WES data. The sequencing of these samples is further described in Supplementary Note 1.

Somatic mutation calling

Mutation calling in whole exomes

Substitutions were called using CaVEMan (v.1.15.1) (Cancer Variants through Expectation Maximization) (https://cancerit.github.io/CaVEMan/)⁶¹. Mutations were called against an unmatched normal with the copy number options set manually to 10 and 2 for the mutant and wild-type copy numbers, respectively. The samples were compared against a normal panel consisting of 75 unrelated normal samples to remove common SNPs. Mutations were further filtered if the reads reporting the mutations had a median alignment score <140 or if >50% of the reads were clipped.

All mutations passing these filters in any sample from a donor were next genotyped in all samples from that donor. We used the bam2R() function of the deepSNV⁶² R package (v.1.40.0) to generate pileups of all sites mutated in any sample using only mapped reads that had base quality and mapping quality >30 and which were mapped in a proper pair, were not PCR duplicates, were the primary alignment and which passed platform quality check (SAM Flags 3847, see https://broadinstitute.github.io/picard/explain-flags.html). After these filters, we required a coverage of at least four times at the site and at least three reads reporting the alternate allele in at least one sample from a donor to call a mutation.

Adjacent substitutions called in the same sample were merged into a double-base substitution call if the number of reads reporting the reference and alternative alleles was not significantly different (Fisher test).

Indels were called using cgpPindel (v.3.5.0) (https://github.com/cancerit/cgpPindel)⁶³ using the same unmatched normal sample that was used to call substitutions. We generated pileups of the indel calls in the same way as described above for substitutions and required a coverage of at least four times and at least three reads reporting the alternate allele in at least one sample from a donor to call a mutation as before.

Binomial filtering of somatic mutation calls

To filter rare germline variants not removed by the comparison with the normal panel we applied an exact binomial test of the number of reads reporting each mutation²⁶. Heterozygous germline variants are expected to be present at a VAF of 0.5 in every sample from a patient. For each mutation, we compared the number of reads reporting the reference and alternate alleles across all samples from that patient. We tested the hypothesis that the read counts for the variants were drawn from a binomial distribution with a probability of success of 0.5, or 0.95 for mutations on the sex chromosomes in men. We applied Benjamini–Hochberg correction for multiple testing and excluded mutations with q > 10^–3. We also used binomial filtering to remove erroneous mutation calls. Recurrent sequencing artefacts will be distributed randomly across samples and can be modeled as being drawn from a binomial distribution. In contrast, true somatic mutations will have a high VAF in some samples whilst being completely absent from others. The latter are best represented by a beta-binomial with a high overdispersion. For every mutation call, we calculated the maximum likelihood overdispersion parameter (ρ) in a grid-based way (ranging the value of ρ from 10^–6 to 10^–0.05)²⁶. Calls with ρ < 0.1 were filtered as probably artifactual.

The sensitivity of the mutation calling postfiltering was estimated at 89% as described in the Supplementary Note 1.

Structural variant calling in the exome data with allele-specific copy number analysis of tumors

Structural variants were called using allele-specific copy number analysis of tumors⁶⁴(https://github.com/VanLoo-lab/ascat) as further described in Supplementary Notes 1 and 2.

Identification of SNV clusters by a hierarchical Dirichlet process

We implemented a nonparametric Bayesian hierarchical Dirichlet process (HDP) to cluster autosomal SBSs with similar VAFs. The full mathematical and implementation details of the model are described in a previous publication⁸. Briefly, clones of cells are present across different microbiopsies and this manifests as clusters of mutations that are found at similar VAFs within a microbiopsy. For every mutation, we have two vectors, one containing the number of reads reporting the alternate allele in each microbiopsy and another containing the total sequencing depth at each microbiopsy. We assume that each mutation can be assigned to exactly one cluster but the number of clusters is unknown. We aim to estimate the number of clusters present across all the microbiopsies dissected from a patient, the location of each cluster in the n-dimensional VAF hypercube and the allocation of mutations to each cluster.

We model the data using an N-dimensional Dirichlet process clustering model, where the distribution of clone sizes and numbers follows a Dirichlet process. This has the advantage that there is no need to prespecify the number of clusters present. Instead, mutations are moved around the clusters and, in each sampling iteration, there is a defined probability that a mutation will initiate a new cluster that was not present in previous iterations. Clusters can also cease to exist if all member mutations are assigned to other clusters. Thus, the number of clusters varies throughout the sampling chain. Details of the implementation are provided in the Supplementary Note 1 and in the code accompanying the manuscript.

Inference of phylogenetic trees

Each cluster of single-base-substitutions identified by the N-dimensional Dirichlet process algorithm represents a branch of the phylogenetic tree for that patient. We applied the statistical piegonhole principle to infer phylogenetic relationships between clusters. Given clusters A and B, if the combined mutant cell fraction (CF) of both is >100% (VAF > 0.5) in the same microdissections and B consistently shows a lower CF than A, then that is strong evidence that B is nested within A, that is, mutation cluster B represents a subclone of clone A. If the combined mutant CF is ≤100%, only weak evidence of nesting exists. If B is found at a higher VAF than A in some microdissections but at lower VAF in others, the clusters are interpreted as being independent clones without nesting. We treated each tip of the phylogenetic tree for each patient as a clone. The length of the branches from the root (germline) was used in the mutation burden calculations.

Some clusters are present at VAFs too low for the pigeonhole principle to be incontrovertible (that is, combined cell fraction <100%). In these cases, there is a risk that the cluster does not represent a single clone as assumed above, but a mixture of clones present at similar VAFs. We reconstructed the phylogenetic trees after pruning away branches where there is doubt about the validity of the pigeonhole principle. We retained nested clusters only if the sum of the cellular fraction estimates exceeds 1. Un-nested clusters (that is, branches of the phylogenetic tree consisting of a single cluster) were retained if the median VAF of the cluster is >0.3.

Mutational signatures

Extraction of mutational signatures

To extract mutational signatures and estimate the exposure of each signature, we used a second hierarchical Dirichlet process⁶⁵, as implemented in the HDP (v.0.1.5) R package (https://github.com/nicolaroberts/hdp). The data were organized into a tree structure where the root contained all the mutations in the dataset. This node had, as children, one node that represented the most recent common ancestor of all the patients and frozen pseudocount nodes for signatures that are to be used as priors in the model. The pseudonodes contained 10,000 pseudocounts each. We used COSMIC signatures 1, 5, 2, 13, 7a, 7b, 7c, 7d, 17a, 17b, 18 and 38 as priors. During the Dirichlet process, mutations from the dataset may join the pseudocount clusters, but the pseudocounts are frozen such that they are unable to leave the initial cluster. The patient ancestor node had, as children one node for each of the patients, and each patient node had as children one node for each branch of the phylogenetic tree. The parameters used are provided in Supplementary Note 1 and the code accompanying the manuscript.

The model extracted nine signature components in addition to the unassigned component (Extended Data Figs. 2 and 3). Among these were components corresponding to signatures SBS7b, SBS1/5, SBS2, SBS7c and SBS13, all of which were included as priors in the model. One component did not correspond to any COSMIC signature but is characterized by mutations at TpA sites and accounted for just under 11% of the mutations in the dataset. This is the signature we attribute to psoralen exposure, which, in the context of psoriasis, is likely to occur during treatment with psoralens and high-dose UV-A (PUVA) (see the main text). Finally, three additional components were extracted: unknown components N1–3 in Extended Data Fig. 2. These together accounted for 2.6% of the mutations in the dataset. They may represent individual variation in repair of UV-damage or they may be artefacts of the signature extraction model. We do not have sufficient confidence in these components to draw conclusions from them and have added them to the unassigned component for subsequent analyses.

Characterization of the psoralen signature

From the whole-exome data, we identified a number of samples that showed a large number of mutations at TpA sites, consistent with the known mutagenic effects of psoralens^31,32,33. To enable further characterization of this signature, we selected 16 microbiopsies from patients showing clear evidence of psoralen exposure for WGS.

To visualize the trinucleotide and pentanucleotide spectrums associated with psoralen exposure, we used the R package MutationalPatterns⁶⁶ (v.3.4.0) together with BSgenome (v.1.60.0). To calculate the transcriptional strand bias, we used the gene definitions from the R package TxDb.Hsapiens.UCSC.hg38.knownGene (v.3.13.0) and the strand_occurrences() and strand_bias_test() functions from MutationalPatterns.

To test for TCD, we carried out a similar analysis to that used originally to describe TCD in liver cancers⁶⁷. We divided protein-coding genes into quintiles by ascending expression in sun-exposed skin from the GTEx dataset⁶⁸ (v.8). We extracted the transcriptional start site (TSS) and the strand of each gene from Gencode (v.27) and defined ten 1 kb bins upstream and downstream of the TSSs. We pooled T > (ACG) and A > (CGT) mutations at TpA or ApT sites from all whole-genome sequenced samples. If the gene is on the (–) strand, the transcribed strand is the reference and we counted the number of T > (ACG) mutations overlapping each 1 kb bin. If the gene is on the (+) strand, the transcribed strand is the complement of the reference and we counted the number of A > (CGT) mutations. This was reversed for the untranscribed strand. We observed a drop in the mutation rate on the transcribed strand upstream of the TSSs, indicating TCR. However, we also found an increased mutation burden on the untranscribed strand, indicative of TCD of this strand. To test the statistical significance of the increased mutation burden, we fit two linear models with and without a parameter indicating whether each position was upstream or downstream of the TSSs and used a likelihood ratio test to test whether the fit of the model was improved. Figure 3c shows the mutation rate in each 1 kb bin relative to the –10 kb bin—the intergenic bin furthest from the TSS.

To test the effect of gene expression levels on PUVA mutagenesis we again used expression data from sun-exposed skin from the GTEx dataset⁶⁸ (v.8). We split protein-coding genes into ten equally sized bins by ascending levels of expression. We used the Bedtools⁶⁹ intersect function (v.2.18) to count the number of mutations overlapping genes in each bin. Figure 3d shows the relative mutation rate in each bin compared with the lowest expression bin. We also looked for an effect of the replication timing and replication strand as described in Supplementary Notes 1 and 5.

To assess the potential functional effects of the psoralen signature relative to other main signatures in the skin (Fig. 3f), we used the context_potential_damage_analysis() and signature_potential_damage_analysis() functions of the MutationalPatterns package in R. Considering only the genes previously reported to be under selection in squamous cell carcinomas and/or normal skin, we counted the number of trinucleotide changes expected to give rise to the different types of mutations (synonymous, missense, nonsense or splice site). For each of the main signatures in the skin (Psoralen, SBS1, SBS5, SBS7a and SBS7b), we normalized the expected fraction of mutations falling in each annotation class by the fractions that would be expected given either a uniform mutation rate or SBS7b. We note that, whereas this analysis gives a hint of how damaging a signature might be, it is based only on a trinucleotide mutational context and other important features, such as gene expression, strandedness or the extended nucleotide context, are not taken into account. The genes we used in this analysis are AJUBA, ARID2, ASXL1, CASP8, CDKN2A, FAT1, KMT2D, NOTCH1, NOTCH2, NOTCH3, PPM1D, RB1, RBM10, TP53 and TP63.

To see whether psoralen exposure was associated with increased levels of clonal spread after controlling for age, we estimated the median VAF for microbiopsies dissected from each patient and regressed this against the age of the patient in a linear model. We then fitted a second linear model which also included as covariate whether the psoralen signature was seen in that patient. We used a likelihood ratio test to see whether this improved the fit of the model (Fig. 3f).

Mutation burden estimation

Linear mixed effect models

We used linear mixed effect models to compare the mutation burdens of lesional and nonlesional skin and test for an effect of disease duration. We used as response variables the estimates of the mutation burdens of SBS7 and SBS1/5 after correcting for VAF and coverage as described in Supplementary Notes 1 and 3. The models include fixed effects for age and the anatomical location of the biopsy from which the clone is derived. Cell clones from the same biopsy are likely to have correlated levels of UV-exposure and some correlation is also likely to exist between biopsies taken from the same patient. To model this, we include random effects for patient and biopsy, with the effect of biopsy being nested within that of the patient. We provide a mathematical description of the models used in Supplementary Note 1 and implement the models in R in Supplementary Note 4.

Selection and driver analyses

Exome-wide driver discovery

We used the dNdScv (v.0.0.1.0) software⁴³ (https://github.com/im3sanger/dndscv) to identify genes enriched in nonsynonymous mutations, indicative of positive selection. We first calculated dN/dS ratios across all coding genes using mutations identified in either lesional or nonlesional skin, identifying positive selection of mutations in nine genes (Fig. 4). We next estimated global dN/dS values for groups of genes belonging to 11 genesets defined a priori as described in Supplementary Note 1.

Fraction of mutated cells

We compared the fraction of cells that carry mutations in any of the nine genes that showed a significant enrichment of mutations between lesional and nonlesional skin. We calculated the fraction of mutated cells separately in lesional and nonlesional biopsies from each individual by multiplying twice the VAF of each mutation in each microbiopsy by the volume of the microbiopsy and dividing that by the total volume of microbiopsies dissected from that patient. For clones that carry more than one mutation in the same gene, we counted only the mutation with higher VAF.