从零开始学GWAS之数据格式

从零开始学GWAS之数据格式
我最近正好在学习GWAS（全基因组关联分析），这⼀套流程呢简单来说就是在全基因组范围内寻找变异序列，在全基因组范围内进⾏相关性分析，通过⽐较发现复杂性状的基因变异。

我学习的过程中参考的是黄学辉教授15年发表在NC上的⽂章“Genomic analysis of hybrid rice varieties reveals numerous superior alleles that contribute to heterosis”。

⾥⾯介绍到分析过程中运⽤到的⼀些软件和⽅法。

原始数据的处理，变异检测（BWA+SAMtools+picard+GATK）这套流程进⾏SNP calling⽣成VCF⽂件我就不介绍了。

在介绍GWAS之前，有⼏个分析过程中产⽣的⽂件，有⼏个格式需要先稍微介绍⼀下。

变异信息存放格式之VCF
VCF是⽤于描述SNP，INDEL，SV的⽂本⽂件。

是GATK表⽰遗传变异的⼀种⽂件格式。

可以分成两个部分来看：
第⼀部分是以##开头的说明⽂件，解释第⼆部分INFO列中可能要出现的⼀些tags和和FORMAT列中对基因型的表⽰。

个⼈觉得最个注释没什么⽤处，我基本都是直接跳过的。

##fileformat=VCFv4.1##FILTER=##FORMAT=##FORMAT=##FORMAT=##INFO=##INFO=##contig=##contig=##contig=##contig=
第⼆部分是重点了，正⽂部分主要包括9列+N样品列
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT B73 F2-Mo17vsB73 Mo17 1 530 . C G 52.17 .
AC=2;AF=0.333;AN=6;BaseQRankSum=-
0.948;DP=535;Dels=0.00;FS=14.393;HaplotypeScore=11.0860;MLEAC=1;MLEAF=0.167;MQ=6.29;MQ0=396;MQRankSum=2.281;QD=0.15;ReadPosRankSum=0.530;SOR=3.223 GT:AD:DP:GQ:PL 0/1:208,3:218:13:13,0,100 0/0:176,0:177:39:0,39,297 0/1:136,2:140:31:45,0,31 1 534 . G A 32.35 .
AC=1;AF=0.167;AN=6;BaseQRankSum=-
0.117;DP=539;Dels=0.00;FS=10.307;HaplotypeScore=15.3371;MLEAC=1;MLEAF=0.167;MQ=6.36;MQ0=397;MQRankSum=3.126;QD=0.15;ReadPosRankSum=0.154;SOR=1.431 GT:AD:DP:GQ:PL 0/1:206,4:214:63:63,0,117 0/0:178,2:182:33:0,33,252 0/0:139,0:143:12:0,12,101 1 542 . C T 32.35 .
AC=1;AF=0.167;AN=6;BaseQRankSum=-
1.405;DP=534;Dels=0.00;FS=11.442;HaplotypeScore=1
2.8859;MLEAC=1;MLEAF=0.167;MQ=6.38;MQ0=391;MQRankSum=2.054;QD=0.15;ReadPosRankSum=-
0.330;SOR=2.221 GT:AD:DP:GQ:PL 0/1:207,10:218:63:63,0,117 0/0:175,3:178:39:0,39,297 0/0:134,3:138:12:0,12,101
CHROM和pos：表⽰变异位点相对reference的位置，⽐如第⼏条染⾊体的第⼏个碱基，如果是indel，pos是indel的第⼀个碱基的位置
ID：如果call出来的SNP存在于dbsnp数据库⾥，就会显⽰相应的dbsnp⾥的rs编号。

不然就是⽤“.”表⽰⼀个novel variant.
REF和ALT：分别代表reference和alter,也就是参考基因组对应的碱基和variant的碱基。

QUAL：表⽰该位点存在variant的可能性，qual值越⼤则variant的可能性越⼤。

FILTER：过滤完了之后，FILTER⼀栏会有过滤记录，通过了过滤标准，那么这些好的变异位点的FILTER⼀栏就会注释⼀个PASS，如果没有通过过滤，就会在FILTER这⼀栏提⽰其他信息。

如果这⼀栏是⼀个“.”的话，就说明没有进⾏过任何过滤。

IFNO：这⼀列表⽰的是variant的详细信息。

GT：表⽰样本的基因型，对于⼀个⼆倍体⽣物，GT值表⽰的是这个样本在这个位点所携带的两个等位基因。

0表⽰跟REF⼀样；1表⽰表⽰跟ALT⼀样；2表⽰第⼆个ALT。

当只有⼀个ALT等位基因的时候，0/0表⽰纯和且跟REF⼀致；0/1表⽰杂合，两个allele⼀个是ALT⼀个是REF；1/1表⽰纯和且都为ALT；
AD：对应两个以逗号隔开的值，这两个值分别表⽰覆盖到REF和ALT碱基的reads数，相当于⽀持REF和⽀持ALT的测序深度。

DP：覆盖到这个位点的总的reads数量，相当于这个位点的深度（并不是多有的reads数量，⽽是⼤概⼀定质量值要求的reads数）。

PL:对应3个以逗号隔开的值，这三个值分别表⽰该位点基因型是0/0，0/1，1/1的没经过先验的标准化Phred-scaled似然值（L）。

如果转换成⽀持该基因型概率（P）的话，由于L=-10lgP，那么P=10^（-L/10）^，因此，当L值为0时，P=10^0^=1。

因此，这个值越⼩，⽀持概率就越⼤，也就是说是这个基因型的可能性越⼤。

GQ：表⽰最可能的基因型的质量值。

表⽰的意义同QUAL。

GQ：表⽰最可能的基因型的质量值。

表⽰的意义同QUAL。

变异信息存放格式之 tped
tped格式：是varient信息+基因型调⽤的⽂本⽂件
tped⽂件时没有header的。

每⼀⾏有4+2N个条⽬，N代表样本数⽬。

第⼀列是染⾊体号
第⼆列是SNP的ID
第三列表⽰单位可以是厘摩或者摩尔根（也可以⽤虚拟值0）
第四列表⽰SNP的位置
后⾯的都表⽰样本列，⼀个⼆倍体的样本会有两列。

两个数分别表⽰等位基因，可以⽤0、1；1、2表⽰。

’12’修饰符表⽰A1等位基因编码为’1’，A2等位基因编码为’2’，⽽’01’则映射A1→0和A2→1。

变异信息存放格式之 tfam
tfam格式: tfam⽂件和tped⽂件相对应，tfam⽂件样品是纵向排列⽽tped（第五列开始）横向排列。

tfam⽂件每⾏包括6个条⽬
第⼀列是样本ID（FID）
第⼆列是样本内部ID（IID;不能为‘0’）
第三列⽗本ID（如果⽗本不在数据集⾥⽤“0”表⽰）
第四列母本ID（如果母本不在数据集⾥⽤“0”表⽰）
第五列表⽰性别代码：（’1’=男性，’2’=⼥性，’0’=未知，也可以⽤其他数字来凑格式）
第六列是表型数据
表型数据
最后还有⼀个表型的输⼊⽂件，⼀般会⽤.pheno的后缀表⽰
表型⽂件在每⼀⾏有三个条⽬，FAMID，INDID和表型值。

缺失的表型值应表⽰为“NA”。

（偷懒⼀点就可以直接⽤tfam⽂件⾥的地1、2、6列）参考资料
/s/blog_12d5e3d3c0101qv1u.html。