Alternative Splice Sites

For the simulation we used the quantification of RefSeq transcripts for the three control samples from [17 (link)] (GSE59335) with Salmon [31 ] as theoretical abundances, and considered genes with only two isoforms containing a skipping exon (SE) or alternative splice site (A5/A3) event and only one associated event. For the benchmarking analysis, we selected a set of positive and a set of negative events for each event type with the same number of randomly chosen events, 277 for SE events and 318 for A5/A3 events. For the positive set we simulated differential splicing by exchanging the theoretical abundance of their associated transcript values. We selected to be positive events only those having an absolute difference of relative abundance greater than 0.2, so that the simulated change was sufficiently large: $$\frac{\mid {\mathit{TPM}}_1-{\mathit{TPM}}_2\mid }{{\mathit{TPM}}_1+{\mathit{TPM}}_2}>0.2$$
where TPM1 and TPM2 are the abundances for the two transcripts in the gene, given in TPM units. For the negative set, we took an equal number of events without exchanging their TPM values. These negative events had a gene expression distribution and a distribution of transcript relative abundance similar to the positive events, and an expected variability between conditions similar to the variability between biological replicates. We used RSEM [39 (link)] to simulate sequencing reads for the two conditions, three replicates each, at various depths (120, 60, 25, 10 and 5 M 100-nt paired-end reads per sample) and at various read lengths (100, 75, 50, and 25 nt, at a depth of 25 M paired-end reads) (Additional file 2: Tables S1–S3). Further details of the simulations are given in the Additional file 3:Supplementary material. Datasets and commands to reproduce these simulations are available at https://github.com/comprna/SUPPA_supplementary_data.

To detect alternative mRNA associations with AUD we used Leafcutter version 0.2.9^{23 (link)}. Leafcutter is a powerful transcriptome-wide splicing method that uses a Dirichlet-multinomial generalized linear regression to identify differentially spliced genes. A differentially spliced gene generally is composed of multiple clusters, each of which includes various alternative splicing events, such as exon-skipping (see Fig. 1), intron retention, alternative acceptor or alternative donor splice sites, which we annotated with the Vertebrate Alternative Splicing and Transcription Database (https://vastdb.crg.eu/wiki/Main_Page). Each splicing event corresponds to a change in percent spliced in (ΔPSI or dPSI) metric. In our AUD analyses, a positive ΔPSI for an exon skipping event would suggest that an individual with AUD is more likely to skip a certain exon than someone without AUD. We utilized the default filtering parameters of Leafcutter that filtered out splicing clusters with < 5 samplers per intron, < 3 samples per group, and required at least 20 reads, which resulted in 18,685 unique genes across human brain regions. Human differential splicing analyses covaried for sex, age, brain pH, PMI, and smoking status. Note leafcutter performs analyses at the cluster level calculating a cluster p-value and then performs a Benjamini–Hochberg False Discovery (BH-FDR) multiple testing correction. Differentially spliced genes/clusters were those that survived a standard BH-FDR adjusted p-value < 0.05. We corrected p-values for multiple testing within brain regions and thus, our analyses do not account for multiple testing across tissues or samples. Since only 21 genes were differentially spliced in primates (BH-FDR < 0.05), we defined significant differential splicing with a nominal p-value threshold < 0.05. When possible, primate differential splicing analyses controlled for age (NAc sample). We assessed linear correlations of the ΔPSI across all significant alternative splicing events that were common across brain regions.
To assess the overlap between human and primate results we used a Fisher’s Exact test at the gene-level and restricted analyses to homologous genes identified by biomaRt^{24 (link)} and only used results from analogous regions of the brain (CEA, NAc, and PFC). In humans, we compared our differential splicing analyses with differentially expressed genes. Differential expression analyses leveraged featureCounts to count aligned RNA-seq reads and used DESeq2^{25 (link)} to determine differential expression. Differential expression analyses used the same covariates and p-value adjustment as differential splicing analyses. Previous differential splicing analyses of these data^{7 (link)} used rMATS^{26 (link)} that focuses on individual splicing events (rather than broader clusters within genes) and leverages a joint likelihood function combining binomial and normal distributions.

We used rMATS (https://rnaseq-mats.sourceforge.net/index.html) to detect and analyze alternative splicing events, including skipped exon, alternative 5′ splice site (A5SS), alternative 3′ splice site (A3SS), mutually exclusive exons (MXE), and retained intron (RI).