RNA-seq data analysis.

Sequences were mapped to the mm10 genome by using the STAR aligner (Dobin et al., 2013). Transcript-level data was counted using HTSEQ (Anders et al., 2015). Differential gene expression was modeled using voom (Law et al., 2014), available in the limma R software package. Normalization factors were generated using the TMM method. The voom-normalized counts were analyzed using the lmFit and eBayes functions in limma. The adjusted P value was estimated using the Benjamini-Hochberg method. Gene set enrichment analysis (GSEA) was performed using gene sets obtained from MSigDB (the Broad Institute) as described previously (Subramanian et al., 2005). The ranking of genes was calculated using the signal-to-noise algorithm, and a P-value for each gene set was estimated by comparing the observed enrichment score to that obtained from a null distribution computed from 1000 permutations of genes within the gene sets. The false-discovery rate (FDR) was estimated as described previously (Subramanian et al., 2005).