常見(jiàn)的群體結(jié)構(gòu)的分析方法有admixture分析隔节、系統(tǒng)發(fā)生數(shù)分析以及主成分分析等鹅经。
1、admixture分析
###過(guò)濾數(shù)據(jù)
常用plink軟件過(guò)濾怎诫，在此就不做介紹了瘾晃，直接開(kāi)始后續(xù)操作。
###dmixture進(jìn)行群體遺傳結(jié)構(gòu)分析（群體數(shù)自己決定）
for K in 3 4 5 6 7; do /home/software/admixture_linux-1.3.0/admixture --cv ld.QC.75_noinclude0-502502-geno02-maf03.bed $K | tee log${K}.out; done
###提取CV值：CV error最小的為最佳K值
grep -h CV log*.out
分析結(jié)束后生成了自己設(shè)定k值的Q文件幻妓，用于在R中繪圖
1）R語(yǔ)言繪圖
admixture的可視化分為兩種
###最佳K值的可視化
ta1 = read.table("ld.QC.75_noinclude0-502502-geno02-maf03.ped.map.3.4.Q") ##用的是最佳K值的那個(gè)Q文件
head(ta1)
barplot(t(as.matrix(ta1)),col = rainbow(3),
xlab = "Individual",
ylab = "Ancestry",
border = NA)
####全部K值的可視化(較復(fù)雜）
利用表格根據(jù)fam文件（三列 1.地區(qū)Asia 2.ID名稱蹦误，與fam文件的一致 3.樣本品種）作三列的order.txt并用制表符分隔形式保存
（###可將order.txt文件的第一列地區(qū)Asia改成真正的個(gè)體名稱，這樣圖中就會(huì)顯示每個(gè)個(gè)體名稱
###可將order.txt文件中的順序進(jìn)行調(diào)整則圖中的順序即為order.txt文件的個(gè)體順序）
Session中上傳工作目錄肉津，需建立一個(gè)文件夾（包括Q文件强胰，fam、bed妹沙、bim文件偶洋，order.txt文件）
##安裝軟件
install packages（Rcolorbrewer）
##（導(dǎo)入含有Q order.txt bed bid fam的文件夾,修改以下程序中的文件名和K值）
sort.admixture <- function(admix.data){
## sort columns according to the cor
k <- length(admix.data)
n.ind <- nrow(admix.data[[1]])
name.ind <- rownames(admix.data[[1]])
admix.sorted <- list()
if (admix.data[[1]][1,1] > admix.data[[1]][1,2]){
admix.sorted[[1]] <- admix.data[[1]]
}else{
admix.sorted[[1]] <- admix.data[[1]][,c(2,1)]
}
for (i in 1:(k-1)){
admix <- matrix(nrow = n.ind, ncol = (i + 2))
cors <- cor(admix.sorted[[i]], admix.data[[i + 1]])
sorted.loc <- c()
for (j in 1:nrow(cors)){
cor <- cors[j,]
cor[sorted.loc] <- NA
sorted.loc <- c(sorted.loc, which.max(cor))
}
sorted.loc <- c(sorted.loc, which(! 1:ncol(cors) %in% sorted.loc))
cat("n_max = ", sorted.loc, "\n")
admix <- admix.data[[i + 1]][,sorted.loc]
rownames(admix) <- name.ind
admix.sorted[[i + 1]] <- admix
}
return(admix.sorted)
}
sort.iid <- function(k.values, groups){
##k.values <- admix.values[[1]]
##groups <- admix.fam
max.col <- which.max(colSums(k.values))
k.values <- cbind(k.values, groups[match(rownames(k.values), as.character(groups$iid)),])
k.values <- transform(k.values, group = as.factor(k.values$fid))
k.means <- tapply(k.values[,max.col], k.values$group, mean)
k.means <- k.means[order(k.means)]
k.sort <- data.frame(id = names(k.means),
order = order(k.means),
mean = k.means)
k.values$order <- k.sort[match(as.character(k.values$group), k.sort$id), 3]
k.values <- k.values[order(k.values$order, k.values[,max.col]),]
return(rownames(k.values))
}
sort.fid <- function(iid.order, fid.order, fam.table){
new.order <- c()
for (fid in fid.order){
new.order <- c(new.order, which(iid.order %in% fam.table[fam.table$fid == fid, "iid"]))
}
return(iid.order[new.order])
}
read.structure <- function(file, type = "structure"){
if (type == "structure"){
k.values <- read.table(file = file, header = F)
rownames(k.values) <- k.values[,1]
k.values[,1:3] <- NULL
}else{
k.values <- read.table(file = file, header = F)
}
return(k.values)
}
add.black.line <- function(data, groups, nline = 1){
# data <- as.matrix(plot.data[[2]])
# groups <- group.name
# nline <- 3
data <- as.matrix(data)
group.name <- unique(groups)
new.data <- matrix(NA, ncol = ncol(data))
black.data <- matrix(NA, nrow = nline, ncol = ncol(data))
new.name <- c(NA)
for (name in group.name){
new.data <- rbind(new.data, black.data)
new.data <- rbind(new.data, data[which(groups == name),])
new.name <- c(new.name, rep(NA,nline))
new.name <- c(new.name, rownames(data)[which(groups == name)])
}
added.data <- new.data[(nline + 2):nrow(new.data),]
rownames(added.data) <- new.name[(nline + 2):nrow(new.data)]
return(added.data)
}
##============================================================================
##
header <- "ld.QC.75_noinclude0-502502-geno02-maf03" ######plink格式的文件波闹，改為自己使用的文件名）
max.k <- 4
admix.fn <- paste(header, 2:max.k, "Q", sep = ".")
fam.fn <- paste(header, "fam", sep = ".")
admix.fam <- read.table(fam.fn, stringsAsFactors = F,
col.names = c("fid", "iid", "pid", "mid", "sex", "pheno"))
admix.values <- lapply(admix.fn, read.table, header = F,
row.names = as.character(admix.fam$iid))
order.fn <- paste(header, "order.txt", sep = ".")
admix.order <- read.table(order.fn, col.names = c("region", "iid", "fid"), stringsAsFactors = F)
id.order <- admix.order$iid
#id.order <- sort.iid(admix.values[[1]], admix.order)
admix.data <- list()
for (i in 1:length(admix.values)){
admix.data[[i]] <- admix.values[[i]][id.order,]
}
species <- as.character(admix.order[,1])
sorted.data <- sort.admixture(admix.data)
## add black line in plot
nline <- 1
plot.data <- list()
group.order <- admix.order[match(id.order, admix.order$iid),3]
for (i in 1:length(sorted.data)){
plot.data[[i]] <- add.black.line(sorted.data[[i]], group.order, nline = nline)
}
## add xlab to plot
plot.id.list <- rownames(plot.data[[1]])
#plot.xlab <- admix.fam[match(x = plot.id.list, table = admix.fam$iid),]
plot.xlab <- admix.order[match(x = plot.id.list, table = admix.order$iid),]
plot.lab <- unique(plot.xlab$fid)
plot.lab <- plot.lab[!is.na(plot.lab)]
plot.at <- c()
start <- 0
for (fid in plot.lab){
xlen <- length(which(plot.xlab$fid == fid))
gap <- start + floor(xlen / 2)
plot.at <- c(plot.at, gap)
start <- start + nline + xlen
}
##=============================================================================
## barplot admixture and structure
library(RColorBrewer)
my.colours <- c(brewer.pal(8, "Dark2"), "mediumblue", "darkred", "coral4",
"purple3", "lawngreen", "dodgerblue1", "paleturquoise3",
"navyblue", "green3", "red1", "cyan",
"orange", "blue", "magenta4", "yellowgreen", "darkorange3",
"grey60", "black")
max.k <- length(plot.data)
n <- dim(plot.data[[1]])[1]
#pdf(file=paste(header, "admix.plot.pdf", sep = "."), width = 16, height = 12)
png(file=paste(header, "admix.plot.png", sep = "."), res=300, width = 2400, height = 2000)
par(mfrow = c(max.k,1), mar=c(0.5,2,0,0), oma=c(6,0,1,0))
par(las=2)
#for (i in 1:(max.k - 1)){
for (i in 1:max.k){
barplot(t(as.matrix(plot.data[[i]])), names.arg = rep(c(""), n),
col = my.colours, border = NA, space = 0, axes = F,
ylab = "")
axis(side = 2, at = 0.5, labels = as.character(i+1), tick = F, hadj = 0)
}
axis(side = 1, at = plot.at, labels = plot.lab, tick = F, lty = 15, cex.axis = 0.5)
#barplot(t(as.matrix(plot.data[[max.k]])), names.arg = rownames(plot.data[[1]]), axes = F,
# col = my.colours, border = NA, space = 0,
# ylab = paste("K=", max.k + 1, sep = ""), cex.axis=0.6(坐標(biāo)軸字體大小), cex.names=0.6)
dev.off()
2）pong對(duì)ADMIXTURE結(jié)果進(jìn)行可視化（真神器）
網(wǎng)址在這兒：https://github.com/ramachandran-lab/pong
1摇幻、簡(jiǎn)便方法，在服務(wù)器做：
pong -m pong_filemap.txt -i nd2pop.txt -n pop_order.txt -l color.txt
運(yùn)行后顯示：
pong server is now running locally & listening on port 4000
Open your web browser and navigate to http://localhost:4000 to see the visualization
在本地瀏覽器打開(kāi)http://localhost:4000（將localhost改成你的服務(wù)器IP）
2君珠、如果想在自己電腦做
推薦在自己電腦下載安裝python肾筐，并勾選配置環(huán)境哆料，會(huì)自動(dòng)配置所需環(huán)境（不然后續(xù)很麻煩）
# win+R打開(kāi)系統(tǒng)界面，輸入cmd并回車
# 檢查是否已安裝pip
python -m pip --version
# 一般是已經(jīng)安裝好了的吗铐，使用pip將pong安裝在自己電腦
pip install pong
# 安裝完成之后就可以用了东亦，這里展示windows的用法（附帶路徑的）
pong -m D:\桌面\pong\pong_filemap.txt -i D:\桌面\pong\ind2pop.txt -n D:\桌面\pong\pop_order.txt -l D:\桌面\pong\color.txt
-m filemap文件，里面是三列文件唬渗，第一列r1u1（字母是隨便編寫典阵，第二個(gè)數(shù)字和K對(duì)應(yīng)，如果想對(duì)K重復(fù)跑幾次镊逝，改變第一個(gè)數(shù)字） 第二列為K壮啊，一行一個(gè)值 第三列為Q文件
-i ind2pop文件，一列撑蒜，為個(gè)體對(duì)應(yīng)的群體ID歹啼，一個(gè)個(gè)體一個(gè)
-n order文件，群體ID進(jìn)行排序座菠，一個(gè)群體一個(gè)ID就好
-l color文件狸眼，這個(gè)參數(shù)可不加，如果嫌棄默認(rèn)顏色丑浴滴，可加拓萌，一行一個(gè)顏色，可為RGB升略，16進(jìn)制以及英文
運(yùn)行完后顯示：
Open your web browser and navigate to http://localhost:4000 to see the visualization
復(fù)制http://localhost:4000到瀏覽器打開(kāi)即可
2微王、系統(tǒng)發(fā)生樹分析
有很多種建樹的方法屡限，我們這里介紹下MEGA用NJ法建樹和IQ-TREE用ML法建樹
過(guò)濾的步驟略去
1）MEGA建樹
###計(jì)算genome
/home/software/plink --bfile XXX-502502-geno02-maf03 --allow-extra-chr --chr-set 26 --genome
編寫perl腳本 .pl結(jié)尾（Linux下新建.pl文件，直接進(jìn)入編輯即可炕倘，當(dāng)然囚霸，底下這個(gè)是現(xiàn)成的，改改可以直接用）
（第一個(gè)open里改成自己的genome文件
第二個(gè)open里改成自己的fam文件
第三個(gè)open里將>后面的改成輸出的文件名.meg
在下面的sample_size里激才，將數(shù)量改成自己使用的數(shù)量）
#!usr/bin/perl
# define array of input and output files
open (AAA,"ld.QC.0-1-cattle-290-chr1-29-snp-indel.filter.pass-502502-geno01-maf00.genome.genome") || die "can't open AAA"; ##用自己上一步的genome文件
open (BBB,"290-cattle-breed.fam.txt") || die "can't open BBB"; ##使用自己的fam文件
open (CCC,">4-22-lzx_Dis.meg"); ##輸出文件名在>后面改
my @aa=<AAA>;
my @bb=<BBB>;
$sample_size=290; ### 個(gè)體數(shù)目，改成自己用的數(shù)目
print CCC "#mega\n!Title: $sample_size pigs;\n!Format DataType=Distance DataFormat=UpperRight NTaxa=$sample_size;\n\n";
foreach ($num1=0;$num1<=$#bb;$num1++){
chomp $bb[$num1];
@arraynum1=split(/\s+/,$bb[$num1]);
print CCC "#$arraynum1[1]\n"; ##個(gè)體的ID名稱
}
print CCC "\n";
@array=();
foreach ($num2=1;$num2<=$#aa;$num2++){
chomp $aa[$num2];
@arraynum1=split(/\s+/,$aa[$num2]);
push(@array,1-$arraynum1[12]);
}
@array2=(0);
$i=$sample_size;
while ($i>0){
push(@array2,$array2[$#array2]+$i);
$i=$i-1;
}
print "@array2";
for ($i=($sample_size-1); $i>=0; $i=$i-1){
print CCC " " x ($sample_size-($i+1));
for ($j=$array2[$sample_size-$i-1]; $j<=$array2[$sample_size-$i]-1; $j++){
print CCC "$array[$j] ";
}
print CCC "\n";
}
close AAA;
close BBB;
close CCC;
#######生成.meg文件利用MEGA軟件可視化
2)IQ-TREE建樹
#####轉(zhuǎn)成map ped
/home/software/plink --allow-extra-chr --chr-set 27 -bfile ld.QC.xll_noinclude0.recode-502502-geno02-maf03 --recode --out ld.QC.xll_noinclude0.recode-502502-geno02-maf03
#####轉(zhuǎn)成fa
nohup python /home/hmy/ped2fa.py ld.QC.xll_noinclude0.recode-502502-geno02-maf03.ped ld.QC.xll_noinclude0.recode-502502-geno02-maf03.fa &
打開(kāi)fa文件修改個(gè)體ID额嘿，當(dāng)然也可以不改
#####過(guò)濾I成N
sed -i "s/I/N/g" ld.QC.xll_noinclude0.recode-502502-geno02-maf03.fa
#####mafft比對(duì)
####多序列比對(duì)：是指把多條（3 條或以上）有系統(tǒng)進(jìn)化關(guān)系的蛋白質(zhì)或核酸序列進(jìn)行比對(duì)瘸恼，盡可能地把相同的堿基或氨基酸殘基排在同一列上。
###這樣做的意義是册养，對(duì)齊的堿基或氨基酸殘基在進(jìn)化上是同源的东帅，即來(lái)自共同祖先（common ancestor）。
nohup mafft --auto ld.QC.xll_noinclude0.recode-502502-geno02-maf03.fa > ld.xll_noinclude0.recode-502502-geno02-maf03_aligned.fa &
#####Iqtree（構(gòu)建進(jìn)化樹）
nohup /home/software/iqtree-1.6.12-Linux/bin/iqtree -s ld.QC.75_xiugai_noinclude0-502502-geno02-maf03.fa -m TEST -st DNA -bb 1000 -nt AUTO &
#####itol畫樹（網(wǎng)頁(yè)搜索）
用bionj文件
3）RAxML建樹：
vcf2phylip.py鏈接：https://github.com/edgardomortiz/vcf2phylip/blob/master/vcf2phylip.py
####1 轉(zhuǎn)為phy格式：
python /home/sll/software/vcf2phylip.py --input FASN.vcf.recode.vcf
####2 建樹：
raxmlHPC-PTHREADS-SSE3 -f a -m GTRGAMMA -p 12345 -x 12345 -# 10 -s FASN.min4.phy -n raxml -T 30
-f a此參數(shù)用于選擇 RAxML 運(yùn)算的算法球拦】勘眨可以設(shè)定的值非常之多。 a 表示執(zhí)行快速 Bootstrap 分析并搜索最佳得分的 ML 樹坎炼。
-x 12345指定一個(gè) int 數(shù)作為隨機(jī)種子愧膀，以啟用快速 Bootstrap 算法。
-p 12345指定一個(gè)隨機(jī)數(shù)作為 parsimony inferences 的種子谣光。
-# 100指定 bootstrap 的次數(shù)檩淋。
-m PROTGAMMALGX 指定核苷酸或氨基酸替代模型。PROTGAMMALGX 的解釋： "PROT" 表示氨基酸替代模型萄金； GAMMA 表示使用 GAMMA 模型蟀悦； X 表示使用最大似然法估計(jì)堿基頻率。
-s ex.phy指定輸入文件氧敢。phy 格式的多序列比對(duì)結(jié)果日戈。軟件包中包含一個(gè)程序來(lái)將 fasta 格式轉(zhuǎn)換為 phy 格式。
-n ex輸出文件的后綴為 .ex 孙乖。
-T 20指定多線程運(yùn)行的 CPUs 浙炼。
3、主成分分析
（gcta）（gcta用的是數(shù)字染色體唯袄，若不是鼓拧，注意替換下）
###make germ
nohup /home/software/gcta_1.92.3beta3/gcta64 --bfile ld.QC.75_noinclude0-502502-geno02-maf03 --make-grm --autosome-num 26 --out ld.QC.75_noinclude0-502502-geno02-maf03.gcta &
# autosome表示只選出常染色體來(lái)運(yùn)行
###pca
nohup /home/software/gcta_1.92.3beta3/gcta64 --grm ld.QC.75_noinclude0-502502-geno02-maf03.gcta --pca 5 --out ld.QC.75_noinclude0-502502-geno02-maf03.gcta.out &
########輸入的文件就是.gcta,這個(gè)文件不是上一步的生成文件，pca后面跟的是要做幾個(gè)主成分的比較
R繪圖
################
a=read.table("pca.eigenvec",header=F)
head(a)
dim(a)
b=read.table("pca.txt",header=F) ######pca.txt為四列的文件越妈，（1列為品種季俩，2列和4列為個(gè)體ID，3列為排序數(shù)字）
head(b)
library("ggplot2")
qplot(a[,3],a[,4],col=b[,1])
Breed=b[,1]
p = ggplot(data = a ,
aes(x = a[,3],
y = a[,4], ############x = a[,3], y = a[,4]代表第3列和第4列比較也就是pc1與pc2比較
group = Breed,
shape = Breed,
color = Breed)
)+geom_point(size=2) +scale_shape_manual(values = seq(0,75))
p + labs(x = "pc1", y = "pc2") #############修改x軸和y軸的坐標(biāo)名稱
