GBS_tutorial/Script_5_Admixture.sh at main · rieseberglab/GBS_tutorial · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
#####################################
#
# Part 5: ADMIXTURE
# January 2025
#
# Code Contributor: Cassandra E. & Yue Y.
# Biodiversity Research Center, UBC
#
#####################################

# This code includes the following:
# 5.1: make BED files from filtered VCF files
# 5.2: run admixture for different k values with loop
# 5.3: bootstrap admixture results
# 5.4: Plot ADMIXTURE result in R


#####################################
#  5.1 VCF -> BED
#####################################
# make Plink file to be able to run Admixture

cd ~/scratch/GBS_workshop/5_ADMIXTURE/

# Copy the filtered VCF file from STEP 3 to the current directory
# cp /scratch/celphin/GBS_workshop/3_SNP_filtering/Cassiope_noMER_r10i.recode.vcf .

# convert the vcf file to Plink
vcftools --vcf Cassiope_noMER_r10i.recode.vcf --plink --out Cassiope_noMER_r10i

#MAKE A BED FILE
module load htslib/1.10.2
plink --file Cassiope_noMER_r10i --allow-no-sex --allow-extra-chr 0 --make-bed --out Cassiope_noMER_r10i

#Save sample names
module load bcftools
bcftools query -l Cassiope_noMER_r10i.recode.vcf > sample_39_names.txt


#####################################
#  5.2 run ADMIXTURE
#####################################
# see https://dalexander.github.io/admixture/admixture-manual.pdf for more details

tmux new-session -s Cassiope
tmux attach-session -t Cassiope
#salloc -c48 --time 03:00:00 --mem 30G --account def-rieseber

cd ~/scratch/GBS_workshop/5_ADMIXTURE

module load StdEnv/2020
module load admixture/1.3.0

mkdir Take1
mkdir bootstrap

cp Cassiope_noMER_r10i* Take1
cp Cassiope_noMER_r10i* bootstrap/

cd Take1
for K in 1 2 3 4 5 6 7 8 9; \
do admixture --cv=10 -s time -j48 -C 0.0000000001  Cassiope_noMER_r10i.bed $K | tee log${K}.out; done


#to get the CV errors and see which K value is the best model
grep -h CV log*

# Take 1
CV error (K=1): 0.54249
*CV error (K=2): 0.33333
CV error (K=3): 0.34385
CV error (K=4): 0.37351
CV error (K=5): 0.37667
CV error (K=6): 0.42367


#####################################
#  5.3 run ADMIXTURE (PERMUTATION 1000) - DO NOT RUN FOR THE WORKSHOP
#####################################

cd /scratch/celphin/GBS_workshop/4_ADMIXTURE/bootstrap

module load StdEnv/2020
module load admixture/1.3.0

#make one last run of K=2-6 with 1000 bootstraps to estimate Standard Errors of the estimates

for K in 2 3 4 5 6 ; \
do admixture -B1000 --cv=10 -s time -j48 -C 0.0000000001 Cassiope_noMER_r10i.bed $K | tee log${K}.out; done

# Yue Note: this is probably going to take a long time, are we running this in the workshop to check standarded errors


#####################################
#  5.4: Plot ADMIXTURE in R
#####################################

cd ~/scratch/GBS_workshop/5_ADMIXTURE

cp sample_39_names.txt Take1/

module load StdEnv/2020
module load r/4.2.1

R

# install.packages("tidyverse")
# install.packages("ggplot2")
# pkg install takes a while (~15min)

library(tidyverse)
library(ggplot2)
get

# set the workimg directory
getwd()
setwd("~/scratch/GBS_workshop/5_ADMIXTURE/Take1")

samplelist <- read_tsv("sample_39_names.txt",col_names = "sample")

all_data <- tibble(sample=character(),
                   k=numeric(),
                   Q=character(),
                   value=numeric())


for (k in 1:6){
  data <- read_delim(paste0("Cassiope_noMER_r10i.",k,".Q"),
                  col_names = paste0("Q",seq(1:k)),
                  delim=" ")

  data$sample <- samplelist$sample
  data$k <- k

  #This step converts from wide to long.
  data %>% gather(Q, value, -sample,-k) -> data
  all_data <- rbind(all_data,data)

  }

all_data


#Plot

plot <- all_data %>%
  ggplot(.,aes(x=sample,y=value,fill=factor(Q))) +
  geom_bar(stat="identity",position="stack") +
  xlab("Sample") + ylab("Ancestry") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 60, hjust = 1)) +
  scale_fill_brewer(palette="Set1",name="K",
                    labels=seq(1:5)) +
  facet_wrap(~k,ncol=1)


#Change width and length!! (according to need)
png(filename = "Cassiope.admixture.k1to6.png", width = 1400, height = 1100, res = 200)
print(plot)
dev.off()


# Using Globus, move to local computer, check final plot


# ---------- END (Yue 2025 Jan 28th)