6598 Nucleic Acids Research, 2011, Vol. 39, No. 15 reciprocal liftOver of the genomic coordinates of the corresponding human annotation. Speciﬁcally, to obtain rhesus macaque protein-coding gene annotation, we downloaded human protein-coding gene annotation from Ensembl (version 57) (35). From here, we obtained the genomic coordinates of exons (50-UTR, coding region, 30-UTR), introns and gene boundaries in the rhesus macaque genome using reciprocal liftOver with default parameters and requiring the length of the matched sequence to be >70 and <130% of the query sequence length (34). Take human region A as one example. For human region A, using default parameters and requiring the length of the ortholog region to be between 70% and 130% of the length of region A, we ﬁrst used liftOver to ﬁnd region A’s ortholog region (region B) in the macaque genome. We then used the same settings to ﬁnd region B’s ortholog region (region C), in the human genome. Finally, we require regions A and C to be located within 100 nt from one another in the human genome. To annotate rhesus macaque miRNA, we downloaded human miRNA annotation from miRBase (Version 12) (36–38) and used reciprocal liftOver, with the same settings used for protein-coding gene annotation, to obtain the corres- ponding coordinates in the rhesus macaque genome. The miRNA expression quantiﬁcation method is described elsewhere (24). The miRNA expression level from ﬁve macaque tissues of this study measured by Illumina sequencing platform are listed in Supplementary Table S14. For other ncRNA (rRNA, tRNA, snoRNA, snRNA, scRNA and misc-RNA), we created a combined human ncRNA annotation by taking the union of ncRNA annotations downloaded from UCSC (34), Refseq (39) and Ensembl (35). We then used reciprocal liftOver, with the same settings used for protein-coding gene annotation, to obtain the corresponding coordinates in the rhesus macaque genome. Mouse and ﬂy genome annotation The mouse genome (NCBI37/mm9) and the ﬂy genome (BDGP R5/dm3) were downloaded from UCSC (34). Mouse and ﬂy transposon/repeat annotations were downloaded from UCSC mouse and ﬂy RepeatMasker track (34). For mouse annotation, we classiﬁed LINE, SINE, LTR and DNA as ‘transposons’ and the rest of the repetitive sequence types, such as simple repeats, as ‘other repeats’. Mouse protein-coding gene annotation was downloaded from Ensembl (version 57) (35). Fly protein-coding gene annotation was downloaded from UCSC FlyBase Genes track (34). Mouse and ﬂy miRNA annotations were downloaded from miRBase (Version 14) (36–38). For other mouse ncRNA (rRNA, tRNA, snoRNA, snRNA and misc-RNA), we created a combined mouse ncRNA annotation by taking the union of ncRNA annotations downloaded from Refseq (39) and Ensembl (35). For other ﬂy ncRNA (rRNA, tRNA, snoRNA, snRNA, ncRNA), we created a combined ﬂy ncRNA annotation by taking the union of ncRNA annotations downloaded from UCSC ﬂyBaseNoncoding Track (39) and Ensembl (35). Small RNA sequence analysis For each sample, we ﬁrst removed all sequences mapped to annotated ncRNA regions deﬁned using the methods described above. The length and the 50-nt distributions of all mapped reads and sequence tags (sequence reads after collapsing identical sequences) are shown in Figures 1A, B, 2A, B, 3A and B, the ﬁrst two columns of Supplementary Figure S3. The length and the 50- nt distributions of sequence reads and tags, excluding those mapped to known ncRNA, are shown in Figures 1C, D, 2C, D, 3C, D, the third and the fourth columns of Supplementary Figure S3. The remaining small RNA sequences with length distributions of between 24 and 32 nt (pilRNA) were used in further analyses. Speciﬁcally, we tested these small RNA with respect to: length distribution, 50-nt bias, genomic context, clustering along the genome and presence of sequence features compatible with the ping-pong biogenesis model. Sequence length distributions were compared with piRNA sequence length distributions found in three mouse piRNA data sets: adult mouse Mili-IP (16), adult mouse Miwi-IP (16) and 16.5dpc mouse Miwi-2 IP (7) (Figure 3 and Supplementary Figure S5). The 50-nt bias was measured as a percentage of sequences with a certain nucleotide at the 50-end among all pilRNA. The signiﬁcance of a 50-nt bias was tested by comparing it to the genome average, using a binomial test. For pilRNA clusters, the 50-nt bias was estimated, following the criteria described below, as a percentage of sequences with a certain nucleotide at the 50-end among all pilRNA that clustered along the genome. As these criteria include minimal percentage of small RNA with 50 uridine, the nucleotide bias within clusters is not an independent indicator of piRNA-like properties. To investigate whether pilRNA cluster along the genome, we used a 3 kb long sliding window with 0.5 kb steps to identify regions containing at least 30 uniquely mapped small RNA sequences and showing at least 60% 50 uridine bias. For the mouse small RNA data sets [GSE21630 (28)] and for all but 4 of the OSS samples for the ﬂy small RNA data sets, we required 10 unique mapped sequences instead of 30, due to relatively low sequencing coverage. Identiﬁed windows located within 1 kb distance from each other were further merged. To identify numbers of windows that expected to be found by chance, we randomized small RNA locations along the genome and repeated the analysis with the same requirements 1000 times. According to the suggested piRNA biogenesis model, the ping-pong model, there are two groups of piRNA populations: primary piRNA and secondary piRNA (1,7,8). Primary piRNAs have 50 bias for uridine and do not have nucleotide bias at position 10. Secondary piRNAs, generated by primary piRNA guided cleavage, have a bias for adenine at position 10 and do not have 50 bias. Based on this knowledge, we tested the ping-pong biogenesis model by checking the complementarity, at dif- ferent offsets, between putative primary piRNA and putative secondary piRNA. Here, small RNA sequences (24–32 nt), excluding ncRNA, were divided into two