detectRUNS is a R package for the detection of runs of homozygosity (ROH/ROHom) and of heterozygosity (ROHet, a.k.a. “heterozygosity-rich regions”) in diploid genomes. ROH/ROHom were first studied in humans (e.g. McQuillan et al. 2008) and rapidly found applications not only in human genetics abut also in animal genetics (e.g. Ferencakovic et al., 2011, in Bos taurus). More recently, the idea of looking also at “runs of heterozygosity” (ROHet or, more appropriately, “heterozygosity-rich regions”) has been proposed (Wiliams et al. 2016)

Biscarini, Filippo

Cozzi, Paolo

Gaspa, Giustino

Marras, Gabriele

Online Research @ Cardiff

detectRUNS: an R package to detect runs ofhomozygosity and heterozygosity in diploidgenomesFilippo Biscarini, Paolo Cozzi, Giustino Gaspa, Gabriele MarrasIBBA-CNR, PTP, Università degli Studi di Sassari, University of Guelph  lippo.biscarini@ibba.cnr.it, gmarras@uoguelph.ca2018-02-06OverviewSample dataDetect runssliding-window-based run detectionconsecutive SNP-based run detection“Runs of heterozygosity” (a.k.a. heterozygosity-rich regions)Summary statistics on detected runsPlots: ROH-based inbreedingImporting data from external  lesReferences## Using detectRUNS 0.9.5OverviewdetectRUNS is a R package for the detection of runs of homozygosity (ROH/ROHom) and ofheterozygosity (ROHet, a.k.a. “heterozygosity-rich regions”) in diploid genomes. ROH/ROHom were rst studied in humans (e.g. McQuillan et al. 2008) and rapidly found applications not only in humangenetics abut also in animal genetics (e.g. Ferencakovic et al., 2011, in Bos taurus). More recently, theidea of looking also at “runs of heterozygosity” (ROHet or, more appropriately, “heterozygosity-richregions”) has been proposed (Wiliams et al. 2016).detectRUNS uses two methods to detect genomic runs:1. sliding-window based method:2. consecutive runs:The sliding-window based method is similar to what is implemented in the computer package Plink(Purcell et al., 2007) [see Bjelland et al., 2013 for a description]. In brief, a sliding window is used toscan the genome, and the characteristics of consecutive windows are used to determine whether aFROHSNP is or not in a run (either ROH/ROHom or ROHet). Parameters for both the sliding window and therun need to be speci ed.The “consecutive runs” method is window-free and directly scans the genome SNP by SNP. It was  rstproposed by Marras et al. (2015). Here, only parameters for the runs need to be speci ed.Besides detecting genomic runs (again, either homozygosity or heterozygosity, either sliding-windowbased or consecutive), and saving results to a data frame of individual runs, detectRUNS can:plot runs along the genome:plot runs per individualplot stacked (piled) runsplot the % of times each SNP is in a run in the population (per chromosome)Manhattan-like plot of the % of times each SNP is in a runplot mean or total run length vs number of runs per individualgenerate summary descriptive statistics on detected runscalculate inbreeding based on ROH ( ), genome-wide and chromosome-wideplot  per chromosomeThe input  les for detectRUNS are Plink ped/map  les. If one wishes to use this R package only forplots and summary statistics, output  les from Plink (.hom  les) can be easily read into detectRUNSthrough a speci c function.detectRUNS can be used with genotype data from any diploid organisms: humans, animals or plants.Sample dataTo illustrate the functionalities of detectRUNS, we use data on sheep (Ovis aries) SNP genotypes fromthe work by Kijas et al. (2016), available on-line through “Dryad” (https://goo.gl/sfAy8k). A subset withtwo breeds (“Jacobs” and “Navajo-Churro”, 100 animals) and two chromosomes (  SNPs from OAA2 and 24) was used.genotypeFilePath <- system.file(   "extdata", "Kijas2016_Sheep_subset.ped", package="detectRUNS") mapFilePath <- system.file(   "extdata", "Kijas2016_Sheep_subset.map", package="detectRUNS")Detect runsFor the detection of genomic runs, detectRUNS uses two main functions:1. slidingRUNS.run: for sliding-window-based detection2. consecutiveRUNS.run: for consecutive-SNP-based detectionInput  les are to be passed as paths to  les (e.g. /home/Documents/experiment/ le.ped/map).FROHFROH4 841sliding-window-based run detectionThe function slidingRUNS.run() accepts in input several parameters: besides the paths (or names) ofped/map  les, there are parameters related to the sliding window and parameters related to thegenomic runs.Sliding-window parameters are: windowSize, threshold (to call a SNP “in run”), minSNP (minimumnumber of homozygous/heterozygous SNP in the window), maxOppWindow (maximum number ofSNP with opposite genotype: heterozygous/homozygous) and maxMissWindow (maximum number ofmissing genotypes).Run-related parameters are: maxGap (maximum gap between consecutive SNPs, in basepairs -bps),minLenghtBps (minimum length of the run, in bps), minDensity (number of SNPs every  kilo-basepairs -kbps), maxOppRun (maximum number of opposite genotypes in the run), maxMissRun(maximum number of missing genotypes in the run).ROHet controls whether runs of homozygosity (ROH/ROHom) or of heterozygosity (heterozygosity-richregions, ROHet) will be detected. It defaults to FALSE (ROH/ROHom).slidingRuns <- slidingRUNS.run(   genotypeFile = genotypeFilePath,    mapFile = mapFilePath,    windowSize = 15,    threshold = 0.05,   minSNP = 20,    ROHet = FALSE,    maxOppWindow = 1,    maxMissWindow = 1,   maxGap = 10^6,    minLengthBps = 250000,    minDensity = 1/10^3, # SNP/kbps   maxOppRun = NULL,   maxMissRun = NULL ) consecutive SNP-based run detectionThe function consecutiveRUNS.run() has a similar structure, obviously without the sliding-windowparameters.consecutiveRuns <- consecutiveRUNS.run(   genotypeFile =genotypeFilePath,   mapFile = mapFilePath,   minSNP = 20,   ROHet = FALSE,   maxGap = 10^6,   minLengthBps = 250000,   maxOppRun = 1, x  maxMissRun = 1 )slidingRUNS.run() detected 1348 ROH; consecutiveRUNS.run() detected 1144 ROH.“Runs of heterozygosity” (a.k.a. heterozygosity-rich regions)By setting ROHet=TRUE, runs of heterozygosity (a.k.a. heterozygosity-rich genomic regions) aredetected instead. Again, the usert can choose whether to use the sliding-window or the consecutivemethod.slidingRuns_het <- slidingRUNS.run(   genotypeFile = genotypeFilePath,    mapFile = mapFilePath,    windowSize = 10,    threshold = 0.05,   minSNP = 10,    ROHet = TRUE,    maxOppWindow = 2,    maxMissWindow = 1,   maxGap = 10^6,    minLengthBps = 10000,    minDensity = 1/10^6, # SNP/kbps   maxOppRun = NULL,   maxMissRun = NULL ) consecutiveRuns_het <- consecutiveRUNS.run(   genotypeFile =genotypeFilePath,   mapFile = mapFilePath,   minSNP = 10,   ROHet = TRUE,   maxGap = 10^6,   minLengthBps = 10000,   maxOppRun = 2,   maxMissRun = 1 )slidingRUNS.run() detected 2297 ROHet; consecutiveRUNS.run() detected 1870 ROHet.Runs of homozygosity (ROH) detected using the sliding-windows method (output fromslidingRUNS.run()) will be used to illustrate summary statistics, plots and inbreeding calculations.Summary statistics on detected runsThe function summaryRuns() takes in input the dataframe with results from runs detection andcalculates a number of basic descriptive statistics on runs. Additional necessary parameters are thepaths to the Plink ped and map  les. Class  and snpInRuns  are optional arguments.summaryList <- summaryRuns(   runs = slidingRuns, mapFile = mapFilePath, genotypeFile = genotypeFilePath,    Class = 6, snpInRuns = TRUE)The returned list includes the following dataframes:summary_ROH_count_chr, summary_ROH_percentage_chr, summary_ROH_count,summary_ROH_percentage, summary_ROH_mean_chr, summary_ROH_mean_class,result_Froh_genome_wide, result_Froh_chromosome_wide, result_Froh_class, SNPinRunWe can, for instance, have a look at the number of runs per class-size (Mbps) in the two breeds: we seethat in Jacobs sheep there are 894 ROH with size up to 6 Mbps.summaryList$summary_ROH_count##       Jacobs Navajo-Churro ## 0-6      894           260 ## 6-12     101            20 ## 12-24     34            12 ## 24-48     14             7 ## >48        3             3Or, the average number of ROH per chromosome and per breed can be obtained.summaryList$summary_ROH_mean_chr##   chrom   Jacobs Navajo-Churro ## 1     2 4.069597      4.492243 ## 2    24 3.697987      2.745963The dataframe “SNPinRun” contains, for each SNP, the proportion of times it falls inside a run in anygiven population/group:head(summaryList$SNPinRun)##   SNP_NAME CHR POSITION COUNT  BREED PERCENTAGE ## 1 s47174.1   2   158066    11 Jacobs    17.1875 ## 2 s30876.1   2   167018    11 Jacobs    17.1875 ## 3 s16671.1   2   354035    12 Jacobs    18.7500 ## 4 s31991.1   2   544617    12 Jacobs    18.7500 ## 5 s33481.1   2   679324    12 Jacobs    18.7500 ## 6 s49858.1   2   901838    12 Jacobs    18.7500The summary information included in the list returned from summaryRuns() is conveniently organizedin data.frames, so that it can easily be visualized, manipulated and written out to text  les (e.g. .csv les).PlotsdetectRUNS produces a number of plots from the dataframe with runs (results from sliding-windowor consecutive scans of the genome for ROH/ROHet).The basic plot, produced by the function plot_Runs(), plots directly all runs detected in each individualagainst their position along the chromosome. Separate plots per chromosome are produced, anddi erent groups/populations are coloured di erently to visualize contrasting patterns.plot_Runs(runs = slidingRuns)Alternatively, runs can still be plotted against their position along the chromosome, but stacked on topof each other: this way, regions of the genome with an excesse of runs can easily be identi ed. In thiscase, separate plots per chromosome and per group/population are produced.plot_StackedRuns(runs = slidingRuns)Finally, the proportion of times each SNP falls inside a run in any given population/group can beplotted against their position along the chromosome, separately per group. The functionplot_SnpsInRuns() requires as arguments, besides the dataframe with detected runs, also the paths tothe original ped (for information on groups) and map (for SNP positions)  les.plot_SnpsInRuns(   runs = slidingRuns[slidingRuns$chrom==2,], genotypeFile = genotypeFilePath,    mapFile = mapFilePath)plot_SnpsInRuns(   runs = slidingRuns[slidingRuns$chrom==24,], genotypeFile = genotypeFilePath,    mapFile = mapFilePath)We can see from the plots above, that in the Jacob sheep breed a region with a “peakk”" of ROH can bespotted approximately halfway on chromosome 2 (OAR2) in the Jacob breed. This corresponds to thestrong GWAS signals found by Kijas et al. (2016) on OAR2 associated with the four-horns phenotype.To identify the position of a runs (ROH in this case) peak, e.g. from plot_SnpsInRuns(), one canconveniently use the function detectRUNS::tableRuns(): this requests as input, besides the runsdataframe, also the paths to the original ped/map  les, and the threshold above which we wantinformation on such “peaks” (e.g. only peaks where SNP are inside runs in more than 70% of theindividuals in that population/group).topRuns <- tableRuns(   runs =  slidingRuns, genotypeFile = genotypeFilePath, mapFile = mapFilePath,    threshold = 0.7)##    Group        Start_SNP          End_SNP chrom nSNP      from        to ## 1 Jacobs OAR2_130215973.1 OAR2_131721845.1     2   21 121796989 123528595The information on the proportion of times each SNP falls inside a run, can also be plotted againstSNP positions in all chromosomes together, similarly to the familiar GWAS Manhattan plots:plot_manhattanRuns(   runs = slidingRuns[slidingRuns$group=="Jacobs",],    genotypeFile = genotypeFilePath, mapFile = mapFilePath): ROH-based inbreedingFrom runs of homozygosity (ROH), individual inbreeding/consanguinity coe cients can be calculatedas:where  is the sum of the length of all ROH detected in an individual, and  is the totallength of the genome that was used.detectRUNS provide functions to calculate individual inbreeding/consanguinityFROH=FROH∑ LROHLgenome∑ LROH Lgenome##     id         group      sum Froh_genome ## 1   H1 Navajo-Churro 16593339  0.05705571 ## 2  H10 Navajo-Churro 45444521  0.15625966 ## 3  H11 Navajo-Churro 18389476  0.06323168 ## 4 H114        Jacobs 71163071  0.24469215 ## 5 H115        Jacobs 64303071  0.22110424 ## 6 H116        Jacobs 62350601  0.21439073The parameter “genome_wide” (which defaults to TRUE) can be used to obtaininbreeding/consanguinity estimates on a per-chromosome basis (by setting “genome_wide=FALSE”)Inbreeding levels can be plotted by group, for example:Importing data from external  lesResults on runs (typically ROH) from external software can be imported into detectRUNS to produceplots, tables and summary statistics. Current options include:ROH from Plink (.hom  le from --homozyg )ROH from BCFtools (output from the roh  option)runs dataframes from detctRUNS written out to  lesThrough the parameter program  the user can select from which source the output  le is coming. Asan illustration, we read in results from detectRUNS saved out to a .csv  le:savedRunFile <- system.file(   "extdata", "Kijas2016_Sheep_subset.sliding.csv", package="detectRUNS") runs <- readExternalRuns(inputFile = savedRunFile, program = "detectRUNS") head(runs)##    group   id chrom nSNP     from       to lengthBps ## 1 Jacobs H114     2   76  1606489  5215635   3609146 ## 2 Jacobs H114     2   44  5798140  8518359   2720219 ## 3 Jacobs H114     2  148 21979508 30849683   8870175 ## 4 Jacobs H114     2   45 35582012 38744795   3162783 ## 5 Jacobs H114     2   15 39198369 40023789    825420 ## 6 Jacobs H114     2   16 52327990 53043790    715800ReferencesBjelland, D. W., K. A. Weigel, N. Vukasinovic, and J. D. Nkrumah. “Evaluation of inbreedingdepression in Holstein cattle using whole-genome SNP markers and alternative measures ofgenomic inbreeding.” journal of dairy Science 96, no. 7 (2013): 4697-4706.Ferencakovic, Maja, Edin Hamzic, Birgit Gredler, Ino Curik, and Johann Sölkner. “Runs ofhomozygosity reveal genome-wide autozygosity in the Austrian Fleckvieh cattle.” AgriculturaeConspectus Scienti cus (ACS) 76, no. 4 (2011): 325-329.Kijas, James W., Tracy Had eld, Marina Naval Sanchez, and Noelle Cockett. “Genome‐wideassociation reveals the locus responsible for four‐horned ruminant.” Animal genetics 47, no. 2(2016): 258-262.Marras, Gabriele, Giustino Gaspa, Silvia Sorbolini, Corrado Dimauro, Paolo Ajmone‐Marsan,Alessio Valentini, John L. Williams, and Nicolò PP Macciotta. “Analysis of runs of homozygosity andtheir relationship with inbreeding in  ve cattle breeds farmed in Italy.” Animal genetics 46, no. 2(2015): 110-121.McQuillan, Ruth, Anne-Louise Leutenegger, Rehab Abdel-Rahman, Christopher S. Franklin,Marijana Pericic, Lovorka Barac-Lauc, Nina Smolej-Narancic et al. “Runs of homozygosity inEuropean populations.” The American Journal of Human Genetics 83, no. 3 (2008): 359-372.Purcell, Shaun, Benjamin Neale, Kathe Todd-Brown, Lori Thomas, Manuel AR Ferreira, DavidBender, Julian Maller et al. “PLINK: a tool set for whole-genome association and population-basedlinkage analyses.” The American Journal of Human Genetics 81, no. 3 (2007): 559-575.Williams, John L., Stephen JG Hall, Marcello Del Corvo, K. T. Ballingall, L. I. C. I. A. Colli, P. A. O. L. O.Ajmone Marsan, and F. Biscarini. “Inbreeding and purging at the genomic Level: the Chillinghamcattle reveal extensive, non‐random SNP heterozygosity.” Animal genetics 47, no. 1 (2016): 19-27.

detectRUNS: Detect runs of homozygosity and runs of heterozygosity in diploid genomes

http://orca.cf.ac.uk/108906/1/detectRUNS_%20an%20R%20package%20to%20detect%20runs%20of%20homozygosity%20and%20heterozygosity%20in%20diploid%20genomes.pdf

detectRUNS: Detect runs of homozygosity and runs of heterozygosity in diploid genomes

Abstract

Similar works

Full text

Available Versions

Online Research @ Cardiff