论文

A new regulator of seed size control in Arabidopsis identified by a genome-wide association study
New Phytologist 2019
Peking University

两篇简书文章重复论文中GWAS分析的过程

• 重复一篇文献的GWAS（一）：基因型数据整理
• 重复一篇文献的GWAS（二）：用GEMMA跑GWAS

重点关注论文中GWAS的分析过程，争取根据两篇简书文章重复出分析过程

测量种子大小 （Seed size measurement）

种子平铺白纸，扫描，使用imageJ软件进行分析

全基因组关联分析

• 表型数据（种子面积）用平均值
• SNP数据下载自1001 Genomes Project website
• BeAGLE软件填充数据
• 数据过滤
  非二等位基因位点
  最小等位基因频率0.05
• 数据过滤
  plink连锁不平衡
• 单变量混合线性模型
  GEMMA软件默认参数
• 后面两句话看不明白
  A linear mixed univariate model was used for marker association tests with a single phenotype to account for population stratification and sample structure.
  The population structure was controlled using the relatedness matrix generated from GEMMA with the -gk 1 parameter.
• 校正p值
  小于0.05显著
• 岭回归分析（？）

分析过程

下载数据

wget http://1001genomes.org/data/GMI-MPI/releases/v3.1/1001genomes_snp-short-indel_only_ACGTN.vcf.gz

初步筛选数据

vcftools --gzvcf 1001genomes_snp-short-indel_only_ACGTN.vcf.gz --remove-indels --min-alleles 2 --max-alleles 2 --recode --keep At_Seed_Size_Samples.txt --out 172sample

• 删除InDel
• 只保留二等位位点
• 挑选指定的样本

基因型填充

java -Xss5m -Xmn25G -Xms50G -Xmx50G -jar ~/mingyan/Bioinformatics_tool/Beagle/beagle.25Nov19.28d.jar nthreads=4 gt=172sample.recode.vcf out=172sample_out ne=172

遇到了报错

Caused by: java.lang.OutOfMemoryError: GC overhead limit exceeded

先跳过这一步

根据最小等位基因来过滤

vcftools --vcf 172sample.recode.vcf --maf 0.05 --recode --out 172sample_maf_filter

给每个snp位点添加ID

bcftools view -h 172sample_maf_filter.recode.vcf > head
bcftools view -H 172sample_maf_filter.recode.vcf | grep "^#" -v | awk '{$3=$1"_"$2;print $0}' | sed 's/ /\t/g' > body
cat head body > 172sample_maf_filter_snpID.vcf

到这一步数据集已经小了很多，试着在这一步对基因型进行填充（先填充在过滤和先过滤再填充的区别？）
基因型填充

java -jar beagle.25Nov19.28d.jar gt=172sample_maf_filter_snpID.vcf out=172sample_out ne=172

基因型填充需要很长时间，这里先用没有填充的做接下来的分析

基于连锁不平衡进行过滤

plink --vcf 172sample_maf_filter_snpID.vcf --recode --out 172sample
plink --file 172sample --indep 50 5 2
plink --file 172sample --extract plink.prune.in --recode --out 172sample_maf_filter_snpID_LD_filter

ped/map转换为fam/bed/bim

plink --file 172sample_maf_filter_snpID_LD_filter --make-bed --out clean_snp

为结果文件添加表型数据，需要在fam文件最后一列（-9）修改为真实表型值，这里我使用https://www.jianshu.com/p/fc628bd1001b 中提到的第二种方法
代码

import pandas as pd
tableA = pd.read_excel("add_Pheno_to_fam.xlsx",sheet_name="Sheet3",converts={'A':str})
tableB = pd.read_excel("add_Pheno_to_fam.xlsx",sheet_name="Sheet2",converts={'A':str})
tableA.head()
tableB.head()
tableA.merge(right=tableB,how="left",left_on="A",right_on="A")
tableC = tableA.merge(right=tableB,how="left",left_on="A",right_on="A")
tableC.to_excel("C.xlsx",index=False)

gemma软件下载
参考文章 使用GEMMA进行复杂性状全基因组关联分析(GWAS)
代码

gemma-0.98.1-linux-static -bfile clean_snp -gk -o kinship
gemma-0.98.1-linux-static -bfile clean_snp -lmm -k ./output/kinship.cXX.txt -o GWAS_results

结果文件是GWAS_results.assoc.txt，用这个文件来画图
代码

rawdf<-read.table("GWAS_results.assoc.txt",
                  header=T,sep="\t")
dim(rawdf)
colnames(rawdf)
df<-data.frame(rs=rawdf$rs,chr=rawdf$chr,pos=rawdf$ps,
               pvalue=rawdf$p_wald)
head(df)
install.packages("rMVP")
library(rMVP)
?MVP.Report
MVP.Report(df)
MVP.Report(df,plot.type="m")

Rectangular-Manhattan.pvalue.jpg
QQplot.pvalue.jpg Circular-Manhattan.pvalue.jpg SNP_Density.pvalue.jpg

这个rMVP包好厉害，一条命令MVP.Report(df)所有图就全部生成了

文章开头提到的两篇简书文章还用rMVP这个R语言包做了全基因组关联分析，另外找时间来重复这部分内容
还有一部分内容是用R语言的glmnet包做岭回归分析分析，这部分内容暂时还没有思路！

参考文章

• 1、 重复一篇文献的GWAS（一）：基因型数据整理
• 2、 重复一篇文献的GWAS（二）：用GEMMA跑GWAS
• 3、 使用GEMMA进行复杂性状全基因组关联分析(GWAS)

本篇文章绝大部分分析代码均来自参考文章1和2，特别感谢文章1和2作者的分享

欢迎大家关注我的公众号
小明的数据分析笔记本
如果需要文章中任何步骤中用到的数据，欢迎在我的公众号留言

公众号二维码.jpg