按照前人的教程，跑了跑GWAS流程，作为初学者，欢迎大家提问，指教。
数据来源：A new regulator of seed size control in Arabidopsis identified by a genome-wide association study
做的是拟南芥191个自交系的籽粒大小的全基因组关联分析，并且对其中一个显著的qtl进行了后续的实验验证。
基因组数据下载
wget http://1001genomes.org/data/GMI-MPI/releases/v3.1/1001genomes_snp-short-indel_only_ACGTN.vcf.gz
表型数据下载：
文章的supplement的pdf文件，可以用R包读取其中三线表中的内容。
vcftools --gzvcf 1001genomes_snp-short-indel_only_ACGTN.vcf.gz \
--remove-indels --min-alleles 2 --max-alleles 2 \
--recode --keep ../sample/sample.txt --out 172sample
#这个vcf文件不只包含文献里面的自交系基因组数据，首先把文献里面的数据基因组数据调出来。remove-indels,min- or max-alleles至少至多多少个等位位点，--recode输出.recode.vcf,--keep file保留file中的自交系，--out输出前缀
#其中sample.txt是表型数据中easyGWAS ID那一列，一行一个ID，有19个材料没有ID，只能用有ID的172个材料。
vcftools --vcf 172sample.recode.vcf --freq --out freq
#--freq,输出等位基因频率，看一下snp基因频率，发现有些snp是的MAF为零。说明这些位点可能在之前的大群体里面是有变异，而在该研究用的这172个自交系里面没有变异。
module load Java/1.8.0_92
module load beagle-lib/2.1.2-foss-2016b
java -Xss5m -Xmx100g -jar $EBROOTBEAGLE/beagle.jar \
nthreads=20 gt=172sample.recode.vcf out=172sample_out ne=172
#基因型填充，从这一步的输出文件来看，SNP个数和前面第一步提取出来的SNP的个数是相同的，说明这一步的作用是：比如有一个SNP在某些材料里面没有被测到，没有基因型，这一步通过单倍型给这些没有分型的位点加上了基因型。
#Xmx*g,用多大内存 nthreads核心数，gt输入的vcf文件，ne=群体大小
vcftools --gzvcf 172sample_out.vcf.gz --maf 0.05 \
--recode --out 172sample_maf_filter
#maf=0.05,maf小于0.05的snp被过滤掉，至此SNP的第一步过滤基本完成。
#后面的步骤主要是要过滤SNP，让标记各自连锁平衡。否则，与真正显著的标记连锁的一大片标记，都会有显著关联。
awk '{FS="\t";OFS="\t"}{if($0~"#"){print $0}else{$3=$1"_"$2;print $0}}' \
172sample_maf_filter.recode.vcf > 172sample_maf_filter_snpID.vcf
#因为后续plink软件的需要，给SNP加上ID。plink中默认用下划线分隔，下划线后面的是SNP的ID号。
plink --vcf 172sample_maf_filter_snpID.vcf --recode --out 172sample
#转成plink格式，生成172sample.log,172sample.map,172sample.ped,172sample.nosex
#其中map和ped参数是基因组和表型数据
plink --file 172sample --indep 50 5 2
#indep是prune时的参数，不太懂模型，第一个是滑动窗口大学，第二个是步长，第三个好像是一个判断是否连锁的阈值吧。
#产生prune.in，prune.out，prune.in文件基本连锁平衡的SNP位点，但是很奇怪，相同的数据和分析流程，我这一步就只剩下了256233个标记了，和我参考的博客差很多，但是没prune之前的标记数和这位师兄做得差不多，我做的是maf过滤后剩下1601218个标记。
plink --file 172sample --extract plink.prune.in \
--recode --out 172sample_maf_filter_snpID_LD_filter
#file输入前缀,--extract根据prune.in提取snp，172sample_maf_filter_snpID_LD_filter.*就是我们需要的数据。但是我们需要二进制的格式，用于后续的关联分析。
plink --file 172sample_maf_filter_snpID_LD_filter \
--make-bed --out clean_snp
#--file,输入文件前缀,--make-bed --out会生成以clean_snp为前缀的bed,bim和fam
#fam存有表型数据,要在fam的第六列加入对应的表型数据
这里我自己写了个python脚本，用来将fam文件和supplement读出来的表格聚合在一起，然后输出一个对应行加上了对应表型值的fam文件。
python3 mergefam_sample.py
###########这是mergefam_sample.py代码
import pandas as pd
sample = pd.read_table("../../sample/sample_list.csv",sep=",")
sample = sample.iloc[:,[1,5]]
#sample有191个样本，easyGWAS这一列为字符
fam = pd.read_table("clean_snp.fam",sep=" ",header=None)
fam.columns = ["ID1","ID2","parent1","parent2","sex","phenotype"]
fam["ID1"] = fam["ID1"].apply(str)
#fam有172个样本，前两列为数字
new_fam = pd.merge(fam,sample,how="inner",left_on="ID1",right_on="easyGWAS ID")
new_fam.phenotype = new_fam["Seed size (mm2)"]
new_fam = new_fam.iloc[:,range(0,6)]
new_fam = new_fam.round(decimals=9)
new_fam.to_csv("clean_snp.fam1",sep=" ",header=False,index=False)
###########这是mergefam_sample.py代码
cp clean_snp* ../../associate
cd clean_snp* ../../associate
rm clean_snp.fam
mv clean_snp.fam1 clean_snp.fam
#转移文件到关联分析的目录下
module load GEMMA/0.98.1
gemma-0.98.1-linux-static -bfile clean_snp -gk 1 -o kinship
#-gk 1或者2,不同计算亲缘关系的方法，一般都用1，官方文档给的说明是，如果一个SNP的maf比较小但是效应值比较大的用2会比较好，如过maf和该SNP的效应值无关就用1。
#-bfile输入为plink的二进制，-o输出前缀,输出一个output的目录.cXX.txt就是亲缘关系的结果
gemma-0.98.1-linux-static -bfile clean_snp -lmm -k ./output/kinship.cXX.txt -o GWAS
#-lmm线性混合模型可以等于1,2,3和4，代表不同的假设检验方法(默认参数为1,文章用的默认参数),-k亲缘关系结果文件，.assoc.txt是结果是需要的结果
#ps SNP物理位置，n_miss SNP缺失个体数，af SNP频率，beta SNP效应值，se beta估计标准误，l_remle 计算该SNP效应时对应的lamda的remle估计值，p_wald wald检验P值。就用最后一列的pvalue挑显著的标记，阈值是0.05/prune之后的标记数。
#最后画图，用r包rMVP，简单好用不愧MVP
rawdf<-read.table("GWAS_results.assoc.txt",header=T,sep="\t")
df<-data.frame(rs=rawdf$rs,chr=rawdf$chr,pos=rawdf$ps,pvalue=rawdf$p_wald)
install.packages("rMVP")
library(rMVP)
MVP.Report(df)
#生成了四张图，qq图，普通和环状的曼哈顿图，标记密度图
参考的几个博客
https://cloud.tencent.com/developer/article/1593313
