Gene expression analysis of alveolar macrophages tolerized to ozone

In this study, we examined the effect of repeated exposure to the air pollutant ozone (O3) on the transcriptome of alveolar macrophages (AMs). Using flow cytometry, we isolated purified AMs from the lungs of female C57BL6/NJ mice exposed to filtered air or 0.8 ppm O3 for four days and then collected 2 days later, with six mice per group. RNA was isolated from AMs and then subjected to RNA-seq analysis. We performed differential gene expression analysis followed by gene ontology enrichment analysis of differentially expressed genes.

Methods
Alveolar macrophages were isolated from mouse whole lung tissue using flow cytometry. Total RNA was extracted from two batches of flow-sorted AMs using QIAGEN RNAeasy kits per the manufacturer’s instructions. RNA integrity was analyzed using an Agilent Bioanalyzer. RIN values ranged from 8.9-9.7, indicating intact RNA. PolyA+ RNA libraries were prepared with the Roche Kapa mRNA stranded library preparation kit as per the manufacturer's instructions. Paired-end sequencing (50 cycles) was performed on an Illumina NovaSeq SP to a depth of >55M read pairs per sample by the UNC High-Throughput Sequencing Facility. 

Raw reads were trimmed and filtered of adapter contamination using cutadapt (Martin, 2011), and further filtered such that at least 90% of bases had a quality score of at least 20 using fastx_toolkit v0.0.14. Reads were then aligned to the reference mouse genome (mm10) (Supplementary Table 3) and GENCODE vM25 transcript annotations using STAR v2.7.7a (Dobin et al., 2013), and transcript abundance was estimated using salmon v1.5.2 (Patro et al., 2017). Differential expression between ozone-exposed vs. filtered air groups was then detected using DESeq2 v1.26.0 (Love et al., 2014) in R v3.6.0, using a design that corrected for flow cytometry batch dates. These batch effects were also removed from the VST-normalized expression values using limma v3.42.2 (Ritchie et al., 2015). All log2 fold-changes reported were shrunken using ashr 2.2-47 (Stephens, 2017). Gene ontology enrichments were then assessed using gprofiler2 v0.2.1 (Raudvere et al., 2019).

Dobin, A. et al. (2013) STAR: ultrafast universal RNA-seq aligner. Bioinformatics, 29, 15–21.

Love, M.I. et al. (2014) Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol., 15, 550.

Martin, M. (2011) Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet.journal, 17, 10.

Patro, R. et al. (2017) Salmon provides fast and bias-aware quantification of transcript expression. Nat. Methods, 14, 417–419.

Raudvere, U. et al. (2019) g:Profiler: a web server for functional enrichment analysis and conversions of gene lists (2019 update). Nucleic Acids Res., 47, W191–W198.

Ritchie, M.E. et al. (2015) Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res., 43, e47.

Stephens, M. (2017) False discovery rates: a new deal. Biostatistics, 18, 275–294.

本研究探讨了反复暴露于空气污染物臭氧（O₃）对肺泡巨噬细胞（alveolar macrophages, AMs）转录组的影响。本研究采用流式细胞术（flow cytometry），从暴露于过滤空气或0.8 ppm臭氧、连续暴露4天后再饲养2天的雌性C57BL6/NJ小鼠肺部分离纯化肺泡巨噬细胞，每组设置6只小鼠。从肺泡巨噬细胞中提取RNA并进行RNA测序（RNA-seq）分析。我们对差异表达基因进行了差异基因表达分析，随后开展了基因本体（gene ontology, GO）富集分析。

方法
采用流式细胞术从小鼠全肺组织中分离肺泡巨噬细胞。使用QIAGEN RNAeasy试剂盒，按照制造商说明书从两批流式分选的肺泡巨噬细胞中提取总RNA。使用安捷伦生物分析仪（Agilent Bioanalyzer）分析RNA完整性，所得RNA完整性数（RIN）介于8.9~9.7之间，表明RNA完整无降解。按照制造商说明书，使用罗氏（Roche）Kapa mRNA链特异性文库制备试剂盒构建PolyA+ RNA文库。由北卡罗来纳大学高通量测序中心（UNC High-Throughput Sequencing Facility）在Illumina NovaSeq SP测序平台上进行50个循环的双端测序，每个样本的测序深度均大于5500万读段对。
使用cutadapt（Martin, 2011）对原始读段进行修剪并过滤适配器污染，随后使用fastx_toolkit v0.0.14进行进一步过滤，确保至少90%的碱基质量值不低于20。使用STAR v2.7.7a（Dobin等, 2013）将读段比对至参考小鼠基因组mm10（补充表3）及GENCODE vM25转录本注释文件，随后使用salmon v1.5.2（Patro等, 2017）估算转录本丰度。使用R v3.6.0环境下的DESeq2 v1.26.0（Love等, 2014）检测臭氧暴露组与过滤空气组之间的差异表达，分析模型校正了流式细胞术的批次日期效应。使用limma v3.42.2（Ritchie等, 2015）去除方差稳定变换（VST）标准化后的表达矩阵中的批次效应。所有报道的log₂倍变化值均使用ashr 2.2-47（Stephens, 2017）进行收缩校正。随后使用gprofiler2 v0.2.1（Raudvere等, 2019）进行基因本体富集分析。

参考文献
Dobin, A. 等(2013) STAR：超快速通用RNA-seq比对工具。《Bioinformatics》，29，15–21。
Love, M.I. 等(2014) 基于DESeq2的RNA-seq数据倍变化与离散度的有偏估计。《Genome Biology》，15，550。
Martin, M. (2011) Cutadapt可从高通量测序读段中移除适配器序列。《EMBnet Journal》，17，10。
Patro, R. 等(2017) Salmon可快速且偏向性感知地估算转录本表达量。《Nature Methods》，14，417–419。
Raudvere, U. 等(2019) g:Profiler：用于功能富集分析与基因列表转换的在线服务器（2019更新版）。《Nucleic Acids Research》，47，W191–W198。
Ritchie, M.E. 等(2015) Limma助力RNA测序与芯片研究的差异表达分析。《Nucleic Acids Research》，43，e47。
Stephens, M. (2017) 错误发现率：新视角。《Biostatistics》，18，275–294。