Bulk RNA-seq library preparation and analysis

OT-I CD8 T cells transduced with MigRI (n = 3)/Klf4 (n = 3) were repeatedly stimulated with OVA peptide in vitro and then harvested. Total RNA was extracted using TRI Reagent according to the manufacturer’s instructions (Molecular Research Center Inc.). cDNA libraries were prepared using the TruSeq stranded mRNA Sample Preparation Kit (Illumina, CA, USA). cDNA libraries were quantified with the KAPA library quantification kit (Kapa Biosystems, MA, USA) according to the manufacturer’s library quantification protocol. Agilent 2100 BioAnalyzer (Agilent, CA, USA) was used to evaluate the quality of these cDNA libraries. Sequencing was performed as paired-end [2 × 150 base pairs (bp)] using an Illumina NovaSeq6000 (Illumina, CA, USA) by Theragen Bio sequencing service (Theragen Bio Co. Ltd., Seongnam, Gyeonggi, Korea). The ends of the reads less than Phred quality score of 20 and the adapter sequences were trimmed. In addition, the reads shorter than 50 bp were removed by using cutadapt (v2.8) (67). Aligner STAR (v.2.7.1.a) (68) was used to map filtered reads to the reference genome (mm10), and RSEM (v.1.3.1) (69) was used to estimate gene expression. FPKM (Fragments per kilobase of transcript per million mapped fragments), RPKM (Reads per kilobase of transcript per million mapped reads), and TPM (Transcripts per million) values were calculated to normalize the sequencing depth among the samples. The TCC (Tag count comparison) R package (v.1.26.0) (70) was used to normalize and compare the tag count data. To identify physiologically meaningful gene expression changes, genes with read counts lower than 10 and/or RPKM lower than 2 were filtered out. Normalization factors were calculated using the iterative DEseq2 method (71). The DEGs were identified on the basis of the P value less than 0.05 and an absolute value of fold change over 2. GO analysis was performed using Metascape (72). For GSEA, GSEA 4.1.0 software was used.