利用GCAT工具做PCA分析

阿凡亮

发布于 2020-04-13 13:19:35

1.8K0

发布于 2020-04-13 13:19:35

在PCA(Principal Component Analysis)分析中，常用的工具有EIGENSOFT工具的smartpca，GCTA工具的PCA模块和R包中做PCA分析的princomp函数或glPCA功能。EIGENSOFT工具只支持linux系统，从安装到使用都很复杂。GCTA工具支持不同平台（wins/linux/mac）,常用于群体遗传相关分析。在群体遗传中，R包从读取vcf文件、PCA分析到可视化，对内存要求较高。

在这里我们主要介绍，针对测序得到的SNP数据（一般为vcf格式），如何利用GCTA工具进行PCA分析。以棉花的SNP数据为例，大体分析思路分为二进制转换、矩阵构建和可视化三个部分。

01

二进制转换

vcf文件格式一般如下：

CHROM POS ID REF ALT QUAL FILTER FORMAT L044A01 940519 7957 T G 23450.75 . GT:AD:DP:GQ:PL 0/0:2,0:2:6:0,6,85A01 940677 7958 A T 240.32 . GT:AD:DP:GQ:PL 0/0:3,0:3:9:0,9,117A01 940677 7959 A AT 6882.34 . GT:AD:DP:GQ:PL 0/0:3,0:3:9:0,9,97A01 940693 7960 G T 22.68 . GT:AD:DP:GQ:PL 0/0:3,0:3:9:0,9,119A01 940773 7961 A C 15095.96 . GT:AD:DP:GQ:PL 0/0:5,0:5:12:0,12,170A01 940958 7962 G A 20460.04 . GT:AD:DP:GQ:PL 0/0:9,0:9:24:0,24,343A01 941624 7963 A G 44831.97 . GT:AD:DP:GQ:PL 0/0:6,0:6:18:0,18,253A01 941652 7964 A G 46031.62 . GT:AD:DP:GQ:PL 0/0:5,0:5:12:0,12,164A01 941823 7965 T C 47012.52 . GT:AD:DP:GQ:PL 0/0:4,0:4:12:0,12,149A01 941953 7966 G A 19960.09 . GT:AD:DP:GQ:PL 0/0:3,0:3:9:0,9,128

在FORMAT一列中，GT表示样本基因型，AD表示覆盖到REF和ALT上碱基的read数；DP表示覆盖到该位点的总read数，GQ表示最可能基因型的质量值，PL表示该位点基因型为0/0，0/1，1/1的似然值L，L越大，说明该基因型的可能性越小。在LO44一列中，0表示跟REF一样，1表示跟ALT一样，2表示跟第二个ALT一样，./.表示未检测到SNP。

输入的vcf文件，如果是单一的染色体文件，需要用bcftools进行合并，具体过程如下：

1. 用bgzip工具将 vcf 文件压缩成 gz 文件

bgzip -c chr1.vcf > chr1.vcf.gz

2. 用bcftools工具为 gz 文件建索引

bcftools index chr1.vcf.gz

3. 利用bcftools工具将所有vcf文件合并

bcftools merge chr1.vcf.gz chr2.vcf.gz chr3.vcf.gz > chr123.vcf.gz

利用vcftools和plink工具将vcf格式文件转换成二进制文件。plink支持各种格式之间的转换，常见格式类型有：

一般格式（PED/MAP）转置格式（TPED/TFAM）二进制格式（BED/BIM/FAM）

bed文件包含SNP数据，bim文件包含SNP位置信息，fam文件包含表型信息。

1. 用vcftools做格式转换

##--plink输出plink可处理的文件格式vcftools --vcf A01.vcf --plink --out A01

生成.map和.ped（.ped文件具体信息可查看单倍型分析软件Haploview的导入格式及使用）

A01.pedA01.map

2. 用plink转换成二进制文件（输入和输出文件不需要加后缀名）

plink --noweb --file A01 --make-bed --out A01_bfile

生成.bed、.bim 和 .fam 3个文件

A01_bfile.bedA01_bfile.bimA01_bfile.fam

02

矩阵构建

利用GCTA软件构建用于PCA分析，主要分为两步：

1

计算近交系数，生成grm矩阵（输入和输出文件不需要加后缀名）

##--bfile读取二进制文件##--make-grm生成grm矩阵##--autosome指常染色体上所有SNP数据./gcta64 --bfile A01_bfile --make-grm --autosome --out A01_grm

生成三个文件

A01_grm.grm.bin  A01_grm.grm.id   A01_grm.grm.N.bin

2

对grm矩阵进行PCA分析

##--grm读取grm矩阵，--pca确定主成分个数./gcta64 --grm A01_grm --pca 3 --out A01_pca

生成两个文件

A01_pca.eigenvalA01_pca.eigenvec

在A01_pca.eigenvec文件加入表头：id ID pc1 pc2 pc3（这里视前面pca分组数而定）。该文件格式如下：

id   ID          PC1         PC2          PC3L044 L044  0.003068930 -0.05424570 -0.081746800D081 D081  0.002805640 -0.03164370 -0.017718400D069 D069  0.002644090  0.03152530  0.024485100B079 B079  0.002628090  0.01081540  0.004666700

03

可视化

利用R对PCA结果进行可视化画图，设置好工作路径，开始画图：

##读取matrix文件a1 <- read.table("A01_pca.eigenvec", header=T)##绘制散点图##pch取值1的时候代表空心圆，取值2的时候代表上三角plot(a1$PC1, a1$PC2, pch=c(rep(1,150), rep(2,109)), col=c(rep("blue", 150), rep("red", 109)), main="PCA", xlab="pc1", ylab="pc2")##添加图例legend("bottomright", c("grp1","grp2"), pch=c(rep(1), rep(2)), col=c(rep("blue"), rep("red")))