This function takes f2_blocks as input, computes outgroup f3 statistics, and then computes PCA and MDS from the f3 statistics.
calc_mds(f2_blocks, poplabels, outgroup)Returns a data frame storing the PCA and MDS results.
#These lines assign file names to variables file_anc, file_mut, poplabels, file_map.
#see https://myersgroup.github.io/relate/getting_started.html#Output for file formats
file_anc <- system.file("sim/msprime_ad0.8_split250_1_chr1.anc.gz", package = "twigstats")
file_mut <- system.file("sim/msprime_ad0.8_split250_1_chr1.mut.gz", package = "twigstats")
#see https://myersgroup.github.io/relate/input_data.html for file formats
poplabels <- system.file("sim/msprime_ad0.8_split250_1.ind.poplabels", package = "twigstats")
file_map <- system.file("sim/genetic_map_combined_b37_chr1.txt.gz", package = "twigstats") #recombination map (three column format)
f2_blocks <- f2_blocks_from_Relate(file_anc = file_anc, file_mut = file_mut, poplabels = poplabels, file_map = file_map, t = 1000)
df <- calc_mds(f2_blocks, poplabels, "P4")
print(head(df))
#> PC1 PC2 PC3 ID POP GROUP SEX method
#> 1 -4.925615 -3.092605 0.3954621 tsk_0 tsk_0 P1 NA PCA
#> 2 -5.251956 -2.649876 1.0559442 tsk_1 tsk_1 P1 NA PCA
#> 11 -3.404097 2.342448 0.1492084 tsk_10 tsk_10 P2 NA PCA
#> 12 -3.528335 4.113717 0.5620541 tsk_11 tsk_11 P2 NA PCA
#> 13 -4.202899 1.147829 -0.9252663 tsk_12 tsk_12 P2 NA PCA
#> 14 -3.734241 0.796142 -0.1207356 tsk_13 tsk_13 P2 NA PCA