2.7. Bioinformatics Analysis of RNA-seq Data

Raw reads were pre-processed to remove low quality regions and adapter sequences. Sequence tag preprocessing was performed according to a previously described protocol with some modification. Reads with adaptors of low quality (>50%) or a high proportion of unknown bases (>5%) were removed. Clean data were mapped to the A. cerana genome (https://www.ncbi.nlm.nih.gov/genome/?term=Apis+cerana) using TopHat software with a maximum allowance of 2 nucleotide mismatches.

The gene expression level was normalized using the fragments per kilobase of transcript per million mapped reads (FPKM) method [38]. Differential expression analysis of two groups was performed using the DESeq R package (1.18.0). The resulting p-values were adjusted using the Benjamini and Hochberg’s approach for controlling the false discovery rate. Genes with corrected p-values < 0.05 and|log2 (Fold change)| ≥ 1 were set as the threshold for significantly differential expression. Further, functional classification of the DEGs within up-regulated or down-regulated clusters was carried out using WEGO software [39], and KEGG pathway annotation by Blastall software against the KEGG (Kyoto Encyclopedia of Genes and Genomes) database (http://www.kegg.jp/). To determine the separation of expression patterns across samples, principle component analysis (PCA) and heatmap and cluster analysis (HCA) were carried out using the OmicShare tools, a free online platform for data analysis (www.omicshare.com/tools). Gene expression data for trend analysis (from 1 h to 16 h of A. cerana treatments) were normalized to log2(CK/CK) log2 (1h/CK), log2 (8h/CK), log2 (16h/CK). The trend profiles were generated by the Short Time-series Expression Miner (STEM, version 1.2.2b) software [40].