A new rhesus macaque assembly and annotation for next-generation sequencing analyses

Genomic DNA sequencing

We obtained genomic DNA from the reference rhesus macaque (animal 17573) [1] and performed whole genome Illumina sequencing on a GAIIx instrument, yielding 107 billion bases of sequence data. We deposited these sequences in the Sequence Read Archive (SRA) under accessions [GenBank:SRX112027, GenBank:SRX113068, GenBank:SRX112904]. In addition, we used a human exome capture kit (Illumina TruSeq Exome Enrichment) to enrich exonic sequence from the reference rhesus macaque genomic DNA. Illumina HiSeq2000 sequencing of exonic fragments from this animal generated a total of 17.7 billion bases of data. We deposited these sequences in the SRA under accession [GenBank:SRX115899].

Contig and scaffold assembly

We assembled the combined set of Sanger (approximately 6× coverage), Illumina whole genome shotgun (approximately 35× coverage) and exome reads using MaSuRCA (then MSR-CA) assembler version 1.8.3 [9]. We pre-screened and pre-trimmed the Sanger data with the standard set of vector and contaminant sequences used by the GenBank submission validation pipeline. The MaSuRCA assembler is based on the concept of super-read reduction whereby the high-coverage Illumina data is transformed into 3-4× coverage by much longer super-reads. This transformation is done by uniquely extending the Illumina reads using k-mers and then combining the reads that extend to the same sequence. We transformed the exome sequence data from the reference animal into a separate set of exome super-reads. We then used these exome super-reads along with Sanger and whole genome shotgun Illumina data in the assembly. The exome super-reads were marked as non-random and therefore were excluded from the contig coverage evaluation step that is designed to distinguish between unique and repeat contigs.

Chromosome assembly steps

A flowchart (Figure 1) illustrates the overall process of assembly and annotation.

We used BLAST + (version 2.2.25) for all BLASTn [10] alignments. We used default parameters for BLASTn alignments with the following exceptions:

We have deposited our new rhesus macaque assembly (MacaM_Assembly_v7) in NCBI’s BioProjects database under accession [GenBank:PRJNA214746].

Chromosome assembly validation

We used genomic DNA from the reference rhesus macaque (animal 17573) to create a 400 bp library according to manufacturer’s instructions (Ion Torrent, Personal Genome Machine). We sequenced this library on an Ion 318 chip and deposited 1.5 billion bases of sequence in the SRA under accession [GenBank:SRR1216390]. To independently assess genome assembly, we aligned these Ion Torrent reads (which were not included in our assembly) against rheMac2, CR_1.0 and MacaM assemblies with TMAP 4.0 [18].

RNA sequencing and transcript assembly

We extracted RNA from 11 samples using standard methods and performed sequencing with an Illumina Genome Analyzer IIx. We sequenced RNA from the cerebral cortex from 6 different animals at 76 bp, (single-end reads) and deposited the sequences in the SRA under accessions [GenBank:SRX099247, SRX101205, SRX101272, SRX101273, SRX101274 and SRX101275]. We sequenced RNA from the caudate nucleus of one animal at 76 bp (paired-end reads) and deposited the sequences at SRA under accession [GenBank:SRX103458]. We sequenced RNA from the caudate nucleus, cerebral cortex, thymus, and testis from a single rhesus macaque, 001T-NHP, at 76 bp for the caudate nucleus and at 100 bp for the other tissues (paired-end reads for all samples) and deposited the sequences in the SRA under the accessions [GenBank:SRX101672, SRX092157, SRX092159 and SRX092158], respectively.

To filter out genomic contamination, we aligned reads against human RefSeq mRNA transcripts using BLASTn [10]. For 76 bp reads, we filtered out sequences if they had an alignment length of 70 bp or less with human transcripts. For 100 bp reads, we filtered out sequences if they had an alignment length of 90 bp or less with human transcripts. For paired-end reads, we also removed a read if its mate was removed. We assembled filtered reads for each sample using Velvet-Oases [19, 20]. We used K-mer values of 29 for the six samples with single reads and 31 for the remaining samples with paired-end reads. We set the coverage cutoff and expected coverage to ‘auto’. We set the minimum contig length to 200 bp. We used default parameters for Oases.

We obtained 369,197 de novo transcripts using the Velvet/Oases pipeline. We deposited transcripts in the NCBI’s Transcriptome Shotgun Assembly database under accessions [GenBank:JU319578 - JU351361; GenBank:JU470459 - JU497303; GenBank:JV043150 - JV077152; GenBank:JV451651 - JV728215] (Additional file 2). We also used reference-guided transcriptome assembly to identify rhesus transcripts. We performed spliced alignment of the rhesus RNA-seq reads to the MaSuRCA assembly using TopHat2 [21] (version 2.0.8b, default parameters). We provided the resulting BAM file of the read alignments as input to Cufflinks2 [22] (version 2.02, default parameters) for reference guided transcriptome assembly.

To identify rhesus orthologs of human genes, we used the BLASTx [23] program to align conceptual translations from the assembled rhesus transcripts against human reference proteins. We used the top hit in the annotation. We did not use a single set of cutoff values to identify orthologs. Instead, alignment lengths and percent similarity were manually inspected (see Annotation, procedure 1 for rationale). We identified a total of 11,712 full-length proteins representing 9,524 distinct genes with full length coding regions using the de novo and reference-guided methods described above.

Annotation

We produced a GTF file (MacaM_Annotation_v7.6.8, Additional file 3) that serves as our annotation of the MacaM assembly. We used the following procedures to generate this file:

We were able to identify complete protein models for 16,052 rhesus genes (Additional file 3) from a human gene target list of 19,063 named protein-coding genes (Additional file 1).

Protein comparison

We downloaded the most recent human assembly (GRCh38) and GFF3 annotation from NCBI on February 6, 2014. To obtain a list of genes that contained only a single isoform, we filtered the GFF3 file so that only Gene IDs linked to a single RefSeq mRNA accession were retained. This procedure resulted in a list of 11,148 Gene IDs. We downloaded the most recent rhesus assembly, rheMac2, and the GFF3 annotation from NCBI on February 6, 2014 [30]. We then created a list of genes that was common to GRCh38, rheMac2, MacaM by determining if the gene name from the GRCh38 (human) single isoform list was also present in the NCBI annotation of rheMac2 and our new annotation of the MacaM assembly. We used the Cufflinks2 [22] tool, gffread, to obtain protein sequences for each of these genes in the two rhesus genomes. We then aligned the rheMac2/NCBI and MacaM proteins against their human protein orthologs using the EMBOSS [31] global alignment tool, needle v.6.3.1. Once the alignments were done, we used custom scripts to compile the results into a summary table (Additional file 4). We used this table to calculate mean values for Identity, Similarity and Gaps.

We previously identified a gene that was misassembled in rheMac2 – the Src homology 2 domain containing E (SHE) gene [3]. To compare models of the SHE gene from different assemblies and annotations, we attempted to find the proteins with this designation from several rhesus macaque annotations. Using BLASTp [32] with the human SHE protein as a query, we found a protein annotated by NCBI in the rheMac2 assembly as gene name “LOC716722”, protein definition “PREDICTED: SH2 domain-containing adapter protein E-like [Macaca mulatta]” under the accession XP_002801853.1. Ensembl annotated the SHE gene in rheMac2 with gene identification ENSMMUG00000022980. We obtained the protein sequence associated with this gene model (ENSMMUT00000032345) for comparison with other rhesus macaque annotations. We accessed the RhesusBase database (http://www.rhesusbase.org/) on June 17, 2014 and searched under Gene Symbol for “SHE” and Gene Full Name for “Src homology 2 domain containing E”. In both cases the message returned was “No gene found!”. In an attempt to see if the SHE gene was annotated in the CR_1.0 genome, we downloaded the CR.pep.fa file, containing proteins derived from this genome from the GIGA database site (http://gigadb.org/dataset/100002) containing sequences and annotations related to CR_1.0 on June 17, 2014. In this file, we identified a protein with the ENSEMBL identifier ENSMMUP00000030263, with a partial match to the rhesus macaque SHE protein we identified in MacaM. We searched the ENBSEMBL website for “ENSMMUP00000030263” and received a message directing us to the rhesus macaque SHE gene. This protein was also submitted to NCBI under accession EHH15290.1 where it is described as “hypothetical protein EGK_01357, partial [Macaca mulatta]”. The various putative rhesus macaque SHE proteins were aligned against the human SHE protein (NP_001010846.1).

RNA expression analysis

We aligned raw reads from three samples: testes, thymus and caudate nucleus of the brain to both the rheMac2 and the MacaM assemblies using TopHat [21] (version 2.0.8) with default parameters. We analyzed those alignment files using Cufflinks2 [22], specifically Cuffdiff, to generate normalized expression (FPKM) values. The samples from which these mRNA sequences are derived are described above in the section on RNA sequencing and transcript assembly (accessions GenBank:SRX092158, GenBank:SRX092159 and GenBank:SRX101672).

We also obtained 60 peripheral blood mononuclear cell (PBMC) samples from rhesus macaques in a social hierarchy experiment performed at the Yerkes National Primate Research Center. We individually subjected 10 macaques that were dominant in the social hierarchy and 10 subordinate macaques to a human intruder as a stressor. Whole blood was then collected at several time points using a BD CPT vacutainer which allow for the collection of PBMCs. RNA was purified from PBMCs using the RNeasy kit (QIAGEN, Valencia CA). We prepared libraries using standard TruSeq chemistry (Illumina Inc., San Diego CA) and sequenced them on an Illumina Hi-Seq 1000 as 2 × 100 base paired-end reads at the Yerkes NHP Genomics Core Laboratory (http://www.yerkes.emory.edu/nhp_genomics_core/). Sequences were deposited at NCBI under accessions [GenBank:SAMN02743270 - SAMN02743329]. We mapped reads with STAR [33] (version 2.3.0e) to both the rheMac2 and MacaM genomes, using the reference annotations as splice junction references. rheMac2 annotations were obtained from UCSC. We discarded un-annotated non-canonical splice junctions, non-unique mappings and discordant paired-end mappings. We performed transcript assembly, abundance estimates and differential expression analysis with Cufflinks2 (version 2.1.1) and Cuffdiff 2 [22]. We determined differentially expressed transcripts for pair-wise experimental group comparisons with an FDR-corrected p-value (q-value) <0.05. We compared differences in unique read counts and mapping percentages between rheMac2 and MacaM assemblies using a paired T-test.

Statement of ethical approval

Materials used in these studies were from animal work performed under Institutional Animal Care and Use Committee approval from the University of Nebraska Medical Center, Oregon Health and Sciences University and Yerkes National Primate Research Center. Animal welfare was maintained by following NIH (Public Health Service, Office of Laboratory Animal Welfare) and USDA guidelines by trained veterinary staff and researchers under Association for Assessment and Accreditation of Laboratory Animal Care certification, insuring standards for housing, health care, nutrition, environmental enrichment and psychological well-being. These met or exceeded those set forth in the Guide for the Care and Use of Laboratory Animals from the National Research Council of the US National Academy of Sciences.

Assembly of MacaM

We combined Sanger sequences from the reference rhesus macaque with newly generated Illumina whole genome and exome sequences from the same animal and created new de novo contig and scaffold assemblies using the MaSuRCA assembler [9]. Although there are many procedures for generating scaffolds, they all produce misassemblies with greater or lesser frequency [34–37]. To identify and correct misassembled scaffolds and assign scaffolds to chromosomes, we used independent mapping data and preliminary scaffold annotation. To properly assign genomic segments, we introduced 395 breaks in the assembled scaffolds. There were a total of 2312 scaffolds assigned to chromosomes.

To guide the splitting of misassembled scaffolds and placement of scaffolds on chromosomes, we used BLASTn [10] to align human exons with rhesus scaffolds, identified rhesus radiation hybrid markers [11, 12] within scaffolds, and used rhesus FISH maps to identify areas of synteny as well as species differences in chromosome structure [5, 13, 14].

We term our new assembly: MacaM.

Annotation of MacaM

Future updates to the assembly and annotation will be made available here: http://www.unmc.edu/rhesusgenechip/index.htm#NewRhesusGenome.

Comparison of MacaM with rheMac2 and CR_1

To independently assess the completeness and accuracy of the three rhesus macaque assemblies, we aligned Ion Torrent reads from the reference rhesus macaque against each assembly. These reads had not been used in any of the assemblies, including MacaM. We were able to align 93, 94 and 98% of these reads to the rheMac2, CR_1.0 and MacaM assemblies, respectively. This suggests that the MacaM assembly is the most complete and accurate of the three available rhesus macaque assemblies.

We were able to annotate genes in MacaM whose exons had been split by misassemblies in rheMac2 [3] but which were contiguous in the MacaM assembly (Figure 2). We used this case (the SHE gene) to assess the annotations for rheMac2 and CR_1.0 (Figure 3). A search for SHE at the RhesusBase database indicated that no gene was found. We were able to identify proteins that appeared to correspond to the SHE protein in MacaM, rheMac2 (both NCBI and Ensembl annotations) and CR_1.0. We aligned these sequences against the human SHE protein sequence. The MacaM sequence was highly similar to the human SHE protein sequence. We found that a portion of the protein sequence encoded by exon 1 was missing in CR_1.0. The protein sequence encoded by exon 3 was missing from the NCBI annotation of rheMac2. For the Ensembl annotation of rheMac2 and CR_1.0, spurious sequence, presumably derived from intronic sequence, was substituted for the correct protein sequence encoded by exon 3. A correct protein sequence encoded by exon 3 was derived from MacaM because the scaffold containing exon 3 was correctly assembled which was not the case for rheMac2 [3] or, apparently, CR_1.0 (Figure 2). This demonstrates an important principle: better genome assemblies permit better annotations.

mRNA expression studies with MacaM

To compare the performance of the rheMac2 and MacaM genomes on a ‘real-world’ mRNA-seq dataset, we examined the effects of an acute stressor on a background of chronic stress in social-housed female rhesus macaques. Social subordination in macaques is a natural stressor that produces distinct stress-related phenotypes and chronically stressed subordinate subjects [40]. We drew blood from both dominant and subordinate animals; following this, we exposed both to an acute stressor (human intruder paradigm [41]) and then drew blood at different time points after exposure. mRNA-seq read-mapping was significantly higher with MacaM as compared to rheMac2 (mean 3.1 × 10⁷ vs. 2.6 × 10⁷, p <0.0001) (Figure 4A). Similarly, mRNA-seq read-mapping percentages were significantly higher using the MacaM assembly than the rheMac2 assembly (mean 85.2% vs. 70.0%, p <0.0001 (Figure 4B). When we compared mRNA expression levels with one set of animals at two time points (Figure 4C), we detected many more Differentially Expressed Genes (DEGs) with the MacM genome than with rheMac2. We observed this same pattern in 8 additional comparisons (Figure 5).

Reviewer #1: Professor Lutz Walter

Zimin and co-workers re-sequenced the Indian rhesus macaque individual that was used for the first genome assembly (rheMac2). They combined their Illumina sequences (whole genome as well as exome sequencing) with the original Sanger sequences (of rheMac2) to produce a new rhesus macaque genome assembly which they called MacaM. These efforts resulted in a weighted average contig size being twice as large as the original rheMac2 assembly and five times larger than the CR_1.0 (Chinese rhesus macaque) assembly and in correction of several misassemblies. Inclusion of Illumina RNA-seq data from different tissues (cerebral cortex, caudate nucleus, thymus, testis) improved gene annotations. Further, the authors used RNA-seq data from peripheral blood mononuclear cells (PBMCs) to demonstrate the superior accuracy and capability of the MacaM assembly to align reads from RNA-seq experiments.

The previous rhesus macaque genome assemblies and gene annotations urgently needed an update, which Robert Norgren’s group provided at the right time. Their data will help researchers to use the rhesus macaque genome much more efficiently.

Authors’ response: We thank Professor Walter for his kind evaluation.

Authors’ response: For the mRNA-seq expression comparison, we initially did not provide data on the cerebral cortex because we provided data from another region of the brain, the caudate nucleus, and wanted to emphasize results from widely varying tissues. However, we can understand the desire to include all four samples in this table. We have therefore added values for the cerebral cortex to Table 6.

Further, RNA-seq from 26 animals should contain a lot of SNP data. Is there also an increase in the ability to align SNPs on the MacaM assembly compared to rheMac2?

Authors’ response: Yes, this is on our list of things to do. We expect to continue to update and improve both the MacaM assembly and annotation based on user input. Future updates to the assembly and annotation will be posted here: http://www.unmc.edu/rhesusgenechip/index.htm#NewRhesusGenome. We have also added this link to the text.4) Statistical analyses of the data shown in Figures 4C and 5 are missing.

Authors’ response: We appreciate Professor Walter’s comment on this issue.

Reviewer #2: Dr. Soojin Yi

In this manuscript the authors provide a new sequencing and assembly of a rhesus macaque using the animal used to generate the widely used RheMac1 assembly. Rhesus macaque is a promising nonhuman primate model species and the utility of an improved rhesus macaque genome assembly for biomedical and comparative studies is tremendous. This paper contains abundant merits and I look forward to using the new assembly for research and education in the near future. The obvious merits of this work include 1) combining Sanger and Illumina reads for assembly 2) validation of assembly with Ion Torrent reads 3) direct comparison of the efficiencies of different assemblies for differential gene expression studies. The RNA sequencing and transcript assembly strategy as well as annotation strategy are excellent.

Authors’ response: We thank Dr. Yi for her generous evaluation of our work.

I have a few technical comments most of which are for clarifications or details that will be helpful for future users.

p.2 and throughout the manuscript: the authors refer to the rheMac3 assembly as CR_1.0. It may help avoiding confusion if the authors used a more widely used term CR_1.0/rheMac3 (as in UCSC browser).

Authors’ response: We agree with Dr. Yi that we should note that CR_1.0 is termed rheMac3 in the UCSC genome browser and have now done so in the text. However, this is a bit of a controversial issue as only the original submitter of an assembly is technically allowed to issue updates. Baylor may wish to update their assembly, rheMac2, and would normally refer to their update as rheMac3. Because BGI has termed their assembly rheMac3, there may be some confusion if/when Baylor attempts to name their next update. We named our assembly MacaM to avoid contributing to this confusion.

p.5/Table 1: I would appreciate more details on the logic underlying the proposed new chromosome nomenclature.

Authors’ response: Chimpanzees, gorillas, and orangutans have the same general chromosomal structure as humans, with one notable exception. The human chromosome 2 appears to be the result of a fusion event that occurred during hominid evolution. Thus, the great apes and rhesus macaques both have two chromosomes that roughly correspond to the short and long arms of human chromosome 2. In the great apes, these two chromosomes are referred to as 2a and 2b. We adopted the same nomenclature for rhesus macaques to make comparisons between different primates easier. We have added additional text to our manuscript to make this point more clearly.

Adding chromosome nomenclature in other ape species may help the readers to see the merit of the proposed nomenclature system.

Authors’ response: We have added chimpanzee, gorilla and orangutan to Table 1 as requested.

p.6 identifying orthologs of human genes: as described currently it is not clear whether the reciprocal-best-blast hit strategy is used; the reciprocal step using human reference proteins against the rhesus macaque transcripts is not specified.

Authors’ response: We did not use a reciprocal-best-blast hit strategy to identify orthologs of human genes.

Also the cutoff values need to be given (e.g. e-value coverage similarity criteria).

Authors’ response: A single set of cutoff values was not used because different types of genes evolve at different rates. Our rationale for this is similar to that used to determine whether sim4cc and GMAP gene models were correct. We have added text to the manuscript to make this clearer.

p.7 p.9 discussions on protein sequence annotation and exon misassembly: The authors use the example of the Src homology 2 domain containing E (SHE) gene to demonstrate the utility of MacaM assembly and how different annotations differ in their qualities.

However I failed to find a mention of this gene until this specific section. It will be nice to see why this specific example has been selected.

Authors’ response: We agree that our discussion of the misassembly of the SHE gene begins a bit abruptly. We selected this gene because we had previously identified it as misassembled in another paper [3]. We have added text to the Protein Comparison section of the Methods to make this clear.

Did they perform a global comparison between MacaM versus rheMac2? If they did it will be nice to see a table of such results (beyond the Table 5 which simply indicates identity/similarity/gaps) something like the mismatch of exon numbers and so on would be useful) as well as the list of curated/improved gene sets. Even though this is an optional comment I feel that detailed discussion on such results would also be a nice addition to support the authors’ conclusion that the MacaM is a better assembly and thus offers a better annotation.

Authors’ response: No, we did not do a global comparison of gene assembly between MacaM and rheMac2. We discovered a number of misassembled genes in rheMac2, described in [3], in the course of attempting to annotate that assembly. I cannot say that we have found them all. Automation only takes one so far in annotation. To confirm that a given gene model is correct in one assembly and wrong in another takes a tremendous amount of manual effort, unfortunately.

p.9 please provide the full name for the gene MAMU-F. People who are not directly working with MHC class genes are not likely to understand.

Authors’ response: This was an oversight on our part. We thank Dr. Yi for pointing it out. We added the full gene name for MAMU-F to the text.

p.9 Comparison of MacaM with rheMac2 and CR_1: I understand what the authors are saying but I would remove or edit the three sentences following (Table 4) as this discussion is not directly validated with the identification of misassemblies in rheMac2.

Authors’ response: We have removed these three sentences, as suggested.

p.10 mRNA expression studies with MacaM: the authors conclude that MacaM assembly provides biologically relevant information while rheMac2 doesn’t. I agree that the results described supports “MacaM assembly provides biologically relevant information” but I am not seeing the support for “while rheMac2 doesn’t”. Perhaps they authors can elaborate.

At any rate these new data sound very interesting and perhaps the authors are preparing to publish them independently. If the authors decide not to include this data in the current manuscript so that they can publish in another paper I feel it is justified.

Authors’ response: We acknowledge that our approach may miss genes present in rhesus macaques but not in humans. Gene prediction programs are improving. So, it may be possible that additional genes will be identified using purely computational methods.

Reviewer #3: Dr. Kateryna Makova

This reviewer provided no comments for publication.