群体分层与主成分分析教程 Population structure & PCA

原理简介
使用千人基因组项目数据与PLINK进行PCA
绘图
参考

prefix="ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/ALL.chr" ;

suffix=".phase3_shapeit2_mvncall_integrated_v5b.20130502.genotypes.vcf.gz" ;

for chr in {1..22}; do
    wget "${prefix}""${chr}""${suffix}" "${prefix}""${chr}""${suffix}".tbi ;
done

wget ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/working/20130606_sample_info/20130606_g1k.ped ;

wget ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/human_g1k_v37.fasta.gz ;
wget ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/human_g1k_v37.fasta.fai ;

gunzip human_g1k_v37.fasta.gz ;

for chr in {1..22}; do
    bcftools norm -m-any --check-ref w -f human_g1k_v37.fasta \
      ALL.chr"${chr}".phase3_shapeit2_mvncall_integrated_v5b.20130502.genotypes.vcf.gz | \
      bcftools annotate -x ID -I +'%CHROM:%POS:%REF:%ALT' | \
        bcftools norm -Ob --rm-dup both \
          > ALL.chr"${chr}".phase3_shapeit2_mvncall_integrated_v5b.20130502.genotypes.bcf ;

    bcftools index ALL.chr"${chr}".phase3_shapeit2_mvncall_integrated_v5b.20130502.genotypes.bcf ;
done

for chr in {1..22}; do
    plink --noweb \
      --bcf ALL.chr"${chr}".phase3_shapeit2_mvncall_integrated_v5b.20130502.genotypes.bcf \
      --keep-allele-order \
      --vcf-idspace-to _ \
      --const-fid \
      --allow-extra-chr 0 \
      --split-x b37 no-fail \
      --make-bed \
      --out ALL.chr"${chr}".phase3_shapeit2_mvncall_integrated_v5b.20130502.genotypes ;
done

mkdir Pruned ;

for chr in {1..22}; do
    plink --noweb \
      --bfile ALL.chr"${chr}".phase3_shapeit2_mvncall_integrated_v5b.20130502.genotypes \
      --maf 0.10 --indep 50 5 1.5 \
      --out Pruned/ALL.chr"${chr}".phase3_shapeit2_mvncall_integrated_v5b.20130502.genotypes ;

    plink --noweb \
      --bfile ALL.chr"${chr}".phase3_shapeit2_mvncall_integrated_v5b.20130502.genotypes \
      --extract Pruned/ALL.chr"${chr}".phase3_shapeit2_mvncall_integrated_v5b.20130502.genotypes.prune.in \
      --make-bed \
      --out Pruned/ALL.chr"${chr}".phase3_shapeit2_mvncall_integrated_v5b.20130502.genotypes ;
done

find . -name "*.bim" | grep -e "Pruned" > ForMerge.list ;

sed -i 's/.bim//g' ForMerge.list ;

plink --merge-list ForMerge.list --out Merge ;

plink --bfile Merge --pca

#Populations

This file describes the population codes where assigned to samples collected for the 1000 Genomes project. These codes are used to organise the files in the data_collections' project data directories and can also be found in column 11 of many sequence index files.

There are also two tsv files, which contain the population codes and descriptions for both the sub and super populations that were used in phase 3 of the 1000 Genomes Project:

<ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/20131219.populations.tsv>
<ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/20131219.superpopulations.tsv>

###Populations and codes

         CHB	Han Chinese             Han Chinese in Beijing, China
         JPT	Japanese                Japanese in Tokyo, Japan
         CHS	Southern Han Chinese    Han Chinese South
         CDX	Dai Chinese             Chinese Dai in Xishuangbanna, China
         KHV	Kinh Vietnamese         Kinh in Ho Chi Minh City, Vietnam
         CHD	Denver Chinese          Chinese in Denver, Colorado (pilot 3 only)
	
         CEU	CEPH                    Utah residents (CEPH) with Northern and Western European ancestry 
         TSI	Tuscan                  Toscani in Italia 
         GBR	British                 British in England and Scotland 
         FIN	Finnish                 Finnish in Finland 
         IBS	Spanish                 Iberian populations in Spain 
	
         YRI	Yoruba                  Yoruba in Ibadan, Nigeria
         LWK	Luhya                   Luhya in Webuye, Kenya
         GWD	Gambian                 Gambian in Western Division, The Gambia 
         MSL	Mende                   Mende in Sierra Leone
         ESN	Esan                    Esan in Nigeria
	
         ASW	African-American SW     African Ancestry in Southwest US  
         ACB	African-Caribbean       African Caribbean in Barbados
         MXL	Mexican-American        Mexican Ancestry in Los Angeles, California
         PUR	Puerto Rican            Puerto Rican in Puerto Rico
         CLM	Colombian               Colombian in Medellin, Colombia
         PEL	Peruvian                Peruvian in Lima, Peru

         GIH	Gujarati                Gujarati Indian in Houston, TX
         PJL	Punjabi                 Punjabi in Lahore, Pakistan
         BEB	Bengali                 Bengali in Bangladesh
         STU	Sri Lankan              Sri Lankan Tamil in the UK
         ITU	Indian                  Indian Telugu in the UK

Should you have any queries, please contact info@1000genomes.org.

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

#读取文件
pca = pd.read_csv("1000_all.eigenvec"," ",header=None)
ped = pd.read_csv("20130606_g1k.ped","\t")
#重命名pca文件的列
pca = pca.rename(columns=dict([(1,"Individual ID")]+[(x,"PC"+str(x-1)) for x in range(2,22)]))

#join两个数据表
pcaped=pca.join(ped,on="Individual ID",how="inner")

#给个体加上super population的标签
pcaped["Superpopulation"]="unknown"
pcaped.loc[pcaped["Population"].isin(["JPT","CHB","CDX","CHS","JPT","KHV","CHD"]),"Superpopulation"]="EAS"
pcaped.loc[pcaped["Population"].isin(["BEB","GIH","ITU","PJL","STU"]),"Superpopulation"]="SAS"
pcaped.loc[pcaped["Population"].isin(["CLM","MXL","PEL","PUR"]),"Superpopulation"]="AMR"
pcaped.loc[pcaped["Population"].isin(["CEU","FIN","GBR","IBS","TSI"]),"Superpopulation"]="EUR"
pcaped.loc[pcaped["Population"].isin(["ACB","ASW","ESN","GWD","LWK","MSL","YRI"]),"Superpopulation"]="AFR"

#绘制PC1与PC2的散点图，以Super population分不同颜色
plt.figure(figsize=(15,15))
sns.scatterplot(data=pcaped,x="PC1",y="PC2",hue="Superpopulation",s=50)

#绘制PC1与PC2的散点图，以Population分不同颜色
sns.scatterplot(data=pcaped,x="PC1",y="PC2",hue="Population",s=50)