mawk '!/#/' xxx.vcf | wc -l
1.新建了conda 2brad 环境
2.安装一系列软件，pear samtools bcftools plink，其中由于测序数据是双端测序150bp产生，2brad中2b型限制性内切酶BsaxI产生等长33bp序列，因此要对原始双端测序的数据进行merge和拆分。其中merge拆分时使用2bRADExtraction.pl 脚本，对应样本名称后可直接产生每个样本的酶切数据。
3.使用bowtie2进行比对，产生sam文件转成sortedbam文件
4.使用bcftools进行snp calling，主要命令
bcftools mpileup -Ou -f ../../02.ref/Gipl.BsaXI.fa *bam | bcftools call -mv -Ob -o calls.bcf
5.使用vcftools 进行过滤
vcftools --vcf calls.vcf --max-missing 0.5 --mac 3 --minQ 30 --recode --recode-INFO-all --out xxx
6.使用plink进行主成分分析
plink --vcf xx.vcf --recode --out xx --const-fid 0 --allow-extra-chr
上述命令产生了.map,.nosex和.ped结尾的三个文件，利用.ped 文件产生.bed 文件
plink --allow-extra-chr --file xx --noweb --make-bed --out xx
上述命令产生bed 和 bim，下一步计算PCA
plink --allow-extra-chr --threads 20 -bfile xx --pca 20 --out xx
产生的两个文件可以利用ggplot2画图
library("ggplot2")
library("tidyverse")
df_pcs<- readxl::read_xlsx("new1.xlsx")
re1a<-read.table ("new1.eigenval",header = F)
#percentage<-eigval$V1/sum(eigval$V1)*100
re1a$por<-re1a$V1
ggplot(df_pcs, aes(x=PC1,y=PC2)) +
geom_point(aes(color= Group, shape= Group),size=1.5)+
labs(x="PC1",y="PC2")+theme_bw()+
theme(panel.border=element_blank(),
panel.grid.major=element_blank(),
panel.grid.minor=element_blank(),
axis.line= element_line(colour = "black"))+
xlab(paste0("PC1(",round(re1a$por[1],2),"%)"))+ ylab(paste0("PC2(",round(re1a$por[2],2),"%)"))
ggplot(df_pcs, aes(x=PC1,y=PC3)) +
geom_point(aes(color= Group, shape= Group),size=1.5)+
labs(x="PC1",y="PC3")+theme_bw()+
theme(panel.border=element_blank(),
panel.grid.major=element_blank(),
panel.grid.minor=element_blank(),
axis.line= element_line(colour = "black"))+
xlab(paste0("PC1(",round(re1a$por[1],2),"%)"))+ ylab(paste0("PC3(",round(re1a$por[3],2),"%)"))
ggplot(df_pcs, aes(x=PC2,y=PC3)) +
geom_point(aes(color= Group, shape= Group),size=1.5)+
labs(x="PC2",y="PC3")+theme_bw()+
theme(panel.border=element_blank(),
panel.grid.major=element_blank(),
panel.grid.minor=element_blank(),
axis.line= element_line(colour = "black"))+
xlab(paste0("PC2(",round(re1a$por[2],2),"%)"))+ ylab(paste0("PC3(",round(re1a$por[3],2),"%)"))
