Analyze Data

Prepare and import genetic data

assignPOP accepts genetic data in GENEPOP¹ and STRUCTURE² formats. Both formats can be either diploid or haploid data.

When importing a GENEPOP file, two forms of GENEPOP can be read into R using the function read.Genepop(). Populatino names can be specified using the argument pop.names, and the order of population names should follow the group order in your GENEPOP file³. These population names can be used to customize the assignment accuracy plot later on. Use the argument haploid to specify your data type (e.g., use haploid=T for haploid data).

YourGenepop <- read.Genepop( "simGenepop.txt", pop.names=c("pop_A","pop_B","pop_C"), haploid = FALSE)

The example input file, simGenepop.txt, contains three simulated populations (or subpopulations) and each population has 30 individuals by 1,000 SNP loci. Depending on the data size, the importing process may take a few secnods to minutes. A progress bar will be running while importing the data. After importing the file, you should see the following message printed in your R console. Use the information to double check whether your data were correctly imported.

## ############### assignPOP v1.x.x ###############
## 
## A GENEPOP format file was successfully imported!
## 
## Imported Data Info: 90 obs. by 1000 loci (diploid)
## Number of pop:3
## Number of inds (pop_A): 30
## Number of inds (pop_B): 30
## Number of inds (pop_C): 30
## 
## DataMatrix: 90 rows by 2001 columns, with 2000 allele variables
## 
## Data output in a list comprising the following three elements:
## YOUR_LIST_NAME$DataMatrix
## YOUR_LIST_NAME$SampleID
## YOUR_LIST_NAME$LocusName

When importing a STRUCTURE file, use the function read.Structure(), and specify the data type using the argument haploid. Your Structure file should include locus names (each name separated by whitespace or tab) in the first row as the column name (header). The first and second column should be your sample ID and population label, respectively. No whitespace should be used within the sample ID and population label. The rest of columns are genetic data. It does not matter whether your file has column names for sample ID and/or population label.

YourStructure <- read.Structure( "simStructure.txt", haploid = TRUE)

After importing your STRUCTURE file, you will also see the printed message in the console. For the sake of simplicity, we only use the GENEPOP dataset to demonstrate the rest of analyses.

Missing data and homozygosity. If an locus has missing alleles or is homomzygous across all individuals in the dataset, the locus will be automatically removed.

Remove low variance loci (optional)

When analyzing genetic data, you can remove low variance loci⁴ across the dataset for further analyses. The default of variance threshold is 0.95, meaning that a locus will be removed from the dataset if its major allele occurs in over 95% of individuals across the populations.

To remove low variance loci, use the following function and provide a new return object name (YourGenepopRd in this case) :

YourGenepopRd <- reduce.allele( YourGenepop, p = 0.95)

After entering the code, R console should print the following message.

## New data matrix has created! :)
## New DataMatrix size: 90 rows by 1387 columns
## 614 columns (alleles) have been removed
## 693 loci remain

In our example, 307 low variance loci⁵ were removed, leaving 693 loci in the new data set, which will be used in further analyses.

Concatenate genetic and non-genetic data

Concatenating genetic and non-genetic data is one of the novel features in this package. After importing your GENEPOP file, you can use the function compile.data() to concatenate a non-genetic dataset and the genetic matrix returned from either read.Genepop(), read.Structure(), or reduce.allele(). Your non-genetic data should be saved in a .csv file (elements separated by commas) or table-like text file (elements separated by spaces) in which the first column must be sample IDs that match the IDs in your GENEPOP or STRUCTURE file, and the rest of columns are non-genetic data, whether it is numeric or categorical. If a sample ID exists in only one of the data types, the individual will be ignored in the new integrated data.

To concatenate datasets, use the following function and provide a new return object name (YourIntegrateData in this case):

YourIntegrateData <- compile.data( YourGenepopRd, "morphData.csv")

The above morphData.csv file has the same individuals as the simGenepop.txt and it contains 4 numeric variables (morphometric measurements) for each individual. After executing the function, it will prompt the following message and wait for your answer to verify the data type (numeric or categorical).

## Import a .CSV file.
## 4 additional variables detected.
## Checking variable data type...
##  D1.2(integer)  D2.3(integer)  D3.4(integer)  D1.4(integer)
##  Are they correct? (enter Y/N):

If a variable is categorical data, it will automatically convert it to new dummy variables using model.matrix() so that the categorical data can be quantified for further analyses. If your data type was incorrectly identified, simply enter N and then follow the prompted dialogue to change your data type. If your data type is correctly identified (entry is Y in this example), then it will print the following message.

## New data set created!!
## It has 90 observations by 1391 variables
## including 693 loci(1386 alleles) plus 4 additional variables(4 columns)

Now the integrated data (693 SNP loci plus 4 morphometric measurements) matrix is saved in the return object named comin.

Prepare and import non-genetic data

In addition to analyzing genetic-only or integrated data, you can perform the same analysis on non-genetic data so that you could compare the assignment results between data types. A non-genetic data set (collected from source populations) should include 1) sample IDs in the first column, 2) population label⁶ in the last column, and 3) non-genetic data in columns between the sample ID and population label columns. The file should be saved in a .csv (comma delimited) or table-like text file that can be read into R using R basic functions read.csv() or read.table(). For example, if you are going to perform assignment test using the sample morphometric data (morphData.csv) provided in the pakcage, you would need to add a population column after the existing last column. It can be edited manually (e.g., using Excel) or by R scripts (see below). If your sample IDs and population names are numeric, make sure to convert them to factor data type using as.factor() function.

#Import morphometric data containing sample IDs and features
morphdf <- read.csv( "morphData.csv", header=TRUE )
#Create a string vector for population label (repeat each name for 30 individuals)
pop_label <- c( rep("pop_A", 30), rep("pop_B", 30), rep("pop_C", 30) ) 
#Add the pop_label to the last column; 'morphdf_pop' is a data frame with population label in the last column
morphdf_pop <- cbind(morphdf, pop_label)

#Handling numeric population name
#A set of population names: 1, 2, and 3
pop_label <- c( rep(1, 30), rep(2, 30), rep(3, 30))
#Add the pop_label
morphdf_pop <- cbind(morphdf, pop_label)
#Convert population name to factor data type 
morphdf_pop$pop_label <- as.factor(morphdf_pop$pop_label)

Now the non-genetic (morphometric) data with its population label is saved in the object morphdf_pop.

Perform resampling cross-validation

We provide two resampling cross-validation methods, Monte-Carlo and K-fold cross-validations, to evaluate baseline data. Monte-Carlo cross-validation helps estimate the mean and variance of assignment accuracy through resampling random training individuals; K-fold cross-validation helps determine membership probability across all individuals through using one group as test individuals and the remaining K-1 groups as training individuals (hence every individual is tested once).

When analyzing genetic data (genetic-only or integrated), users can specify what proportions of loci should be used as training loci, and determine whether the proportion of loci is chosen either randomly or based on locus F_ST estimated within training individuals. This helps evaluate whether using a subset of loci performs as well as using all loci. Most informative loci also can be identified via the function check.loci (see Identify informative loci)

assign.MC( YourGenepopRd, train.inds=c(0.5, 0.7, 0.9), train.loci=c(0.1, 0.25, 0.5, 1),
           loci.sample="fst", iterations=30, model="svm", dir="Result-folder/")

assign.kfold( YourGenepopRd, k.fold=c(3, 4, 5), train.loci=c(0.1, 0.25, 0.5, 1), 
             loci.sample="random", model="lda", dir="Result-folder2/")

Calculate assignment accuracy

When resampling cross-validations are done, you can calculate assignment accuracies using the following functions.

accuMC <- accuracy.MC(dir = "Result-folder/") #Use this function for Monte-Carlo cross-validation results
accuKF <- accuracy.kfold(dir = "Result-folder2/") #Use this function for K-fold cross-validation results

The functions read through each assignment result in the designated folder (arg. dir) and calculate the assignment accuracies for overall and individual populations. The results will be saved in a text file named Rate_of…txt, and a return object (accuMC and accuKF in the above examples). The output text file also can be read into R using the basic function read.table() as follows.

accuMC <- read.table("Rate_of....txt", header=T)

Compared with K-fold cross-validation, the results of Monte-Carlo cross-validation, with reasonable iterations (e.g., resample > 30 times), are more suitable for estimating mean and variance of assignment accuracy. As such, using the function accuracy.MC() for Monte-Carlo results is more useful than using accuracy.kfold(). The results generated by accuracy.MC() (or accuracy.kfold()) can be visualized in an assignment accuracy plot; the results generated by the K-fold cross-validation (assign.kfold()) can be directly used to make a membership probability plot. See Visualize Results page for more details.

Identify informative loci

In some cases, using a subset of high F_ST loci may produce similar assignment accuracy with using all available loci. Identification of these loci not only could help reduce time and cost in preparing samples in the future but also could help identify loci that might be associated with functional genes.

The function check.loci() reads through each training locus file associated with the assignment test, counts the frequency of occurrence of each locus, and writes the results to a text file.

check.loci(dir = "Result-folder/", top.loci = 20)

This function should only be used for the cross-validation results based on resampling high F_ST training loci (loci.sample = "fst") in the assign.MC() analysis.

By default, this function outputs top 20 loci - most frequently used high F_ST training loci across your assignment tests, but an arbitrary number can be specified in the argument top.loci.

When the function is executed, as shown below, you will be prompted to choose which assignment results you would like to check. The options are the groups of different proportions or numbers of training individuals that you specified in the argument train.inds of the function assign.MC().

## 3 levels of training individuals are found.
## Which levels would you like to check? (separate levels by a whitespace if multiple)
## Options: 0.5, 0.7, 0.9, or all
## enter here: (You will enter your answer here)

Because a locus F_ST value estimated based on a training set likely varies among subsets of training individuals, we recommend running the function check.loci() across all levels of proportions or numbers of training individuals to evaluate if the loci consistently are those high F_ST loci. If so, then they could be recognized as truly informative loci (high F_ST is not due to sampling artifacts).

The output file (High_Fst_Locus_Freq.txt) includes a list of locus names ordered by F_ST value (highest F_ST in the top row). The number in parentheses right after each locus name indicates the number of tests that locus appeared in that rank. A snapshot of output content is shown below.

Loci occur in top 20 high Fst across all training data
top.1(1): Locus_171(360), 
top.2(1): Locus_442(360), 
top.3(1): Locus_475(360), 
top.4(1): Locus_481(360), 
top.5(1): Locus_696(360), 
top.6(1): Locus_697(360), 
top.7(1): Locus_729(360), 
top.8(1): Locus_745(360), 
top.9(1): Locus_812(360), 
top.10(1): Locus_941(360), 
top.11(1): Locus_992(360), 
top.12(5): Locus_113(248), Locus_114(76), Locus_115(24), Locus_320(8), Locus_139(4), 
top.13(7): Locus_245(240), Locus_181(72), Locus_147(24), Locus_115(8), Locus_560(8), Locus_137(4), Locus_160(4), 

The above results show that the top.1-top.11 loci are the top 11 highest F_ST loci are across all 360 tests. Five loci (top.12(5): Locus_113, Locus_114, Locus_115, Locus_320, and Locus_139) appear to be the 12th highest F_ST locus across the tests, and Locus_113 occurs most frequently (248 out of 360 tests) in this rank.

Prepare and import data of unknown individuals

To predict sources of unknown individuals, users should save the data of unknown individuals in another GENEPOP or STRUCTURE file, and then import it into R using the same functions as we used for the baseline data.

For the GENEPOP format, data should be saved with only one group - only one ‘pop’ label in the file (e.g., simGenepopX.txt). For the STRUCTURE format, the column of population label should also be included, but the label can be anything (this column will be ignored when performing the assignment test).

#Import a GENEPOP file containing unknown individuals, and ignore the argument `pop.names` 
YourGenepopUnknown <- read.Genepop( "simGenepopX.txt" )

#Import a STRUCTURE file containing unknown individuals.
YourStructureUnknown <- read.Structure( "simStructureX.txt" )

You can also concatenate genetic and non-genetic data of unknown individuals using the function compile.data as follows.

YourIntegrateUnknown <- compile.data(YourGenepopUnknown, "morphDataX.csv")

The non-genetic (morphDataX.csv) and genetic (simGenepopX.txt) data should have identical sample IDs. In your non-genetic data, sample IDs should be saved in the first column, and there would be no population label in the last column.

If you are analyzing non-genetic data alone, use R basic functions, read.csv() or read.table(), to import your file.

#Import non-genetic data containing sample IDs and features
NonGeneticUnknown <- read.csv( "morphDataX.csv", header=TRUE )

Perform assignment test on unknown individuals

After importing your baseline data and unknown individuals (test data), the assignment test can be performed using the function assign.X. Make sure your training and test data have the same data type and feature names. Use the arguments x1 and x2 to specify the baseline and test data,respectively.

# 1.Perform assignment test using genetic data and naive Bayes
assign.X( x1=YourGenepopRd, x2=YourGenepopUnknown, dir="Result-folder3/", model="naiveBayes")

# 2.Perform assignment test using integrated data and decision tree
assign.X( x1=YourIntegrateData, x2=YourIntegrateUnknown, dir="Result-folder4/", model="tree")

# 3.Perform assignment test uisng non-genetic data and random forest
assign.X( x1=morphdf_pop, x2=OtherUnknown, dir="Result-folder5/", model="randomForest")

# Use `?assign.X` to see other argument settings.

The assignment result will be saved in a text file named AssignmentResult.txt under the designated folder. This file includes sample IDs of your unknown individuals, predicted populations and the probabilities. When the assignment is done, it will prompt a message asking you whether you want to make a membership probability plot. Simply enter Y to visualize the results right away.

Handling unmatched alleles between baseline and test data. In some cases, alleles found in test data may not exist in baseline data, or vice versa. The function assign.X will check if alleles of a locus are identical between your baseline and test data. If they are different, it will ask you if you want to use the loci in common to continue the analysis. If answer ‘Y’, the assignment test will ignore the loci that have different alleles between baseline and test data.

PCA for data transformation

The assignment functions, assign.MC(), assign.kfold() and assign.X(), use PCA for dimensionality reduction. It converts your independent variables to Principle Components and retain the PCs that have eigenvalues greater than 1 (the Kaiser-Guttman criterion) as new features (arg. pca.PCs="kaiser-guttman"). Or, you can specify an integer in the argument pca.PCs (e.g., pca.PCs = 10) to choose a specific number of PCs to be retained.

The PCA will always perform on the genetic data, whether you are analyzing genetic-only or integrated data. When analyzing integrated data, you have three options to perform PCA on the data. First, you can perform PCA across all features (set arg. pca.method="mixed") - each PC includes information of genetic and non-genetic data. This is default setting. Second, you can perform PCA on genetic and non-genetic data independently (set arg. pca.method="independent") - transformed data includes genetic PCs and non-genetic PCs. Third, you can choose to use original non-genetic data as features without performing PCA on them (set. arg. pca.method="original"). This option uses genetic PCs and original non-genetic data to build predictive models.

When analyzing non-genetic data alone, you can determine whether PCA is performed on the data or not. You set the argument pca.method=TRUE to perform PCA, or pca.method=FALSE to not perform PCA.

Data standardization

When analyzing non-genetic or integrated data, you have options to standardize (zero mean and unit variance) the independent variables or to keep the original data. If non-genetic data has varying ranges of values and is not standardized, it may be problematic when one feature has a variance that is orders of magnitude larger than others, as it might make the classifier unable to learn from other features correctly. To center and scale the data, set the argument scaled=TRUE in the function assign.MC(), assign.kfold() or assign.X(). If you are analyzing genetic data, then standardization is not required because all the loci (alleles) have been converted to a binary-like value 0, 0.5, or 1.