Autogenous - FST estimates
Luciano V Cosme
2026-04-06
Population label key: The code labels used in the analysis below correspond to the following manuscript labels: AUT = AUTO, MAN = NON-AUTO, NEW = NON-AUTO-FIELD.
Libraries
## Warning in system("timedatectl", intern = TRUE): running command 'timedatectl'
## had status 1 and error message 'Function not implemented'library(here)
library(adegenet)
library(StAMPP)
library(reshape2)
library(data.table)
library(scales)We can use different data sets to run our fst estimates.
Clean env and memory
# Remove all objects from the environment
rm(list = ls())
# Run the garbage collector to free up memory
gc()## used (Mb) gc trigger (Mb) max used (Mb)
## Ncells 2558351 136.7 4255384 227.3 3629285 193.9
## Vcells 4281062 32.7 10146329 77.5 8388606 64.01. Estimate pairwise FST for all populations
We can convert to raw format
plink \
--keep-allele-order \
--bfile output/ldna/files/file1 \
--recodeA \
--out output/fst/pop_pairwise \
--silent;
grep 'samples\|variants\|remaining' output/fst/pop_pairwise.log## 110353 variants loaded from .bim file.
## 110353 variants and 60 people pass filters and QC.
## Among remaining phenotypes, 28 are cases and 22 are controls. (10 phenotypesLook at https://rdrr.io/cran/StAMPP/man/stamppFst.html for details of Fst estimations
pop_pairwise <-
read.PLINK(
here(
"output", "fst", "pop_pairwise.raw"
),
quiet = FALSE,
chunkSize = 1000,
parallel = require("parallel"),
n.cores = 4
)
summary(pop_pairwise)Now lets convert the genlight object to Stampp format, and estimate pairwise Fst values
The command below would also work, but you can simplify it and put only the numbers: pop_pairwise2 <- stamppFst(pop_pairwise, nboots=10, percent=95 + nclusters=8)
This chunk will take a couple minutes to run.
# Convert
pop_pairwise2 <- stamppConvert(pop_pairwise, type="genlight")
# 10 bootstraps, 95% CI, 8 threads
pop_pairwise2 <- stamppFst(pop_pairwise2, 10, 95, 8)Save the result
Now lets look at the object
## Length Class Mode
## Fsts 9 -none- numeric
## Pvalues 9 -none- numeric
## Bootstraps 16 data.frame listIf you want you can save the fst values as csv.
# Convert to data frame
pop_pairwise2_df <- data.frame(pop_pairwise2)
# Save it
write.csv(pop_pairwise2_df, file = here("output", "fst", "pop_pairwise2_df.csv"))Check the Fst values
## Fsts.MAN Fsts.AUT Fsts.NEW Pvalues.MAN Pvalues.AUT Pvalues.NEW
## MAN NA NA NA NA NA NA
## AUT 0.12987132 NA NA 0 NA NA
## NEW 0.04311201 0.1273737 NA 0 0 NA
## Bootstraps.Population1 Bootstraps.Population2 Bootstraps.1 Bootstraps.2
## MAN MAN AUT 0.12895659 0.12921884
## AUT MAN NEW 0.04237085 0.04259314
## NEW AUT NEW 0.12638433 0.12664232
## Bootstraps.3 Bootstraps.4 Bootstraps.5 Bootstraps.6 Bootstraps.7
## MAN 0.12927785 0.12934050 0.12957487 0.12989073 0.13012637
## AUT 0.04264201 0.04290901 0.04311632 0.04313661 0.04321334
## NEW 0.12690892 0.12697258 0.12707649 0.12711577 0.12746235
## Bootstraps.8 Bootstraps.9 Bootstraps.10 Bootstraps.Lower.bound.CI.limit
## MAN 0.1303125 0.13053896 0.13076437 0.12895659
## AUT 0.0432618 0.04328578 0.04364687 0.04237085
## NEW 0.1275964 0.12792371 0.12810409 0.12638433
## Bootstraps.Upper.bound.CI.limit Bootstraps.p.value Bootstraps.Fst
## MAN 0.13053896 0 0.12987132
## AUT 0.04328578 0 0.04311201
## NEW 0.12792371 0 0.12737374We can subset the object
# Create an object for columns with 'Fsts.'
fst_columns <- grep("Fsts\\.", names(pop_pairwise2_df), value = TRUE)
fst_df <- pop_pairwise2_df[, fst_columns]
# Create an object for columns with 'Pvalues.'
pvalue_columns <- grep("Pvalues\\.", names(pop_pairwise2_df), value = TRUE)
pvalue_df <- pop_pairwise2_df[, pvalue_columns]
# You can view the first few rows of these dataframes to confirm
head(fst_df)## Fsts.MAN Fsts.AUT Fsts.NEW
## MAN NA NA NA
## AUT 0.12987132 NA NA
## NEW 0.04311201 0.1273737 NA## Pvalues.MAN Pvalues.AUT Pvalues.NEW
## MAN NA NA NA
## AUT 0 NA NA
## NEW 0 0 NAAll the p-values are zero, indicating statistical significance
We will convert the data into a matrix.
# First we can rename the columns of the fst values object
# Use gsub to remove 'Fsts.' from the column names
names(fst_df) <- gsub("Fsts\\.", "", names(fst_df))
# Create matrix
aa <- as.matrix(fst_df)
aa[upper.tri(aa)] <- t(aa)[upper.tri(t(aa))]
head(aa)## MAN AUT NEW
## MAN NA 0.1298713 0.04311201
## AUT 0.12987132 NA 0.12737374
## NEW 0.04311201 0.1273737 NALets check if the matrix is symmetric.
## [1] TRUEWe will also add NA on the upper left side of the matrix.
Now we have to convert the matrix to a data frame to plot it with ggplot.
## Warning: The melt generic in data.table has been passed a matrix and will
## attempt to redirect to the relevant reshape2 method; please note that reshape2
## is superseded and is no longer actively developed, and this redirection is now
## deprecated. To continue using melt methods from reshape2 while both libraries
## are attached, e.g. melt.list, you can prepend the namespace, i.e.
## reshape2::melt(aa). In the next version, this warning will become an error.## Var1 Var2 value
## MAN:3 MAN:3 Min. :0.04311
## AUT:3 AUT:3 1st Qu.:0.08524
## NEW:3 NEW:3 Median :0.12737
## Mean :0.10012
## 3rd Qu.:0.12862
## Max. :0.12987
## NA's :6Now lets plot the data with ggplot. You can click in the little square on the top left of the plot to open it on a new window. It will have the right proportions.
pairfst.f <- ggplot(pairfst.long, aes(Var1, Var2)) +
geom_tile(aes(fill = value), colour = "white") +
scale_fill_gradient(
low = "white",
high = "#71b6ff",
name = "Fst",
na.value = "white",
limits = c(0, 0.5)
) +
scale_x_discrete(position = "top") +
theme_bw() +
geom_text(aes(label = ifelse(
is.na(value), "", formatC(value, digits = 2, format = "f")
)), size = 6) +
theme(
axis.text.x = element_text(angle = 90, hjust = 1, size = 20),
axis.title = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_blank(),
panel.background = element_blank(),
axis.text.y = element_text(hjust = 0, size = 20)
)
pairfst.f
While the Fst between the wild populations, NEW and MAN, is 0.04, the
values between them the the autogenous is 0.13.
2. Estimate Fst between males and females
Note: This section (sex-stratified pairwise Fst) was not used in the final version of the manuscript.
We do not have sex information for mosquitoes of MAN, but we do have it for NEW and AUT
We can import the fam file to identify males and females
# Specify the path using the 'here' function
fam_file_path <- here("output", "ldna", "files", "file1.fam")
# Now, read the file using the path - using data.table library
fam_data <- fread(fam_file_path, header = FALSE)
# Rename columns for clarity
setnames(fam_data, c("FID", "IID", "PID", "MID", "Sex", "Phenotype"))
head(fam_data)## FID IID PID MID Sex Phenotype
## <char> <int> <int> <int> <int> <int>
## 1: MAN 306 0 0 0 -9
## 2: MAN 307 0 0 0 -9
## 3: MAN 308 0 0 0 -9
## 4: MAN 309 0 0 0 -9
## 5: MAN 310 0 0 0 -9
## 6: MAN 311 0 0 0 -9Separate males and females per population
# Create a list of data frames, each representing a different FID
fam_list <- split(fam_data, fam_data$FID)
# Initialize lists to store male and female data separately for each FID
males_per_fid <- list()
females_per_fid <- list()
# Loop through each FID and subset males and females
for(fid in names(fam_list)) {
males_per_fid[[fid]] <- fam_list[[fid]][Sex == 1, ]
females_per_fid[[fid]] <- fam_list[[fid]][Sex == 2, ]
}Now we can save a list of mosquitoes for which we know the sex information
# Combine male and female data into one data frame
combined_data <- rbindlist(c(males_per_fid, females_per_fid), use.names = TRUE, fill = TRUE)
# Select only the FID and IID columns
plink_data <- combined_data[, .(FID, IID)]
# Export the combined data to a text file
fwrite(plink_data, file = here("output", "fst", "male_female.txt"), sep = "\t", quote = FALSE)Now we can subset the data
plink \
--keep-allele-order \
--bfile output/ldna/files/file1 \
--keep output/fst/male_female.txt \
--out output/fst/new_aut \
--make-bed \
--silent;
grep 'samples\|variants\|remaining' output/fst/new_aut.log## 110353 variants loaded from .bim file.
## --keep: 50 people remaining.
## Total genotyping rate in remaining samples is 0.976527.
## 110353 variants and 50 people pass filters and QC.
## Among remaining phenotypes, 28 are cases and 22 are controls.We can import the fam file to identify males and females
# Specify the path using the 'here' function
fam_file_path <- here("output", "fst", "new_aut.fam")
# Now, read the file using the path - using data.table library
fam_data2 <- fread(fam_file_path, header = FALSE)
# Rename columns for clarity
setnames(fam_data2, c("FID", "IID", "PID", "MID", "Sex", "Phenotype"))
head(fam_data2)## FID IID PID MID Sex Phenotype
## <char> <int> <int> <int> <int> <int>
## 1: AUT 384 0 0 1 2
## 2: AUT 385 0 0 1 2
## 3: AUT 386 0 0 1 2
## 4: AUT 387 0 0 1 2
## 5: AUT 388 0 0 1 2
## 6: AUT 389 0 0 1 2Now we can change the FID to separate the males and females for each population, so we end up with 4 groups of mosquitoes
# Add the suffix '_M' or '_F' to the IID column based on the Sex column
fam_data2$FID <-
ifelse(
fam_data2$Sex == 1,
paste0(fam_data2$FID, "_M"),
ifelse(fam_data2$Sex == 2, paste0(fam_data2$FID, "_F"), fam_data2$FID)
)
# Now you can check the first few rows to confirm the change
head(fam_data2)## FID IID PID MID Sex Phenotype
## <char> <int> <int> <int> <int> <int>
## 1: AUT_M 384 0 0 1 2
## 2: AUT_M 385 0 0 1 2
## 3: AUT_M 386 0 0 1 2
## 4: AUT_M 387 0 0 1 2
## 5: AUT_M 388 0 0 1 2
## 6: AUT_M 389 0 0 1 2Now we can save the file and replace the original .fam file
# Save the modified dataframe to a text file
fwrite(fam_data2, file = here("output", "fst", "new_aut.fam"), sep = "\t", quote = FALSE, col.names = FALSE)Check it
## AUT_M 384 0 0 1 2
## AUT_M 385 0 0 1 2
## AUT_M 386 0 0 1 2
## AUT_M 387 0 0 1 2
## AUT_M 388 0 0 1 2
## AUT_M 389 0 0 1 2
## AUT_M 390 0 0 1 2
## AUT_M 391 0 0 1 2
## AUT_M 392 0 0 1 2
## AUT_M 393 0 0 1 2Now we can convert it to raw format and repeat the fst estimates
We can convert to raw format
# For consistance we need to use the same thresholds we used to create the file with the 3 populations. So, we set genotype missingness to 20% and MAF to 5%
plink \
--keep-allele-order \
--bfile output/fst/new_aut \
--recodeA \
--geno 0.2 \
--maf 0.05 \
--out output/fst/sex \
--silent;
grep 'samples\|variants\|remaining' output/fst/sex.log## 110353 variants loaded from .bim file.
## 357 variants removed due to missing genotype data (--geno).
## 12221 variants removed due to minor allele threshold(s)
## 97775 variants and 50 people pass filters and QC.
## Among remaining phenotypes, 28 are cases and 22 are controls.Clean env and memory
# Remove all objects from the environment
rm(list = ls())
# Run the garbage collector to free up memory
gc()## used (Mb) gc trigger (Mb) max used (Mb)
## Ncells 3153258 168.5 8277314 442.1 8277314 442.1
## Vcells 5340422 40.8 42217756 322.1 52772192 402.7sex_pairwise <-
read.PLINK(
here(
"output", "fst", "sex.raw"
),
quiet = FALSE,
chunkSize = 1000,
parallel = require("parallel"),
n.cores = 4
)
summary(sex_pairwise)Now lets convert the genlight object to Stampp format, and estimate pairwise Fst values
The command below would also work, but you can simplify it and put only the numbers: pop_pairwise2 <- stamppFst(pop_pairwise, nboots=10, percent=95 + nclusters=8)
This chunk will take a couple minutes to run.
# Convert
sex_pairwise2 <- stamppConvert(sex_pairwise, type="genlight")
# 10 bootstraps, 95% CI, 8 threads
sex_pairwise2 <- stamppFst(sex_pairwise2, 10, 95, 8)Save the result
Load the result
Now lets look at the object
## Length Class Mode
## Fsts 16 -none- numeric
## Pvalues 16 -none- numeric
## Bootstraps 16 data.frame listIf you want you can save the fst values as csv.
# Convert the 'Fsts' and 'Pvalues' matrices to dataframes
fsts_df <- as.data.frame(sex_pairwise2$Fsts)
pvalues_df <- as.data.frame(sex_pairwise2$Pvalues)
head(fsts_df)## AUT_M AUT_F NEW_M NEW_F
## AUT_M NA NA NA NA
## AUT_F 0.03971471 NA NA NA
## NEW_M 0.16839622 0.1541147 NA NA
## NEW_F 0.13607526 0.1224449 0.03188417 NACheck the p-values
## AUT_M AUT_F NEW_M NEW_F
## AUT_M NA NA NA NA
## AUT_F 0 NA NA NA
## NEW_M 0 0 NA NA
## NEW_F 0 0 0 NAAll the p-values are zero, indicating statistical significance
We will convert the data into a matrix.
# First we can rename the columns of the fst values object
# Use gsub to remove 'Fsts.' from the column names
names(fsts_df) <- gsub("Fsts\\.", "", names(fsts_df))
# Create matrix
aa <- as.matrix(fsts_df)
aa[upper.tri(aa)] <- t(aa)[upper.tri(t(aa))]
head(aa)## AUT_M AUT_F NEW_M NEW_F
## AUT_M NA 0.03971471 0.16839622 0.13607526
## AUT_F 0.03971471 NA 0.15411469 0.12244490
## NEW_M 0.16839622 0.15411469 NA 0.03188417
## NEW_F 0.13607526 0.12244490 0.03188417 NALets check if the matrix is symmetric.
## [1] TRUEWe will also add NA on the upper left side of the matrix.
Now we have to convert the matrix to a data frame to plot it with ggplot.
## Var1 Var2 value
## AUT_M:4 AUT_M:4 Min. :0.03188
## AUT_F:4 AUT_F:4 1st Qu.:0.06040
## NEW_M:4 NEW_M:4 Median :0.12926
## NEW_F:4 NEW_F:4 Mean :0.10877
## 3rd Qu.:0.14960
## Max. :0.16840
## NA's :10Now lets plot the data with ggplot. You can click in the little square on the top left of the plot to open it on a new window. It will have the right proportions.
pairfst.f <- ggplot(pairfst.long, aes(Var1, Var2)) +
geom_tile(aes(fill = value), colour = "white") +
scale_fill_gradient(
low = "white",
high = "#71b6ff",
name = "Fst",
na.value = "white",
limits = c(0, 0.5)
) +
scale_x_discrete(position = "top") +
theme_bw() +
geom_text(aes(label = ifelse(
is.na(value), "", formatC(value, digits = 2, format = "f")
)), size = 6) +
theme(
axis.text.x = element_text(angle = 90, hjust = 1, size = 20),
axis.title = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_blank(),
panel.background = element_blank(),
axis.text.y = element_text(hjust = 0, size = 20)
)
pairfst.f
The highest differentiation is between the males across population. For example, NEW males and AUT males have an Fst value of 0.17 while the females have 0.12. The differences within sex is low, 0.04 for AUT and 0.03 for NEW. The Fst estimates for females between the populations is lower than males, for females of AUT and NEW the Fst value is 0.12, while between males is 0.17. Now, between sex, for AUT males vs NEW females the fst value is 0.14, while AUT females vs NEW males is 0.15, which are quite similar.
output_path <- here("output", "fst", "figures", "sex_estimates.pdf")
ggsave(output_path, pairfst.f, height = 6, width = 6, dpi = 300)Now we can compare the males vs females using non-overlapping sliding window. First we compare the two populations, NEW and MAN vs AUT. Then, we compare the males and females for each population we have sex information (NEW and AUT), next we can compare the males of AUT vs the males of NEW, and do the estimates for the females as well. We can use vcftools to perform the estimates for us.
3. Sliding windows Fst estimates between the populations
Clean env and memory
# Remove all objects from the environment
rm(list = ls())
# Run the garbage collector to free up memory
gc()## used (Mb) gc trigger (Mb) max used (Mb)
## Ncells 3154795 168.5 8277314 442.1 8277314 442.1
## Vcells 5240036 40.0 33774205 257.7 52772192 402.7We will compare NEW vs MAN, NEW vs AUT, MAN vs AUT, then we can plot the estimates using a facet plot and add a smooth line to see if there are genomic regions with high differenciation
First we need to create a vcf with Plink. We can use the original file since we will compare all populations.
Convert to vcf
plink \
--keep-allele-order \
--bfile output/ldna/files/file1 \
--recode vcf \
--out output/fst/albo \
--silent;
grep 'samples\|variants\|remaining' output/fst/albo.log## 110353 variants loaded from .bim file.
## 110353 variants and 60 people pass filters and QC.
## Among remaining phenotypes, 28 are cases and 22 are controls. (10 phenotypesWe can use bcftools to check the sample names
## MAN_306
## MAN_307
## MAN_308
## MAN_309
## MAN_310
## MAN_311
## MAN_312
## MAN_313
## MAN_314
## MAN_315
## AUT_384
## AUT_385
## AUT_386
## AUT_387
## AUT_388
## AUT_389
## AUT_390
## AUT_391
## AUT_392
## AUT_393
## AUT_394
## AUT_395
## AUT_396
## AUT_397
## AUT_398
## AUT_400
## AUT_401
## AUT_402
## AUT_403
## AUT_404
## AUT_405
## AUT_406
## AUT_407
## AUT_408
## AUT_409
## AUT_410
## AUT_411
## AUT_412
## NEW_413
## NEW_414
## NEW_415
## NEW_416
## NEW_417
## NEW_418
## NEW_419
## NEW_420
## NEW_421
## NEW_422
## NEW_423
## NEW_424
## NEW_425
## NEW_426
## NEW_427
## NEW_428
## NEW_429
## NEW_430
## NEW_431
## NEW_432
## NEW_433
## NEW_434Now we create a list of mosquitoes for each population
# Define the VCF file path
VCF_FILE="output/fst/albo.vcf"
# Extract the list of individuals
INDIVIDUALS=$(bcftools query -l $VCF_FILE)
# Filter and create files for each population
echo "$INDIVIDUALS" | grep '^MAN_' > output/fst/MAN.txt
echo "$INDIVIDUALS" | grep '^AUT_' > output/fst/AUT.txt
echo "$INDIVIDUALS" | grep '^NEW_' > output/fst/NEW.txtCheck one file to see if is correct
## MAN_306
## MAN_307
## MAN_308
## MAN_309
## MAN_310
## MAN_311
## MAN_312
## MAN_313
## MAN_314
## MAN_315Now we can estimate Fst between the populations using vcftools. We can use overlapping windows for finer resolution
MAN vs NEW
# Since we are using the SNP chip data, we can use a window of 1Mb and 100kb steps
vcftools --vcf output/fst/albo.vcf \
--weir-fst-pop output/fst/MAN.txt \
--weir-fst-pop output/fst/NEW.txt \
--fst-window-size 1000000 \
--fst-window-step 100000 \
--out output/fst/man_new##
## VCFtools - 0.1.17
## (C) Adam Auton and Anthony Marcketta 2009
##
## Parameters as interpreted:
## --vcf output/fst/albo.vcf
## --fst-window-size 1000000
## --fst-window-step 100000
## --weir-fst-pop output/fst/MAN.txt
## --weir-fst-pop output/fst/NEW.txt
## --keep output/fst/MAN.txt
## --keep output/fst/NEW.txt
## --out output/fst/man_new
##
## Warning: Expected at least 2 parts in INFO entry: ID=PR,Number=0,Type=Flag,Description="Provisional reference allele, may not be based on real reference genome">
## Keeping individuals in 'keep' list
## After filtering, kept 32 out of 60 Individuals
## Outputting Windowed Weir and Cockerham Fst estimates.
## Weir and Cockerham mean Fst estimate: 0.034658
## Weir and Cockerham weighted Fst estimate: 0.043112
## After filtering, kept 110353 out of a possible 110353 Sites
## Run Time = 1.00 secondsCheck the result
## CHROM BIN_START BIN_END N_VARIANTS WEIGHTED_FST MEAN_FST
## 1 1 1000000 22 -0.00301908 -0.00497006
## 1 100001 1100000 21 -0.00274564 -0.00477669
## 1 200001 1200000 28 0.00575673 -0.000393366
## 1 300001 1300000 36 0.0253471 0.0213285
## 1 400001 1400000 45 0.0529739 0.0519991
## 1 500001 1500000 42 0.0606128 0.0570881
## 1 600001 1600000 49 0.0572354 0.0581592
## 1 700001 1700000 71 0.0597176 0.0623646
## 1 800001 1800000 74 0.0611665 0.0633218MAN vs AUT
# Since we are using the SNP chip data, we can use a window of 1Mb and 100kb steps
vcftools --vcf output/fst/albo.vcf \
--weir-fst-pop output/fst/MAN.txt \
--weir-fst-pop output/fst/AUT.txt \
--fst-window-size 1000000 \
--fst-window-step 100000 \
--out output/fst/man_aut##
## VCFtools - 0.1.17
## (C) Adam Auton and Anthony Marcketta 2009
##
## Parameters as interpreted:
## --vcf output/fst/albo.vcf
## --fst-window-size 1000000
## --fst-window-step 100000
## --weir-fst-pop output/fst/MAN.txt
## --weir-fst-pop output/fst/AUT.txt
## --keep output/fst/MAN.txt
## --keep output/fst/AUT.txt
## --out output/fst/man_aut
##
## Warning: Expected at least 2 parts in INFO entry: ID=PR,Number=0,Type=Flag,Description="Provisional reference allele, may not be based on real reference genome">
## Keeping individuals in 'keep' list
## After filtering, kept 38 out of 60 Individuals
## Outputting Windowed Weir and Cockerham Fst estimates.
## Weir and Cockerham mean Fst estimate: 0.12225
## Weir and Cockerham weighted Fst estimate: 0.12987
## After filtering, kept 110353 out of a possible 110353 Sites
## Run Time = 1.00 secondsNEW vs AUT
# Since we are using the SNP chip data, we can use a window of 1Mb and 100kb steps
vcftools --vcf output/fst/albo.vcf \
--weir-fst-pop output/fst/NEW.txt \
--weir-fst-pop output/fst/AUT.txt \
--fst-window-size 1000000 \
--fst-window-step 100000 \
--out output/fst/new_aut##
## VCFtools - 0.1.17
## (C) Adam Auton and Anthony Marcketta 2009
##
## Parameters as interpreted:
## --vcf output/fst/albo.vcf
## --fst-window-size 1000000
## --fst-window-step 100000
## --weir-fst-pop output/fst/NEW.txt
## --weir-fst-pop output/fst/AUT.txt
## --keep output/fst/NEW.txt
## --keep output/fst/AUT.txt
## --out output/fst/new_aut
##
## Warning: Expected at least 2 parts in INFO entry: ID=PR,Number=0,Type=Flag,Description="Provisional reference allele, may not be based on real reference genome">
## Keeping individuals in 'keep' list
## After filtering, kept 50 out of 60 Individuals
## Outputting Windowed Weir and Cockerham Fst estimates.
## Weir and Cockerham mean Fst estimate: 0.10553
## Weir and Cockerham weighted Fst estimate: 0.12737
## After filtering, kept 110353 out of a possible 110353 Sites
## Run Time = 0.00 secondsNow we can import the results into R MAN vs NEW
# Define the file path using here
file_path <- here("output", "fst", "man_new.windowed.weir.fst")
# Import the file using read_delim
man_new <- read_delim(file_path,
delim = "\t", # Specify the tab delimiter
col_names = TRUE,
show_col_types = FALSE) # Indicates that the first row contains column names
# View the first few rows of the data
head(man_new)## # A tibble: 6 x 6
## CHROM BIN_START BIN_END N_VARIANTS WEIGHTED_FST MEAN_FST
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 1 1000000 22 -0.00302 -0.00497
## 2 1 100001 1100000 21 -0.00275 -0.00478
## 3 1 200001 1200000 28 0.00576 -0.000393
## 4 1 300001 1300000 36 0.0253 0.0213
## 5 1 400001 1400000 45 0.0530 0.0520
## 6 1 500001 1500000 42 0.0606 0.0571MAN vs AUT
# Define the file path using here
file_path <- here("output", "fst", "man_aut.windowed.weir.fst")
# Import the file using read_delim
man_aut <- read_delim(file_path,
delim = "\t", # Specify the tab delimiter
col_names = TRUE,
show_col_types = FALSE) # Indicates that the first row contains column names
# View the first few rows of the data
head(man_aut)## # A tibble: 6 x 6
## CHROM BIN_START BIN_END N_VARIANTS WEIGHTED_FST MEAN_FST
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 1 1000000 22 0.245 0.192
## 2 1 100001 1100000 21 0.245 0.188
## 3 1 200001 1200000 27 0.213 0.161
## 4 1 300001 1300000 36 0.192 0.150
## 5 1 400001 1400000 45 0.178 0.153
## 6 1 500001 1500000 42 0.134 0.125NEW vs AUT
# Define the file path using here
file_path <- here("output", "fst", "new_aut.windowed.weir.fst")
# Import the file using read_delim
new_aut <- read_delim(file_path,
delim = "\t", # Specify the tab delimiter
col_names = TRUE,
show_col_types = FALSE) # Indicates that the first row contains column names
# View the first few rows of the data
head(new_aut)## # A tibble: 6 x 6
## CHROM BIN_START BIN_END N_VARIANTS WEIGHTED_FST MEAN_FST
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 1 1000000 22 0.158 0.126
## 2 1 100001 1100000 21 0.155 0.120
## 3 1 200001 1200000 28 0.133 0.102
## 4 1 300001 1300000 37 0.127 0.101
## 5 1 400001 1400000 44 0.155 0.129
## 6 1 500001 1500000 41 0.133 0.113Combine the data
# Add a new column 'comparison' to each tibble
man_new$comparison <- "MAN_NEW"
new_aut$comparison <- "NEW_AUT"
man_aut$comparison <- "MAN_AUT"
# Combine the tibbles using rbind
combined_data <- rbind(man_new, new_aut, man_aut)
# View the first few rows of the combined data
head(combined_data)## # A tibble: 6 x 7
## CHROM BIN_START BIN_END N_VARIANTS WEIGHTED_FST MEAN_FST comparison
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 1 1 1000000 22 -0.00302 -0.00497 MAN_NEW
## 2 1 100001 1100000 21 -0.00275 -0.00478 MAN_NEW
## 3 1 200001 1200000 28 0.00576 -0.000393 MAN_NEW
## 4 1 300001 1300000 36 0.0253 0.0213 MAN_NEW
## 5 1 400001 1400000 45 0.0530 0.0520 MAN_NEW
## 6 1 500001 1500000 42 0.0606 0.0571 MAN_NEWSave the data
Now we plot
# Calculate the middle position
combined_data$mid_pos <- with(combined_data, (BIN_START + BIN_END) / 2)
# Reorder the levels of the comparison factor
combined_data$comparison <- fct_relevel(combined_data$comparison, "MAN_NEW", "MAN_AUT", "NEW_AUT")
# Define a set of pastel colors for three chromosomes
pastel_colors <- c("#ebd99f", "#B3CDE3", "#CCEBC5")
# Create the facet plot
window_fst <-
ggplot(combined_data,
aes(
x = mid_pos,
y = WEIGHTED_FST,
color = as.factor(CHROM)
)) +
geom_point(size = 0.1) +
facet_grid(comparison ~ CHROM, scales = "free_x") +
scale_x_continuous(labels = label_number(scale = 1e-6, suffix = "Mb")) +
scale_color_manual(values = pastel_colors) +
# Smooth line for MAN_NEW in black
geom_smooth(
data = subset(combined_data, comparison == "MAN_NEW"),
aes(group = CHROM),
se = FALSE,
span = 0.3,
color = "black"
) +
# Smooth line for MAN_AUT in red
geom_smooth(
data = subset(combined_data, comparison == "MAN_AUT"),
aes(group = CHROM),
se = FALSE,
span = 0.3,
color = "red"
) +
# Smooth line for NEW_AUT in red
geom_smooth(
data = subset(combined_data, comparison == "NEW_AUT"),
aes(group = CHROM),
se = FALSE,
span = 0.3,
color = "red"
) +
theme_bw() +
labs(x = "Position", y = "Weighted Fst", title = "") +
theme(
strip.text.x = element_text(angle = 0, hjust = .5),
strip.text.y = element_text(angle = 90, hjust = .5),
panel.spacing.x = unit(1.01, "lines"),
legend.position = "none"
)
window_fst## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
output_path <- here("output", "fst", "figures", "pairwise_estimates_windows.pdf")
ggsave(output_path, window_fst, height = 7, width = 9, dpi = 300)## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'Chromosome 2 peaks seems interesting. We can find out what are the windows that have such ligh fst values
# Lets find the windows for which the fst values are > 0.5
chr2 <- combined_data |>
dplyr::filter(
CHROM == "2"
) |>
dplyr::filter(
WEIGHTED_FST >= 0.5
)
# View the results
print(chr2)## # A tibble: 30 x 8
## CHROM BIN_START BIN_END N_VARIANTS WEIGHTED_FST MEAN_FST comparison mid_pos
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <fct> <dbl>
## 1 2 410300001 411300000 26 0.560 0.395 NEW_AUT 4.11e8
## 2 2 410400001 411400000 30 0.542 0.376 NEW_AUT 4.11e8
## 3 2 410500001 411500000 33 0.520 0.365 NEW_AUT 4.11e8
## 4 2 410600001 411600000 36 0.567 0.407 NEW_AUT 4.11e8
## 5 2 410700001 411700000 35 0.574 0.427 NEW_AUT 4.11e8
## 6 2 410800001 411800000 37 0.591 0.429 NEW_AUT 4.11e8
## 7 2 410900001 411900000 42 0.589 0.424 NEW_AUT 4.11e8
## 8 2 411000001 412000000 45 0.569 0.408 NEW_AUT 4.12e8
## 9 2 411100001 412100000 42 0.591 0.440 NEW_AUT 4.12e8
## 10 2 411200001 412200000 36 0.540 0.405 NEW_AUT 4.12e8
## # i 20 more rows4. Sliding windows Fst estimates between the sexes
Note: This section (sliding window Fst between sexes) was not used in the final version of the manuscript.
Clean env and memory
# Remove all objects from the environment
rm(list = ls())
# Run the garbage collector to free up memory
gc()## used (Mb) gc trigger (Mb) max used (Mb)
## Ncells 3366667 179.8 8277314 442.1 8277314 442.1
## Vcells 6449465 49.3 27019364 206.2 52772192 402.7We need to create a list of males and females for each population
# Specify the path using the 'here' function
fam_file_path <- here("output", "fst", "new_aut.fam")
# Now, read the file using the path - using data.table library
fam_data2 <- fread(fam_file_path, header = FALSE)
# Rename columns for clarity
setnames(fam_data2, c("FID", "IID", "PID", "MID", "Sex", "Phenotype"))
head(fam_data2)## FID IID PID MID Sex Phenotype
## <char> <int> <int> <int> <int> <int>
## 1: AUT_M 384 0 0 1 2
## 2: AUT_M 385 0 0 1 2
## 3: AUT_M 386 0 0 1 2
## 4: AUT_M 387 0 0 1 2
## 5: AUT_M 388 0 0 1 2
## 6: AUT_M 389 0 0 1 2We need the files to be in this format
## NEW_413
## NEW_414
## NEW_415
## NEW_416
## NEW_417
## NEW_418
## NEW_419
## NEW_420
## NEW_421
## NEW_422We can subset for each sex and create a file that we will use with vcftools
NEW males
new_males <- fam_data2 |>
dplyr::filter(FID == "NEW_M")
# Create ids to match what we have in the vcf
new_males$ID <- paste0(sub("_M", "_", new_males$FID), new_males$IID)
# Save it
write.table(new_males$ID, here("output", "fst","new_males_IDs.txt"), row.names = FALSE, col.names = FALSE, quote = FALSE)
head(new_males)## FID IID PID MID Sex Phenotype ID
## <char> <int> <int> <int> <int> <int> <char>
## 1: NEW_M 413 0 0 1 1 NEW_413
## 2: NEW_M 414 0 0 1 1 NEW_414
## 3: NEW_M 415 0 0 1 1 NEW_415
## 4: NEW_M 416 0 0 1 1 NEW_416
## 5: NEW_M 417 0 0 1 1 NEW_417
## 6: NEW_M 418 0 0 1 1 NEW_418Check it
## NEW_413
## NEW_414
## NEW_415
## NEW_416
## NEW_417
## NEW_418
## NEW_419
## NEW_420
## NEW_421
## NEW_422Now do it for the remaining files
NEW females
new_females <- fam_data2 |>
dplyr::filter(FID == "NEW_F")
# Create ids to match what we have in the vcf
new_females$ID <- paste0(sub("_F", "_", new_females$FID), new_females$IID)
# Save it
write.table(new_females$ID, here("output", "fst","new_females_IDs.txt"), row.names = FALSE, col.names = FALSE, quote = FALSE)
head(new_females)## FID IID PID MID Sex Phenotype ID
## <char> <int> <int> <int> <int> <int> <char>
## 1: NEW_F 424 0 0 2 1 NEW_424
## 2: NEW_F 425 0 0 2 1 NEW_425
## 3: NEW_F 426 0 0 2 1 NEW_426
## 4: NEW_F 427 0 0 2 1 NEW_427
## 5: NEW_F 428 0 0 2 1 NEW_428
## 6: NEW_F 429 0 0 2 1 NEW_429AUT males
aut_males <- fam_data2 |>
dplyr::filter(FID == "AUT_M")
# Create ids to match what we have in the vcf
aut_males$ID <- paste0(sub("_M", "_", aut_males$FID), aut_males$IID)
# Save it
write.table(aut_males$ID, here("output", "fst","aut_males_IDs.txt"), row.names = FALSE, col.names = FALSE, quote = FALSE)
head(aut_males)## FID IID PID MID Sex Phenotype ID
## <char> <int> <int> <int> <int> <int> <char>
## 1: AUT_M 384 0 0 1 2 AUT_384
## 2: AUT_M 385 0 0 1 2 AUT_385
## 3: AUT_M 386 0 0 1 2 AUT_386
## 4: AUT_M 387 0 0 1 2 AUT_387
## 5: AUT_M 388 0 0 1 2 AUT_388
## 6: AUT_M 389 0 0 1 2 AUT_389AUT females
aut_females <- fam_data2 |>
dplyr::filter(FID == "AUT_F")
# Create ids to match what we have in the vcf
aut_females$ID <- paste0(sub("_F", "_", aut_females$FID), aut_females$IID)
# Save it
write.table(aut_females$ID, here("output", "fst","aut_females_IDs.txt"), row.names = FALSE, col.names = FALSE, quote = FALSE)
head(aut_females)## FID IID PID MID Sex Phenotype ID
## <char> <int> <int> <int> <int> <int> <char>
## 1: AUT_F 400 0 0 2 2 AUT_400
## 2: AUT_F 401 0 0 2 2 AUT_401
## 3: AUT_F 402 0 0 2 2 AUT_402
## 4: AUT_F 403 0 0 2 2 AUT_403
## 5: AUT_F 404 0 0 2 2 AUT_404
## 6: AUT_F 405 0 0 2 2 AUT_405Now we can estimate fst with vcftools
NEW males vs NEW females
# Since we are using the SNP chip data, we can use a window of 1Mb and 100kb steps
vcftools --vcf output/fst/albo.vcf \
--weir-fst-pop output/fst/new_males_IDs.txt \
--weir-fst-pop output/fst/new_females_IDs.txt \
--fst-window-size 1000000 \
--fst-window-step 100000 \
--out output/fst/new_sex##
## VCFtools - 0.1.17
## (C) Adam Auton and Anthony Marcketta 2009
##
## Parameters as interpreted:
## --vcf output/fst/albo.vcf
## --fst-window-size 1000000
## --fst-window-step 100000
## --weir-fst-pop output/fst/new_males_IDs.txt
## --weir-fst-pop output/fst/new_females_IDs.txt
## --keep output/fst/new_males_IDs.txt
## --keep output/fst/new_females_IDs.txt
## --out output/fst/new_sex
##
## Warning: Expected at least 2 parts in INFO entry: ID=PR,Number=0,Type=Flag,Description="Provisional reference allele, may not be based on real reference genome">
## Keeping individuals in 'keep' list
## After filtering, kept 22 out of 60 Individuals
## Outputting Windowed Weir and Cockerham Fst estimates.
## Weir and Cockerham mean Fst estimate: 0.021818
## Weir and Cockerham weighted Fst estimate: 0.031432
## After filtering, kept 110353 out of a possible 110353 Sites
## Run Time = 0.00 secondsAUT males vs AUT females
# Since we are using the SNP chip data, we can use a window of 1Mb and 100kb steps
vcftools --vcf output/fst/albo.vcf \
--weir-fst-pop output/fst/aut_males_IDs.txt \
--weir-fst-pop output/fst/aut_females_IDs.txt \
--fst-window-size 1000000 \
--fst-window-step 100000 \
--out output/fst/aut_sex##
## VCFtools - 0.1.17
## (C) Adam Auton and Anthony Marcketta 2009
##
## Parameters as interpreted:
## --vcf output/fst/albo.vcf
## --fst-window-size 1000000
## --fst-window-step 100000
## --weir-fst-pop output/fst/aut_males_IDs.txt
## --weir-fst-pop output/fst/aut_females_IDs.txt
## --keep output/fst/aut_males_IDs.txt
## --keep output/fst/aut_females_IDs.txt
## --out output/fst/aut_sex
##
## Warning: Expected at least 2 parts in INFO entry: ID=PR,Number=0,Type=Flag,Description="Provisional reference allele, may not be based on real reference genome">
## Keeping individuals in 'keep' list
## After filtering, kept 28 out of 60 Individuals
## Outputting Windowed Weir and Cockerham Fst estimates.
## Weir and Cockerham mean Fst estimate: 0.025784
## Weir and Cockerham weighted Fst estimate: 0.039617
## After filtering, kept 110353 out of a possible 110353 Sites
## Run Time = 1.00 secondsNEW males vs AUT males
# Since we are using the SNP chip data, we can use a window of 1Mb and 100kb steps
vcftools --vcf output/fst/albo.vcf \
--weir-fst-pop output/fst/new_males_IDs.txt \
--weir-fst-pop output/fst/aut_males_IDs.txt \
--fst-window-size 1000000 \
--fst-window-step 100000 \
--out output/fst/new_aut_males##
## VCFtools - 0.1.17
## (C) Adam Auton and Anthony Marcketta 2009
##
## Parameters as interpreted:
## --vcf output/fst/albo.vcf
## --fst-window-size 1000000
## --fst-window-step 100000
## --weir-fst-pop output/fst/new_males_IDs.txt
## --weir-fst-pop output/fst/aut_males_IDs.txt
## --keep output/fst/new_males_IDs.txt
## --keep output/fst/aut_males_IDs.txt
## --out output/fst/new_aut_males
##
## Warning: Expected at least 2 parts in INFO entry: ID=PR,Number=0,Type=Flag,Description="Provisional reference allele, may not be based on real reference genome">
## Keeping individuals in 'keep' list
## After filtering, kept 26 out of 60 Individuals
## Outputting Windowed Weir and Cockerham Fst estimates.
## Weir and Cockerham mean Fst estimate: 0.12558
## Weir and Cockerham weighted Fst estimate: 0.16704
## After filtering, kept 110353 out of a possible 110353 Sites
## Run Time = 0.00 secondsNEW females vs AUT females
# Since we are using the SNP chip data, we can use a window of 1Mb and 100kb steps
vcftools --vcf output/fst/albo.vcf \
--weir-fst-pop output/fst/new_females_IDs.txt \
--weir-fst-pop output/fst/aut_females_IDs.txt \
--fst-window-size 1000000 \
--fst-window-step 100000 \
--out output/fst/new_aut_females##
## VCFtools - 0.1.17
## (C) Adam Auton and Anthony Marcketta 2009
##
## Parameters as interpreted:
## --vcf output/fst/albo.vcf
## --fst-window-size 1000000
## --fst-window-step 100000
## --weir-fst-pop output/fst/new_females_IDs.txt
## --weir-fst-pop output/fst/aut_females_IDs.txt
## --keep output/fst/new_females_IDs.txt
## --keep output/fst/aut_females_IDs.txt
## --out output/fst/new_aut_females
##
## Warning: Expected at least 2 parts in INFO entry: ID=PR,Number=0,Type=Flag,Description="Provisional reference allele, may not be based on real reference genome">
## Keeping individuals in 'keep' list
## After filtering, kept 24 out of 60 Individuals
## Outputting Windowed Weir and Cockerham Fst estimates.
## Weir and Cockerham mean Fst estimate: 0.097099
## Weir and Cockerham weighted Fst estimate: 0.12133
## After filtering, kept 110353 out of a possible 110353 Sites
## Run Time = 1.00 secondsNow we can import the results into R
NEW male vs female
# Define the file path using here
file_path <- here("output", "fst", "new_sex.windowed.weir.fst")
# Import the file using read_delim
new_sex <- read_delim(file_path,
delim = "\t", # Specify the tab delimiter
col_names = TRUE,
show_col_types = FALSE) # Indicates that the first row contains column names
# View the first few rows of the data
head(new_sex)## # A tibble: 6 x 6
## CHROM BIN_START BIN_END N_VARIANTS WEIGHTED_FST MEAN_FST
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 1 1000000 22 0.00698 0.00657
## 2 1 100001 1100000 21 0.00863 0.00819
## 3 1 200001 1200000 28 0.0149 0.0113
## 4 1 300001 1300000 35 0.0198 0.0159
## 5 1 400001 1400000 42 0.0248 0.0168
## 6 1 500001 1500000 39 0.0304 0.0197AUT male vs female
# Define the file path using here
file_path <- here("output", "fst", "aut_sex.windowed.weir.fst")
# Import the file using read_delim
aut_sex <- read_delim(file_path,
delim = "\t", # Specify the tab delimiter
col_names = TRUE,
show_col_types = FALSE) # Indicates that the first row contains column names
# View the first few rows of the data
head(aut_sex)## # A tibble: 6 x 6
## CHROM BIN_START BIN_END N_VARIANTS WEIGHTED_FST MEAN_FST
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 1 1000000 19 0.0268 0.0152
## 2 1 100001 1100000 19 0.0268 0.0152
## 3 1 200001 1200000 24 0.0296 0.0150
## 4 1 300001 1300000 33 0.0212 0.00665
## 5 1 400001 1400000 37 0.0278 0.0123
## 6 1 500001 1500000 35 0.0282 0.0131NEW males vs AUT males
# Define the file path using here
file_path <- here("output", "fst", "new_aut_males.windowed.weir.fst")
# Import the file using read_delim
new_aut_males <- read_delim(file_path,
delim = "\t", # Specify the tab delimiter
col_names = TRUE,
show_col_types = FALSE) # Indicates that the first row contains column names
# View the first few rows of the data
head(new_aut_males)## # A tibble: 6 x 6
## CHROM BIN_START BIN_END N_VARIANTS WEIGHTED_FST MEAN_FST
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 1 1000000 21 0.182 0.141
## 2 1 100001 1100000 20 0.179 0.135
## 3 1 200001 1200000 27 0.160 0.118
## 4 1 300001 1300000 36 0.159 0.116
## 5 1 400001 1400000 43 0.196 0.143
## 6 1 500001 1500000 40 0.171 0.125NEW females vs AUT females
# Define the file path using here
file_path <- here("output", "fst", "new_aut_females.windowed.weir.fst")
# Import the file using read_delim
new_aut_females <- read_delim(file_path,
delim = "\t", # Specify the tab delimiter
col_names = TRUE,
show_col_types = FALSE) # Indicates that the first row contains column names
# View the first few rows of the data
head(new_aut_females)## # A tibble: 6 x 6
## CHROM BIN_START BIN_END N_VARIANTS WEIGHTED_FST MEAN_FST
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 1 1000000 21 0.158 0.127
## 2 1 100001 1100000 20 0.157 0.123
## 3 1 200001 1200000 26 0.149 0.112
## 4 1 300001 1300000 35 0.133 0.102
## 5 1 400001 1400000 43 0.151 0.120
## 6 1 500001 1500000 40 0.134 0.106Combine the data
# Add a new column 'comparison' to each tibble
new_sex$comparison <- "NEW_sex"
aut_sex$comparison <- "AUT_sex"
new_aut_males$comparison <- "Males"
new_aut_females$comparison <- "Females"
# Combine the tibbles using rbind
combined_data <- rbind(new_sex, aut_sex, new_aut_males, new_aut_females)
# View the first few rows of the combined data
head(combined_data)## # A tibble: 6 x 7
## CHROM BIN_START BIN_END N_VARIANTS WEIGHTED_FST MEAN_FST comparison
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 1 1 1000000 22 0.00698 0.00657 NEW_sex
## 2 1 100001 1100000 21 0.00863 0.00819 NEW_sex
## 3 1 200001 1200000 28 0.0149 0.0113 NEW_sex
## 4 1 300001 1300000 35 0.0198 0.0159 NEW_sex
## 5 1 400001 1400000 42 0.0248 0.0168 NEW_sex
## 6 1 500001 1500000 39 0.0304 0.0197 NEW_sexNow we plot
# Calculate the middle position
combined_data$mid_pos <- with(combined_data, (BIN_START + BIN_END) / 2)
# Reorder the levels of the comparison factor
combined_data$comparison <- fct_relevel(combined_data$comparison, "NEW_sex", "AUT_sex", "Males", "Females")
# Define a set of pastel colors for three chromosomes
pastel_colors <- c("#ebd99f", "#B3CDE3", "#CCEBC5")
# Create the facet plot
window_fst <-
ggplot(combined_data,
aes(
x = mid_pos,
y = WEIGHTED_FST,
color = as.factor(CHROM)
)) +
geom_point(size = 0.1) +
facet_grid(comparison ~ CHROM, scales = "free_x") +
scale_x_continuous(labels = label_number(scale = 1e-6, suffix = "Mb")) +
scale_color_manual(values = pastel_colors) +
geom_smooth(
data = subset(combined_data, comparison == "NEW_sex"),
aes(group = CHROM),
se = FALSE,
span = 0.3,
color = "black"
) +
geom_smooth(
data = subset(combined_data, comparison == "AUT_sex"),
aes(group = CHROM),
se = FALSE,
span = 0.3,
color = "red"
) +
geom_smooth(
data = subset(combined_data, comparison == "Males"),
aes(group = CHROM),
se = FALSE,
span = 0.3,
color = "orange"
) +
geom_smooth(
data = subset(combined_data, comparison == "Females"),
aes(group = CHROM),
se = FALSE,
span = 0.3,
color = "orange"
) +
theme_bw() +
labs(x = "Position", y = "Weighted Fst", title = "") +
theme(
strip.text.x = element_text(angle = 0, hjust = .5),
strip.text.y = element_text(angle = 90, hjust = .5),
panel.spacing.x = unit(1.01, "lines"),
legend.position = "none"
)
window_fst## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
output_path <- here("output", "fst", "figures", "sex_pairwise_estimates_windows.pdf")
ggsave(output_path, window_fst, height = 7, width = 9, dpi = 300)## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'5. Check windows of interest
Chromosome 2 peaks seems interesting. We can find out what are the windows that have such ligh fst values
# Lets find the windows for which the fst values are > 0.5
chr2 <- combined_data |>
dplyr::filter(
CHROM == "2"
) |>
dplyr::filter(
WEIGHTED_FST >= 0.5
)
# View the results
print(chr2)## # A tibble: 23 x 8
## CHROM BIN_START BIN_END N_VARIANTS WEIGHTED_FST MEAN_FST comparison mid_pos
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <fct> <dbl>
## 1 2 127200001 128200000 3 0.555 0.316 Males 1.28e8
## 2 2 407500001 408500000 9 0.519 0.415 Males 4.08e8
## 3 2 410200001 411200000 24 0.564 0.357 Males 4.11e8
## 4 2 410300001 411300000 25 0.650 0.450 Males 4.11e8
## 5 2 410400001 411400000 29 0.636 0.426 Males 4.11e8
## 6 2 410500001 411500000 32 0.613 0.412 Males 4.11e8
## 7 2 410600001 411600000 34 0.660 0.471 Males 4.11e8
## 8 2 410700001 411700000 34 0.660 0.470 Males 4.11e8
## 9 2 410800001 411800000 36 0.673 0.470 Males 4.11e8
## 10 2 410900001 411900000 41 0.673 0.455 Males 4.11e8
## # i 13 more rowsFind the overlapping windows between males and females
# Find rows with the same BIN_START but different comparison values
overlapping_starts <- chr2 %>%
group_by(BIN_START) %>%
filter(n() > 1 & n_distinct(comparison) > 1) %>%
ungroup()
# View the results
print(overlapping_starts)## # A tibble: 6 x 8
## CHROM BIN_START BIN_END N_VARIANTS WEIGHTED_FST MEAN_FST comparison mid_pos
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <fct> <dbl>
## 1 2 410800001 411800000 36 0.673 0.470 Males 4.11e8
## 2 2 410900001 411900000 41 0.673 0.455 Males 4.11e8
## 3 2 411100001 412100000 41 0.681 0.476 Males 4.12e8
## 4 2 410800001 411800000 35 0.516 0.410 Females 4.11e8
## 5 2 410900001 411900000 40 0.515 0.411 Females 4.11e8
## 6 2 411100001 412100000 40 0.511 0.420 Females 4.12e8We have 3 windows. Starts on 410800001 and end on 412100000
## [1] 1299.999Now we can import the bin file that has the chromosomal scale and filter out by position
# Import the function
source(
here(
"notebooks", "helpers", "import_bim.R")
)
# Import the data
snps <- import_bim(here("output", "quality_control", "file7.bim"))
# Check it
head(snps)## # A tibble: 6 x 6
## Scaffold SNP Cm Position Allele1 Allele2
## <chr> <chr> <int> <dbl> <chr> <chr>
## 1 1 AX-581444870 0 97856 C T
## 2 1 AX-583035083 0 305518 A G
## 3 1 AX-583035102 0 308124 A G
## 4 1 AX-583033342 0 315059 C G
## 5 1 AX-583035163 0 315386 A G
## 6 1 AX-583033356 0 315674 C TNow we can select the window and get the SNP ids
# Filter the data for Scaffold 2 and Position range
window_snps <- snps %>%
filter(Scaffold == "2", Position >= 410800001, Position <= 412100000)
# View the results
print(window_snps) # 46 SNPs## # A tibble: 46 x 6
## Scaffold SNP Cm Position Allele1 Allele2
## <chr> <chr> <int> <dbl> <chr> <chr>
## 1 2 AX-579630382 0 411040898 C T
## 2 2 AX-579630409 0 411053076 G C
## 3 2 AX-579630426 0 411076185 G A
## 4 2 AX-579630432 0 411095416 T G
## 5 2 AX-579630434 0 411100836 G C
## 6 2 AX-579630438 0 411101050 T G
## 7 2 AX-579630465 0 411118768 T C
## 8 2 AX-579630489 0 411120151 T C
## 9 2 AX-579630502 0 411120383 A C
## 10 2 AX-579632017 0 411129367 T C
## # i 36 more rowsNow we can filter out it from our gene data
Import the data
## # A tibble: 6 x 8
## SNP Chromosome Position_chr Scaffold Position Gene_ID Start End
## <chr> <chr> <dbl> <chr> <dbl> <chr> <dbl> <dbl>
## 1 AX-581444870 1 97856 1.1 97856 LOC10939~ 64508 102530
## 2 AX-583033342 1 315059 1.1 315059 LOC10943~ 310594 316608
## 3 AX-583035163 1 315386 1.1 315386 LOC10943~ 310594 316608
## 4 AX-583033356 1 315674 1.1 315674 LOC10943~ 310594 316608
## 5 AX-583035257 1 442875 1.10 91677 LOC10940~ 66918 192829
## 6 AX-583035268 1 462944 1.10 111746 LOC10940~ 66918 192829The expression data
gene_expression <- read_delim(here("data", "files","MANvsAUTO_sig_mRNAs.csv"), delim = ",", col_names = TRUE, show_col_types = FALSE) |>
dplyr::select(
gene,log2FoldChange
) |>
dplyr::rename(
Gene_ID = gene
)
head(gene_expression)## # A tibble: 6 x 2
## Gene_ID log2FoldChange
## <chr> <dbl>
## 1 LOC115262812 -10.5
## 2 LOC109401291 -9.19
## 3 LOC109397830 8.73
## 4 LOC115264022 8.56
## 5 LOC115260314 8.51
## 6 LOC115258723 8.39First we can find how many genes we have for this window
fst_cluster_chr2 <- snps_genes_chr[snps_genes_chr$SNP %in% window_snps$SNP, ]
# How many genes in the scaffolds
length(unique(fst_cluster_chr2$Gene_ID))## [1] 1## [1] 20Check it
## # A tibble: 20 x 8
## SNP Chromosome Position_chr Scaffold Position Gene_ID Start End
## <chr> <chr> <dbl> <chr> <dbl> <chr> <dbl> <dbl>
## 1 AX-579630382 2 411040898 2.42 2888381 LOC1094~ 2.77e6 3.13e6
## 2 AX-579630409 2 411053076 2.42 2900559 LOC1094~ 2.77e6 3.13e6
## 3 AX-579630426 2 411076185 2.42 2923668 LOC1094~ 2.77e6 3.13e6
## 4 AX-579630432 2 411095416 2.42 2942899 LOC1094~ 2.77e6 3.13e6
## 5 AX-579630434 2 411100836 2.42 2948319 LOC1094~ 2.77e6 3.13e6
## 6 AX-579630438 2 411101050 2.42 2948533 LOC1094~ 2.77e6 3.13e6
## 7 AX-579630465 2 411118768 2.42 2966251 LOC1094~ 2.77e6 3.13e6
## 8 AX-579630489 2 411120151 2.42 2967634 LOC1094~ 2.77e6 3.13e6
## 9 AX-579630502 2 411120383 2.42 2967866 LOC1094~ 2.77e6 3.13e6
## 10 AX-579632017 2 411129367 2.42 2976850 LOC1094~ 2.77e6 3.13e6
## 11 AX-579632033 2 411129637 2.42 2977120 LOC1094~ 2.77e6 3.13e6
## 12 AX-579630663 2 411156954 2.42 3004437 LOC1094~ 2.77e6 3.13e6
## 13 AX-579630698 2 411157516 2.42 3004999 LOC1094~ 2.77e6 3.13e6
## 14 AX-579632196 2 411163842 2.42 3011325 LOC1094~ 2.77e6 3.13e6
## 15 AX-579632353 2 411204548 2.42 3052031 LOC1094~ 2.77e6 3.13e6
## 16 AX-579630892 2 411205467 2.42 3052950 LOC1094~ 2.77e6 3.13e6
## 17 AX-579632415 2 411217508 2.42 3064991 LOC1094~ 2.77e6 3.13e6
## 18 AX-579632471 2 411234189 2.42 3081672 LOC1094~ 2.77e6 3.13e6
## 19 AX-579632599 2 411284443 2.42 3131926 LOC1094~ 2.77e6 3.13e6
## 20 AX-579631153 2 411284926 2.42 3132409 LOC1094~ 2.77e6 3.13e6Is this gene also DE?
fst_cluster_chr2_de <- fst_cluster_chr2[fst_cluster_chr2$Gene_ID %in% gene_expression$Gene_ID, ]
# no it is not DE expressedNow we can check if any SNP from the LD cluster 14 is on this Fst peak
aut_ch2 <- readRDS(here("output", "ldna", "pop", "chr2", "AUT_clusters_snps.rds"))
str(aut_ch2$`3215_0.78`)## chr [1:4652] "AX-579450486" "AX-579456609" "AX-579459157" "AX-579459504" ...We can increase the window, for example check 100kb flanking windows
# Filter the data for Scaffold 2 and Position range
scaffold_242 <- snps_genes_chr %>%
filter(Scaffold == "2.42")
# View the results
length(unique(scaffold_242$Gene_ID)) # 92 genes on this scaffold## [1] 92Are any of these 92 genes DE?
fst_ld_de_genes <- scaffold_242[scaffold_242$Gene_ID %in% gene_expression$Gene_ID, ]
head(fst_ld_de_genes)## # A tibble: 1 x 8
## SNP Chromosome Position_chr Scaffold Position Gene_ID Start End
## <chr> <chr> <dbl> <chr> <dbl> <chr> <dbl> <dbl>
## 1 AX-579624053 2 409015071 2.42 862554 LOC11526~ 839135 863407