2020.11.9【GWAS/WGS流程】丨全基因组关联分析绘图全流程_穆易青的博客-CSDN博客_gwas全基因组关联分析
一、数据筛选及格式转换
！！！注意：如果要同时使用，正确做法是先SNP（geno）后样本(mind)
1.去除多等位基因，indel
nohup /home/software/bcftools/bcftools view -m 2 -M 2 --type "snps" 75.vcf.gz -O z -o 75.recordsnps2vcf.gz &
2.格式转换（map ped to bam bed fam）
nohup/home/software/plink --allow-extra-chr --chr-set 29 --file Neogen_China_CHN_BOVG100V1_20201104 --make-bed --out Neogen_China_CHN_BOVG100V1_20201104 &
3.LD过滤 #（做进化树需要用到这一步）
# 基于连锁不平衡的SNP修剪（窗宽50、删除LD大于0.2的SNP对中一个、将窗口向前移动25个SNP）
nohup /home/software/plink --allow-extra-chr --chr-set 29 -bfile Neogen_China_CHN_BOVG100V1_20201104 --indep-pairwise 50 25 0.2 --out ld.Neogen_China_CHN_BOVG100V1_20201104-502502 &
# 提取过滤好的
nohup /home/software/plink --allow-extra-chr --chr-set 29 -bfile Neogen_China_CHN_BOVG100V1_20201104 --extract ld.Neogen_China_CHN_BOVG100V1_20201104-502502.prune.in --make-bed --out ld.Neogen_China_CHN_BOVG100V1_20201104-502502 &
4.最终QC进行质控，位点过滤（做群体结构需要这一步 ）剔除高缺失率（--geno ）和极低等位基因频率( --maf )的SNP
nohup /home/software/plink --allow-extra-chr --chr-set 29 -bfile ld.Neogen_China_CHN_BOVG100V1_20201104-502502 --geno 0.05 --maf 0.03 --make-bed --out ld.QC.Neogen_China_CHN_BOVG100V1_20201104-502502-geno005-maf03 &
# --geno 0.05 大于95%的个体都具有的变异位点才保留，其他去除；或者说保留检出率高于0.95的SNP。
# --maf 0.03 次等位基因频率，频率较低的第二等位基因的频率（防止假阳性）；将maf < 0.03的SNP筛选出来并过滤掉，仅包括MAF >= 0.03的SNP，默认值为0.03。
# --mind 0.10 去除基因型丢失率大于10%的个体样本；
# --hwe 0.0001 保留符合Hardy-Weinbery 的变异位点。
5.格式转换
# bed转成map ped(待定)
nohup /home/software/plink --allow-extra-chr --chr-set 29 --bfile ld.QC.Neogen_China_CHN_BOVG100V1_20201106-502502-geno02-maf03 --recode --out ld.QC.Neogen_China_CHN_BOVG100V1_20201106-502502-geno02-maf03 &
# 格式转换 将bed bim fam 转vcf文件
nohup /home/software/plink --allow-extra-chr --chr-set 29 --bfile test --recode vcf-iid --out test &
二、群体结构之构建进化树（22.10.25 注意数据不要ld过滤）
1.iqtree 构建进化树(基于每一个snp画，较准确但用时长）
# 转成fa
nohup python /home/hmy/ped2fa.py Neogen_China_CHN_BOVG100V1_20201106-geno02-maf03.ped ld.QC.Neogen_China_CHN_BOVG100V1_20201106-502502-geno02-maf03.fa &
# 过滤I成N
nohup sed -i "s/I/N/g" Neogen_China_CHN_BOVG100V1_20201104-geno02-maf03.fa &
# mafft比对
nohup mafft --auto Neogen_China_CHN_BOVG100V1_20201104-geno02-maf03.fa > Neogen_China_CHN_BOVG100V1_20201104-geno02-maf03_aligned.fa &
# iqtree 计算进化树
nohup /home/software/iqtree-1.6.12-Linux/bin/iqtree -s Neogen_China_CHN_BOVG100V1_20201104-geno02-maf03_aligned.fa -m TEST -st DNA -bb 1000 -nt AUTO &
# 9.iTOL画树( 有时画图需替换ID名称）
https://itol.embl.de/upload.cgi
输入XXX.bionj 文件
# ped2fa.py脚本
#!/usr/bin/env python
import sys
input_ped = sys.argv[1]
output_fa = sys.argv[2]
a1=open(input_ped , "r")
x = a1.readline()
b1 = open(output_fa , "w")
i=0
while x:
x = x.strip().split()
c= (">",x[0])
columns= ''.join(c)
b1.write(columns + "\n")
seq = ''.join(x[6:])
sequence = seq.replace("0", "N")
leng =int(len(sequence))
len1 =leng-1
all = ("A","T","C","G")
while i<len1:
b1.write(sequence[i])
i = i+2
b1.write("\n")
x = a1.readline()
i=0
a1.close
b1.close
2.孙老师脚本画（基于亲缘关系画the neighbor-joining phylogenetic tree was constructed by genetic distance (1-IBS)）
# 画进化树ibs
nohup /home/software/plink --allow-extra-chr --chr-set 29 --noweb --file com-328-100k-need_noneed.ped-geno02-maf03 --genome --out com-328-100k-need_noneed.ped-geno02-maf03 &
(使用质控后的ped和map文件，产生com-328-100k-need_noneed.ped-geno02-maf03.genome文件)
# 使用孙老师脚本（注意更改该脚本中的样品数！！！）
nohup perl 3____PLINK_genome_MEGA.pl &
# 使用mega软件查看进化树，mega软件中可删除样品看进化树！！！！"3-22-lzx_Dis.meg"在mega软件打开。
#3____PLINK_genome_MEGA.pl脚本：
#!usr/bin/perl
# define array of input and output files
open (AAA,"0-1-cattle-293-chr1-29-snp-indel.filter.pass-geno01-maf005.genome") || die "can't open AAA";
open (BBB,"/home/ysq/20210301-cattle-vcf-hebing/final-vcf/NEW-3-21/1-cattle-293-chr1-29-snp-indel.filter.pass-geno01-maf005.fam") || die "can't open BBB";
open (CCC,">3-22-lzx_Dis.meg");
my @aa=<AAA>;
my @bb=<BBB>;
$sample_size=293; ### 样品个数
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"; ##
}
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;
三、群体结构之admixture(注意：重测序数据需替换染色体NCto1)
！！！在构建admixture前需要考虑样本中每个群体中个体的数目，在总样本分析好后，与删除个别个体再次分析后，两次得到图形的颜色占比会发生变化。
# 使用bed文件,其中重测序数据中SNP没有名称，需要补上名称，使用awk:
awk '{print $1,$1":"$4,$3,$4,$5,$6}' XXX.bim > XXX-change.bim
sed -i 's/\t/ /g' XXX-change.bim # 去除vcf文件中的tab键
sed -i 's/ /\t/g' XXX-change.bim # 空格换成tab键
# 画图admixture可以使用gz压缩文件,做出bed
nohup /home/software/plink --allow-extra-chr --chr-set 27 --vcf XXX.vcf.gz --make-bed --out XXX
# ld过滤,maf过滤等
nohup /home/software/plink --allow-extra-chr --chr-set 29 -bfile test --indep-pairwise 50 25 0.2 --out ld.test-502502 &
nohup /home/software/plink --allow-extra-chr --chr-set 29 -bfile test --extract ld.test-502502.prune.in --make-bed --out ld.test-502502 &
#转成map ped
nohup /home/software/plink --allow-extra-chr --chr-set 27 --bfile test --recode --out test &
#格式转换 将bed bim fam 转vcf文件
nohup /home/software/plink --allow-extra-chr --chr-set 27 --bfile test --export vcf --out test &(.log/ .nosex/ .vcf文件)
#使用bed文件，计算K
##Error：Invalid chromosome code! Use integers.
# (判断因为染色体为NC格式，更改成1、2、3格式)替换染色体为数字后，转换为bed文件
# 计算K值
for K in 1 2 3 4 5; do /home/software/admixture_linux-1.3.0/admixture --cv combine_100k-328_hebing-bcftools.bed $K | tee log${K}.out; done
grep -h CV log*.out
#提取vcf文件中的个体list
nohup /home/software/bcftools/bcftools query -l test.vcf > test-id.txt &
nohup /home/software/vcftools/src/cpp/vcftools --vcf test.vcf --plink --out test & （vcf to map不会影响fid、iid，推荐）
nohup /home/software/plink --allow-extra-chr --chr-set 29 --file test --make-bed --out test &（map to bed不会影响fid、iid，推荐）
nohup /home/software/plink --allow-extra-chr --chr-set 29 --bfile test --recode vcf-iid --out test &（bed to vcf 不会影响fid、iid，推荐）
制作order.txt文件（也可为.csv文件）
三列： 1.地区Asia 2.ID名称 3.样本品种
1.fam文件中的IID（第二列）更换为自己想要的个体名称，可批量替换。2.制作一个ld.QC.Neogen_China_CHN_BOVG100V1_20201104-502502-geno02-maf03.order.txt文件(order.txt文件需要另存为制表符分隔格式文件)，可将order.txt文件的第三列改成真正的个体名称，这样图中就会显示每个个体名称。
3.可将第一列换成123...
4.表格中无空格
Africa Dp-8 Dp
Africa Dp-9 Dp
asia Han-1 STH
asia Han-10 STH
asia Han-2 STH
asia SRS1227526 Tan
asia SRS1227527 Tan
Middle_East ERS154526 wild
Middle_East ERS154528 wild
Middle_East ERS154529 wild
Middle_East ERS154530 wild
### R脚本画图
##### 1.未排序直接画图
tbl=read.table("hapmap3.3.Q")
barplot(t(as.matrix(tbl)), col=rainbow(3),xlab="Individual #", ylab="Ancestry", border=NA)
##### 2.排序后画图
#将Q文件，bed文件，fam文件，order.txt文件均放入同一个文件夹下
#!/usr/bin/env Rscript
## sort
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.com-328-100k-64-50k-need_noneed-502502-geno02-maf03"
max.k <- 10 ###文件夹下有多少Q文件，此处数值就写几
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)
#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, cex.names=0.6)
dev.off()
>画一个Q文件的脚本
tbl=read.table("/USER/zhusitao/Project/xj/reseq/result/11.admixture/QC.7.Q");
pdf("/USER/zhusitao/Project/xj/reseq/result/11.admixture/Q7.pdf")
barplot(t(as.matrix(tbl)),col= rainbow(7),xlab="Individual", ylab="Ancestry",border = NA,space = 0);
dev.off()
四、群体结构之主成分分析PCA
# 使用常染色体的vcf数据进行画图（染色体为NC格式，更改成1、2、3格式）
# 生成用于PCA分析的matrix
nohup /home/software/plink --allow-extra-chr --chr-set 27 --vcf ld.test3-502502-chr-33.vcf --make-bed --out ld.test3-502502-chr-33 &
nohup /home/software/gcta_1.92.3beta3/gcta64 --bfile ld.test3-502502-chr-33 --make-grm --autosome-num 27 --out ld.test3-502502-chr-33.gcta &
# --autosome是只分析常染色体
nohup /home/software/gcta_1.92.3beta3/gcta64 --grm ld.test3-502502-chr-33.gcta --pca 3 --out ld.test3-502502-chr-33.gcta.out &
# 得到的out.eigenvec即可用于下一步R作图使用。
制作ld.test3-502502-chr-33.txt表格,共四列,第1列为品种名,24列为样品名,3列为顺序数字
SAMM MLN-1 1 MLN-1
SAMM MLN-2 3 MLN-2
SAMM MLN-3 4 MLN-3
SAMM MLN-4 5 MLN-4
CM SRR11657659 11 SRR11657659
CM SRR11657660 12 SRR11657660
### R脚本画图
a=read.table("test.eigenvec",header=F)
head(a)
dim(a)
b=read.table("ld.test3-502502-chr-33.txt",header=F)
head(b)
library("ggplot2")
qplot(a[,3],a[,5],col=b[,1])
# 比较a文件3和5列,需要比较34,45
Breed=b[,1]
ggplot(data = a ,
aes(x = a[,3],
y = a[,4], #此处也做相应更改
group = Breed,
shape = Breed,
color = Breed)
)+geom_point(size=2) +scale_shape_manual(values = seq(0,75))
# 通常而言，PCA图会和系统发育树以及群体结构一起解释，相互验证。
###PCA.explain.R
m <- as.matrix(read.table("test.eigenval",header=F))
explainm=m/sum(m)
explainm[1:3] #此处上面步骤 --pca 3，一般3个可以解释
plot(explain)
PCA方差解释率显示
m <- as.matrix(read.table("test.eigenval",header=F))
explainm=m/sum(m)
explainm[1:3]
plot(explain)
五、LD衰减图
如果要理解LD衰减图，我们就必须先理解连锁不平衡（Linkagedisequilibrium，LD）的概念。连锁不平衡是由两个名词构成，连锁+不平衡。从一个类似的概念入手，更容易理解LD的概念，那就是基因的共表达。换句话来说，当位于某一座位的特定等位基因与另一座位的某一等位基因同时出现的概率大于群体中因随机分布的两个等位基因同时出现的概率时，就称这两个座位处于连锁不平衡状态（linkage disequilibrium）。如果两个SNP标记位置相邻，那么在群体中也会呈现基因型步调一致的情况。如果两个基因座是相关的，我们将会看到某些基因型往往共同遗传，即某些单倍型的频率会高于期望值。
这种不同基因座间的相关性，用一个数值来衡量就是D值。类似相关系数是标准化后的协方差，LD系数（r2）则是标准化后的D值，这个数值在0~1波动。r2=0就是两个位点完全不相关，群体中单倍型分布是随机的（观测值=期望值）。r2=1就是两个位点完全相关，某些基因型（A）只与特定的基因型（B）共同出现。
一般而言，两个位点在基因组上离得越近，相关性就越强，LD系数就越大。反之，LD系数越小。也就是说，随着位点间的距离不断增加，LD系数通常情况下会慢慢下降。这个规律，通常就会使用LD衰减图来呈现。LD衰减图就是利用曲线图来呈现基因组上分子标记间的平均LD系数随着标记间距离增加而降低的过程。
六、Haploview-连锁不平衡分析
当位于某一座位的特定等位基因与另一座位的某一等位基因同时出现的概率大于群体中因随机分布的两个等位基因同时出现的概率时，就称这两个座位处于连锁不平衡的状态。目前，连锁不平衡分析是群体进化的经典分析条目，分析的软件主要有plink和haploview。Haploview是一款分析单倍型的软件，依托于java形成的可视化界面，以下是整理的相关介绍及使用方法。
安装JAVA、配置环境，安装Haploview
提取位置文件
转换成map ped
将vcf文件转换为map ped文件、并制作info文件
awk '!/^#/{print $2"\t"$4 }' test.map > test.info
在软件中输入ped info文件后，出现错误：invalid affected status
haploview出现"invalid affected status"的解决方法 haploview弹出这种错误是因为haploview的缺失值默认为0，而plink文件的缺失值一般用”-9“表示，当ped文件的缺失值为”-9“时，haploview不认得该字符，就弹出错误提示。 解决方法为将文件中所有的”-9“替换为0。
解决：在linux界面中，用vim命令编辑ped文件后，输入以下字符
sed -i 's/-9/0/g' test.ped
”/-9/0/“表示将该文件的”-9“字符替换为”0“字符，g表示第一行到最后一行所有的”-9“字符都要进行替换。
在软件中输入ped与map文件
网友评论