*2.9. Bioinformatic Analysis*

The reads were first trimmed for the first three nucleotides (-u 3) and adapters overlapping at least five nucleotides with the read (-O 5) using cutadapt v2.9 [27]. Reads shorter than 15 nucleotides were discarded (-m 15). Trimmed reads were mapped using bwa-mem v0.7.17-r1188 in three steps, with the minimum score output (-T) and seed length (-k) set to 15 [28]. The reads were mapped successively against a custom list of transcripts containing ribosomal RNAs (rRNAs) from rFam 14.1, mature miRNAs from miRbase 22.1, ncRNAs from ENSEMBL release 97, piRNAs from piRNA-DB v1.7.5, and cDNAs from ENSEMBL release 97 as the references. In this order, unmapped reads from each step were mapped against the next reference to assure unique attribution per RNA type. The reads were counted per transcript with Samtools v1.10 [29], and 374 miRNAs with more than 10 reads across all samples were analyzed using R package DESeq2 v1.28.1 for the differential gene expression [30]. Differences were classified as significant with a threshold of absolute value of fold change (FC) > 2 and FDR < 0.05. Principal component analysis was conducted and visualized using R package pcaExplorer v2.14.2 based on variance stabilized transformed read counts of miRNAs [31]. Heatmaps were generated using R package pheatmap v1.0.12 with default clustering parameters. KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway analyses were performed using DIANA-miRPath v3.0 with a FDR threshold < 0.05 and other default settings [32]. miRNAs function, family, disease, and regulatory proteins were analyzed using TAM 2.0 by masking cancer-related terms and keeping the other default settings [33]. Network analysis was conducted by miRTargetLink Human with strong evidence, and the resulting genes were uploaded on STRING v11.0 for the Reactome pathways [34].
