RNA-seq data analysis

RNA-seq reads were aligned to the mm10 genome using Rsubread 1.35.19 align and using Rsubread’s inbuilt mm10 RefSeq gene annotation to aid the identification of exon junctions (59). Read counts were obtained for Entrez Gene IDs using featureCounts and Rsubread’s inbuild annotation (60). Gene annotation was downloaded from ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/GENE_INFO (March 2020). Differential expression analyses were undertaken using the edgeR v3.28.0 and limma v3.42.0 software packages (61, 62). Unexpressed genes were filtered using edgeR’s filterByExpr function with min.count = 20 and min.total.count = 25. Mitochondrial genes, ribosomal RNA, genes on the X or Y chromosomes, genes of type “other,” and obsolete gene IDs were also filtered. Library sizes were normalized by edgeR’s trimmed mean of M-values (TMM) method.

Differential expression was assessed using the limma-trend method with arrayWeights and a duplicate correlation to account for correlation between mice (63). Gene counts were transformed to log₂ counts per million. Duplicate correlation was calculated with the individual mouse as the block using the duplicateCorrelation function from limma. The array weights per sample were estimated with the arrayWeights function from limma with method = “genebygene.” Robust empirical Bayes was used to protect against hypervariable genes (64). Differential expression was assessed using t tests relative to a threshold (TREAT)–moderated t tests (64) relative to a fold change threshold of 1.2. Genes with a TREAT false discovery rate below 0.05 were considered to be differentially expressed.

Gene set enrichment was tested using limma’s fry function, which provides fast approximation to rotation testing (66). Gene set enrichment was visualized using limma’s barcode enrichment plot. Heatmaps of the filtered and normalized logCPM value were plotted with the coolmap function from the limma package. Mean difference plots (MD plots) were plotted with the plotMD function.