DOI: 10.6084/m9.figshare.9900401
title: Genome-wide genic copy numbers for nine Macaca species based on 1kb windows

description: This is the genome-wide copy numbers (CNs) for every 1kb window in the genomes of nine Macaca species, including the Chinese rhesus (<i>M. mulatta lasiota</i>, CR), cynomolgus (<i>M. fascicularis</i>, CE), Tibetan (<i>M. thibetana</i>, TM), stump-tailed (<i>M. arctoides</i>, SM), southern pig-tailed (<i>M. nemestrina</i>, PM), Japanese (<i>M. fuscata</i>,<i> </i>JM), Taiwanese (<i>M. cyclopis</i>, TwM), Barbary (<i>M. sylvanus</i>, BM), and lion-tailed (<i>M. Silenus</i>, LM) macaques.<br>Methods: We employed FastQC (v0.11.8) (http://www.bioinformatics.babraham.ac.uk/projects/fastqc) to do quality control checks on raw data of resequencing genomes, then used Trimmomatic (v0.36) (Bolger et al. 2014) to filter and trim the reads. The cleaned reads were aligned to the Mmul_8 reference genome (Zimin et al. 2014) using BWA mem (Li and Durbin, 2009).Two steps were implemented to estimate CNs with the fastCN pipeline, which uses read depth information. The program fastCN was designed to efficiently estimate genome copy number from short read data (https://github.com/KiddLab/fastCN) (Pendleton et al. 2018). This method is built upon the mrsFAST aligner (Hach et al. 2014), and divides reads into 36-bp subreads and determines all possible matching locations on the reference genome with fewer than two substitutions, then reports per-bp read depth in an efficient compressed binary format. First, we performed the GC correction step using custom-defined control regions. The aim of GC correction is to remove the GC bias introduced by PCR during library preparation and sequencing. However, due to the lack of suitable control regions across macaques, we created a two-step process to retrieve copy number invariable (or control) regions in the diploid genome. The second step included the estimation of genome-wide CNs based on RDs using aligned BAM files based on windows containing one-kilobase (kb) unmasked, non-gap positions. RD values were converted to CNs using a correction factor (CF) calculated from average RD of the control regions. The calculation function is as follows:CF = RDctl / 2CN = RD / CFwhere CF stands for the correction factor, RD represents the read depth of specific genomic window, and RDctl is the mean read depth of the control region. Unplaced contigs were merged as ‘chrUn’ in data processing to decrease the CPU time. Bolger AM, Lohse M, Usadel B (2014) Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30: 2114-2120.Hach F, Sarrafi I, Hormozdiari F, Alkan C, Eichler EE, et al. (2014) mrsFAST-Ultra: a compact, SNP-aware mapper for high performance sequencing applications. Nucleic acids research 42: W494-W500Li H, Durbin R (2009) Fast and accurate short read alignment with Burrows–Wheeler transform. bioinformatics 25: 1754-1760.Pendleton AL, Shen F, Taravella AM, Emery S, Veeramah KR, et al. (2018) Comparison of village dog and wolf genomes highlights the role of the neural crest in dog domestication. BMC biology 16: 64Zimin AV, Cornish AS, Maudhoo MD, Gibbs RM, Zhang X, et al. (2014) A new rhesus macaque assembly and annotation for next-generation sequencing analyses. Biology direct 9: 1-15.