1.使用GCTA對群體進行PCA分析
相關結果數(shù)據(jù)見上篇簡書
#!/usr/bin/bash
###GCTA PCA分析
cd /public/home/*/data/test
#step1.準備二進制文件，如果是vcf文件车吹，可以直接使用plink軟件
#plink --vcf chr1D.merged.naive.SNP.rmhardfilter.vcf.gz --geno 0.1 --hwe 0.0001 --make-bed --out Out
plink --vcf chr1D.merged.naive.SNP.rmhardfilter.vcf.gz --pca --allow-extra-chr --out Out
#step2.生成用于PCA分析的matrix #--autosome 只分析常染色體 --pca后面的數(shù)字是需要的分組數(shù)
gcta64 --bfile Out --make-grm --autosome --out RESULT
gcta64 --grm RESULT --pca 5 --out out
#step3.使用R語言繪圖
tbl=read.table("out.eigenvec",header=TRUE)
head(tbl)
names(tbl) = c("Fid","ID","PC1","PC2","PC3")
plot(tbl$PC1,tbl$PC2,pch=c(rep(1,112),rep(2,103)),col=c(rep("blue",112),rep("red",103)),main="PCA",xlab="pc1",ylab="pc2")
legend("bottomright",c("TEXT1","TEXT2"),pch=c(rep(1),rep(2)),col=c(rep("blue"),rep("red")))
#plink --vcf *.vcf --pca -out all_genotypegvcf_plink
#vcftools --vcf *.vcf --plink --out all_genotypegvcf_vcftools
#plink --noweb --file all_genotypegvcf_vcftools --make-bed --out all_genotypegvcf_plink
#gcta64 --bfile all_genotypegvcf_plink --make-grm --autosome 22 --out all_genotypegvcf_plink
#gcta64 --grm all_genotypegvcf_plink --pca 20 --out all_genotypegvcf_plink
2.使用admixture進行群體的structure分析
#!/usr/bin/bash
cd /public/home/*/data/test
#將vcf文件轉化為ped文件
vcftools --gzvcf chr1D.merged.naive.SNP.rmhardfilter.vcf.gz --plink --out xj
#plink過濾ped文件生成bed文件
plink --noweb --file xj --geno 0.05 --maf 0.05 --hwe 0.0001 --make-bed --out QC
#admixture 計算
for K in 1 2 3 4 5 6 7;do admixture --cv QC.bed $K|tee log${K}.out;done
# K最小為最佳的分群
grep -h CV log*.out
# Rplot
tbl = read.table("QC.7.Q");
pdf("Q7.pdf")
barplot(t(as.matrix(tbl)),col= rainbow(7),xlab="Individual", ylab="Ancestry",border = NA,space = 0);
dev.off()
#plink --vcf *.vcf --out snp.clean --allow-extra-chr --make-bed
#plink --bfile snp.clean --out prunData --recode 12 --allow-extra-chr
#admixture --cv prunData.ped 2 >> log.txt
#for K in 2 3 4 5 6 7 8 9; do (admixture --cv prunData.ped $K | tee log${K}.out)& done
