GWAS流程

QC质控

格式转化

• 初始文件为二进制类型文件bed、bim和fam，需要使用plink转化为ped、map格式

# temp是文件前缀，是一个临时文件
plink --bfile HapMap_3_r3_1 --recode --out temp

• 查看SNP个数和样本个数

# test.map是SNP信息， test.ped是样本信息
wc -l test.map test.ped

缺失质控

查看是否有缺失

plink --file test --missing

• 单个SNP缺失的个体数在plink.lmiss.中
• 个体缺失位点的统计在plink.imiss中

个体缺失质控

• 过滤在总样本SNP中缺失率大于2%的SNP（SNP在部分样品中缺少分型）

plink --bfile your_bed_file --geno 0.02 --make-bed --out filetemp

• 使用下面的命令进行验证，检验是否过滤

plink --bfile filetemp --recode --out temp
wc -l temp.map temp.ped

SNP缺失质控

• 过滤样本中总SNP缺失率大于2%的个体（样本质量不合格）

plink --bfile filetemp --mind 0.02 --make-bed --out filetemp2

• 使用下面的命令进行验证，检验是否过滤

plink --bfile filetemp2 --recode --out test
wc -l temp.map temp.ped

性别质控

主要是在人类研究中使用

• 检查性别冲突

plink --bfile filetemp2 --check-sex

• 生成plink.sexcheck文件
• 查看有问题数据的信息

grep "PROBLEM" plink.sexcheck

• 将相关错误的ID提取出来（家系ID，个体ID）

grep 'PROBLEM' plink.sexcheck | awk '{print $1,$2}'
>sex_discrepancy.txt

• 使用remove去掉个体

plink --bfile filetemp2 --remove sex_discrepancy.txt --make-bed --out filetemp3

MAF质控

• MAF < 0.02时，代表SNP代表的基因型分型比较极端（如b极少，B极多），我的理解是变异的基因所在的群体在自然选择下不占优势，或者虽是优良突变但是该突变的群体繁衍次数较少等，此时MAF呈假阳性，该SNP对性状的解释基本没贡献，需要过滤掉
• 以人的染色体为例，选取1~22号常染色体，提取第二列的ID

awk '$1 >= 1 && $1 <= 22 { print $2 }' filetemp2.bim > snp_1_22.txt

• 查看常染色体上的样品SNP个数

wc -l snp_1_22.txt

• 提取常染色体上的位点

plink --bfile filetemp2 --extract snp_1_22.txt --make-bed --out filetemp3

• 共有165个基因型，共有1398544个SNP

• 去掉MAF小于0.05的位点

plink --bfile filetemp3 --maf 0.05 --make-bed --out filetemp4

• 325318个位点被删掉了，剩余1073226个位点

哈温质控

• 计算所有位点的HWE的P值

plink --bfile filetemp3 --hardy

• 设定过滤标准1e-4 （一般为1e-4或1e-6， 若样本量很小，则设为1e-4或更低，避免过度过滤）

plink --bfile filetemp3 --hwe 1e-4 --make-bed --out filetemp4

杂合度质控（大规模饲养的缢蛏是否需要？）

排除近亲繁殖、样品污染、技术误差等因素

• 基因型个体的杂合度过高或者过低，都不正常，我们需要根据杂合度进行过滤。偏差可能表明样品受到污染，近亲繁殖。我们建议删除样品杂合率平均值中偏离±3 SD的个体。
• 计算杂合度

plink --bfile filetemp4 --het --out R_check

• 查看在3倍标准差以外的个体

Rscript heterozygosity_outliers_list.R

# heterozygosity_outliers_list.R
# 先要设置工作路径
het = read.table("R_check.het", head=TRUE)
het$HET_RATE = (het$"N.NM." - het$"O.HOM.")/het$"N.NM."
het_fail = subset(het, (het$HET_RATE < mean(het$HET_RATE)-3*sd(het$HET_RATE)) | (het$HET_RATE > mean(het$HET_RATE)+3*sd(het$HET_RATE)));
het_fail$HET_DST = (het_fail$HET_RATE-mean(het$HET_RATE))/sd(het$HET_RATE);

write.table(het_fail, "fail-het-qc.txt", row.names=FALSE)

• 去掉这些个体：先对数据进行清洗，去掉引号，然后提取家系和个体ID

sed 's/"//g' fail-het-qc.txt |awk '{print $1,$2}' > het_fail_ind.txt

• 使用remove去掉这两个个体

plink --bfile filetemp4 --remove het_fail_ind.txt --make-bed --out filetemp5

亲缘关系质控

• 这里讲亲子关系的个体移除，不是必须要的，比如我们分析的群体里面有亲

  子关系的个体，想要进行分析，不需要做这一步的筛选。

plink --bfile HapMap_3_r3_10 --filter-founders --make-bed --out HapMap_3_r3_11

• 52个被移除

质控结果

• 查看经上面质控后的结果

plink --bfile filetemp4 --recode --out test

wc -l test*

PCA分析

• 使用GCTA对PCA分析

# 使用yang的方法计算遗传关系矩阵
gcta --bfile filetemp4 --make-grm --make-grm-alg 0 --out kinship_yang
# 使用遗传关系矩阵进行主成分分析
gcta --grm kinship_yang --pca 20 --out pca_re

• 协变量文件和PCA文件合并

# 剪切PCA文件
awk '{print $3,$4,$5}' pca_re.eigenvec > pca.txt
# 将剪切的部分粘贴到协变量文件中
paste cov.txt pca.txt | sed 's/[[:space:]]\+/ /g' > cov2.txt

• 可以查看列数，检验是否合并了

head -n 1 cov2.txt | awk '{print NF}'

GEMMA分析

• 将表型数据单独提取一列

awk '{print $3}' PHE1.phe > p.txt

• 生成G矩阵

# 参数gk一般选1
gemma -bfile filetemp4 -gk 2 -p p.txt

• 使用MLM模型进行分析

# 注意这里的参数p后面是上面提取的一列表型数据
gemma -bfile filetemp4 -k output/result.sXX.txt -lmm 1 -p p.txt -c test2.covar

plink --bfile HapMap_3_r3_9 --recode --out fileread
plink --file fileread --recode HV --out fileinfo