--- title: "AncestryMapper Overview" author: "Eoghan O'Halloran, Tiago Magalhães, Darren J. Fitzpatrick" date: "`r Sys.Date()`" output: rmarkdown::html_vignette: toc: yes vignette: > %\VignetteIndexEntry{AncestryMapper Overview} %\VignetteEngine{knitr::rmarkdown} \usepackage[utf8]{inputenc} --- Introduction ============ AncestryMapper is an R package that implements the methods described in "HGDP and HapMap Analysis by Ancestry Mapper Reveals Local and Global Population Relationships" Magalhães TR, Casey JP, Conroy J, Regan R, Fitzpatrick DJ, et al. PLoS ONE 7(11): e49438. (2012) Ancestry Mapper assigns genetic ancestry to an individual and allows the study of relationships between local and global populations. The method gives each individual an Ancestry Mapper Id (AMid), a genetic identifier comprising genetic coordinates that correspond to its relationship to various reference populations. The AMid metrics have intrinsic biological meaning and provide a tool to measure genetic similarity between world populations. Package Functions ================= The package consists of four functions: - **calculateAMidsArith**: calculates and assigns Ancestry Mapper Ids (AMids) to each individual using a new, more precise, arithmetic method - **calculateAMids**: calculates and assigns Ancestry Mapper Ids (AMids) identical to the older versions of AncestryMapper - **plotAMids**: produces a heatmap representation of AMids - **createMedoid**: constructs a reference from sample inputs - **refAdd**: adds user-supplied references or data to a reference file, containing information for order, labelling and color **calculateAMidsArith** -------------- For each individual, **calculateAMidsArith** computes the genetic distances amongst that individual and the set of selected references. As input, the function requires a tPED file, a standard file format used by the PLINK software suite. For details on the format see . It also requires the file 'Corpheno', containing the references with columns in the order of population, reference and order. An example of which can be found with the package. Pheno_Pop Pheno_Data Pheno_Region Colors_Pop Colors_Region Colors_Data Order UNIQID Fam --------- ---------- ------------ ---------- ------------- ----------- ----- ------ --- Example_Pop Fedorova2013 EUR 999999 black black 140 Example_ID Example_Fam Neanderthal.SP SP EUR 999999 black black 140 AltaiNea AltaiNea Denisovan.SP SP EUR 999999 black black 140 AltaiDen DenisovaPinky A file containing the record of the major and minor alleles as described by dbSNP is needed. An example file for use with the toy data can be found as 'MinMaxFreq.rda' included with the package. Only SNPs with entries in the `MinMaxFreq.rda' used will be used in the analysis. A larger MinMaxFreq file containing SNPs used on common Illumina and Affymetrix ChIPs in addition to others is available at: Additionally the user can produce their own 'MinMaxFreq' reference file. This should be a data frame with 2 unnamed columns consisting of a column for each allele and row names corresponding to the rsID. This object must be named 'MinMaxFreq' and saved as an rda file. row.names V1 V2 ------ --- --------- rs6663840 A G rs548726 C T rs10803320 C T **calculateAMidsArith** returns a dataframe containing the genetic distance of each individual to the references. This provides the raw distance measures (Starting with the prefix C\_) and indices (Normalized values, starting with the prefix I\_). The genetic distance is computed as the Euclidean distance normalized by the number of SNPs, between each individual and all the references used. AMids for a single individual from any dataset can be computed provided there is a reasonable overlap between the set of SNPs for that individual and all the references used. The AMids can take values from 0 to 2. In our experience, the values are in the range 0.4 to 1.1. The normalized values of the distances are such that the highest reference is scored as 100, the lowest as 0 and all others adjusted accordingly. These indices place the individual in the genomic map, thus, they provide a global overview on the number of relevant references for each individual. **calculateAMids** -------------- For each individual, calculateAMids computes the genetic distances amongst that individual and the set of HGDP references (or a set provided by the user). As input, the function requires a PED formatted file. PED formatting is the standard file format required by the PLINK software suite. For details on the format see . It also requires a file containing the ids of the individuals to be used as references, and the population they correspond to. As output, returns a dataframe containing the genetic distance of each individual to the all HGDP references. We provide the raw distance measures (starting with the prefix C_) and indices (normalized Values, starting with the prefix I_). As with calculateAMidsArith, genetic distance is computed as the Euclidean distance normalized by the number of SNPs, between each individual and the references. The user can include new references in AMids by editing the file 'HGDP_References.txt', inserting the population and the corresponding invidual’s name. A demo file using 500 SNPs is included using data from the HGDP called 'HGDP_500SNPs.ped', this is compatiable with the 'CorPheno' file also included in the package. The small number of SNPs and samples will result in a highly noisy example. For more information on using **calculateAMids** see, ?**calculateAMids**() ###Producing an Input File for calculateAMids### The PED file should include individuals that will be taken as the population references, which will be used to calculate the ancestry mapper indexes (AMIds) for the user dataset. In our original work we used as references the 51 populations included in the Human Genome Diversity Project. The HGDP dataset can be obtained at . To merge a custom ped file with a ped file with the references, users could use PLINK. The commands \-\-bmerge or \-\-merge are used to merge two ped files. In most cases there will be strand inconsistencies, that can be rectified by flipping snps, using the command \-\-flip. SNPs that are CG AT are impossible to determine which strand they are in and as such be removed. Ancestry mapper requires the ped files to be in the 1/2 coding system. The individual Ids are taken as the second column of the ped file; these ids should be unique. We have produced a bed file with the references for the 51 HGDP populations, with 630,597 snps; the file is named HGDP_51RefAM_AutosomalSnps_630597_ACGT and can be obtained at . **plotAMids** --------- **plotAMids** is used to visualize the relationship amongst individuals and the references. **plotAMids** takes as input the dataframe of genetic distances returned by **calculateAMidsArith**. The user can also provide a file with phenotypes for each individual which will be visible in the plot. Colors used were taken from the BlBrewer and RedBl packages. **createMedoid** --------- **createMedoid** constructs an arithmetic reference from a tPED containing only one population with at least 10 individuals. The tPED should contain only individuals of a specific population and be formatted and prepared as described in 'Producing a Sample Input File'. **refAdd** --------- **refAdd** adds user-supplied references or data to the CorPheno file, giving population, color and ordering information to cluster the results of samples by population. Producing a Sample Input File ==================== AncestryMapper's references are produced using sequences in the same strand orientation as that used by dbSNP. Thus, all data input files for the functions **calculateAMidsArith** and **createMedoid** need also to be in the same orientation as dbSNP. If your data is already in the same strand orientation as dbSNP you can skip this section. The format for any sample data to be analysed by the function **calculateAMidsArith** or used to create a reference with **createMedoid** is a PLINK tPED. All files used in this tutorial can be found in the 'extdata' folder with the package. An example BED file, 'All-00', containing the dbSNP alleles for 1000 SNPs is included. A full version of the 'All-00' BED file covering 53,509,352 SNPs from the dbSNP database can be found at: It is faster to use PLINK BED format for the following steps; at the end the BED file needs to be converted to tPED format. It is faster to extract only those SNPs being used in your dataset from the All-00 BED and working with the output. The most efficient way to get data in the same orientation as that used by dbSNP is to to merge the 00-All file with your BED file using the \-\-merge command in PLINK; both files should be in the ACGT format. ```{r, eval=FALSE} $plink --bfile All-00 --bmerge Prep.bed Prep.bim Prep.fam --make-bed --out PrepMerge ``` If there are SNPs in different strand orientations, a list of SNPs affected will be output in a file ending in '-merge.misssnp'. Flip these SNPs using the file ending in '-merge.missnp' in your sample file before merging the two files again. ```{r, eval=FALSE} $plink --bfile Prep --flip PrepMerge-merge.missnp --make-bed --out PrepFlipped $plink --bfile All-00 --bmerge PrepFlipped.bed PrepFlipped.bim PrepFlipped.fam --make-bed --out PrepMerge ``` At the end of this merging, the 'individual' from the All-00 file should be removed using the \-\-remove command in PLINK. A file containing the ID and family ID for the All-00 individual is included as `IndRm' ```{r, eval=FALSE} $plink --bfile PrepMerge --remove IndRm --make-bed --out PrepMerge ``` SNPs that are CG or AT are invisible to the strand issue. The current AncestryMapper functions exclude them automatically from analysis. Sites with any missingness are also automatically excluded from analysis. Thus, if you wish to use any SNPs with samples missing you will need to impute replacements or exclude individuals who have missing genotypes. The final step is to transform the BED file to a tPED file and to output a list of the SNPs used. ```{r, eval=FALSE} $plink --bfile PrepMerge --transpose --recode --write-snplist --out PrepFinal ``` It is important to keep a file with the SNPs used in the sample with the ending '.snplist' as **calculateAMidsArith** will search for all the '.snplist' files corresponding to the path of the tPEDs. Example Data ======== Example Data in the form of samples and references are provided containing 1000 SNPs with data from numerous populations and datasets containing 147 populations and 591 individuals. Due to the low number of SNPs, the plots will look noisy. Additional Data ======== Full-sized medoids and a larger dbSNP reference ('00-All') for use with real user samples are currently being hosted at: Tutorial ======== Below is a short tutorial #### The first step is to call the paths to the example data files distributed with the package. ```{r, eval=FALSE} library(AncestryMapper) #Path to folder containing population references Refs <- system.file('data', package = 'AncestryMapper') #Path to folder containing samples in tPED format tpeds <- system.file('extdata', package = 'AncestryMapper') #Path to CorPheno file Corpheno <- system.file('extdata', 'CorPheno', package = 'AncestryMapper') #Path to dbSNP allele data file All00Frq <- system.file('data', 'MinMaxFreq.rda', package = 'AncestryMapper') ``` #### Calculate the genetic distance of the samples in the PLINK tPED files to the references. ```{r, eval=FALSE} genetic.distance <- calculateAMidsArith(pathTotpeds = tpeds, NameOut = 'Example', pathToAriMedoids = Refs, pathAll00 = All00Frq) ``` #### The next step is to plot the results. In this example we are using the the normalised values for each individual. (columnPlot = 'I') For more information on the function arguments see ?plotAMids() ```{r, eval=FALSE} plotAMids(AMids = genetic.distance, phenoFile = Corpheno, columnPlot = 'I') ``` #### To plot samples by population or place them in a certain order, the samples need to be added to the CorPheno file. If samples are not in the CorPheno file, the samples will be plotted under 'Undefined'. If one or more references aren't in the CorPheno, then they will be plotted at the start of the y axis before the entries with defined orders. An entry with any ID will suffice for references if the user does not wish to add or use the sample IDs from the data used to create the reference. The CorPheno can easily be modified at the text level to change things such as population orders, colors or add new entries. Alternatively the **refAdd** function may be more convenient. A tutorial on using the **refAdd**. For more information see ?refAdd() #### Users can also create their own population references from tPEDs, containing individuals from a specific population. This is done using the **createMedoid** function. The required arguments are: - pathTotpeds: The path to the tPED file(s), one tPED per population - pathAll00: The path to the file containing the dbSNP alleles Other options also exist and can be found with ?createMedoid() Example: ```{r, eval=FALSE} #Path to folder containing samples for each population in tPED format tpeds <- system.file("extdata", package = "AncestryMapper") #Path to dbSNP allele data file All00Frq <- system.file("data", "MinMaxFreq.rda", package = "AncestryMapper") createMedoid(pathTotpeds = tpeds, pathAll00 = All00Frq) ``` The resulting rda file will be named with the prefix of 'medoidArithmetic_', followed by the name of the tPED used and the number of SNPs. E.g. medoidArithmetic_Demo_1000_ChipMan.rda In order for the new reference to be used, it will need to in the path specified. For more information on the presence of 'ChipMan' at the end of the file name, see ?**createMedoid**() ####