文献笔记五十四：全基因组关联分析鉴定拟南芥中控制种子大小的调节因子

论文

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显著
岭回归分析（？）

分析过程

下载数据

代码语言：javascript复制

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

初步筛选数据

代码语言：javascript复制

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
只保留二等位位点
挑选指定的样本

基因型填充

代码语言：javascript复制

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

遇到了报错

代码语言：javascript复制

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

先跳过这一步

根据最小等位基因来过滤

代码语言：javascript复制

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

给每个snp位点添加ID

代码语言：javascript复制

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

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

代码语言：javascript复制

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

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

基于连锁不平衡进行过滤

代码语言：javascript复制

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

代码语言：javascript复制

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

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

代码语言：javascript复制

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) 代码

代码语言：javascript复制

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，用这个文件来画图代码

代码语言：javascript复制

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")