Raw RNAseq gene counts were obtained as noted above. After eliminating lowly expressed transcripts with <5 reads in more than 5 samples across the dataset, the count data was normalized using the RUVSeq (Risso et al., 2014) R Bioconductor package (Gentleman et al., 2004) to account for variance. As described in the RUVSeq manual, the normalization was accomplished in the following three-step procedure: 1) in silico negative control genes were determined by first-pass differential expression analysis by the edgeR (Robinson et al., 2010) and DESeq2 (Love et al., 2014) R Bioconductor packages, taking genes with FDR-adjusted P values >0.75, as calculated by both methods (approximately 7000 genes were unaffected by the condition of interest); 2) the negative control genes were then used in the RUVs function of the RUVSeq package, for calculation of variance factors; and, 3) the second-pass differential expression analysis (5% FDR and log2 fold change >1) for determining disease-dysregulated genes was performed using the original counts, adjusting for RUVs-calculated variance factors by multi-factor GLM models implemented in the edgeR and DESeq2 packages.
