绘制曼哈顿图的时候，需要指定一个显著的阈值线
常规的阈值是0.05/n或0.01/n。n是基因型的标记数量（snp的数量）。
但是因为存在连锁不平衡，很多时候按照上面的这个阈值，会筛不出来位点。这种情况称为假阴线，所以就需要调整阈值。
调整阈值的软件有：GEC或simpleM.
安装使用GEC
从http://pmglab.top/gec/#/download下载kggsee.jar
java -Xmx10g -jar ~/soft/GEC/kggsee.jar --var-gec --nt 12 --filter-maf-le 0.02 --chrom-col chr --pos-col pos --p-col p_value --p-value-cutoff 0.05 --vcf-ref ~/blink/demo_data/vcf/myData.vcf --pfile EarDia_GWAS_result.txt --out test
参数讲解：
--var-gec 必备参数
--nt 12 #12个cpu
--filter-maf-le 0.02 #maf的过滤阈值,默认是0.05
--chrom-col chr #指定pfile里的染色体的列名，默认是CHR
--pos-col pos #指定pfile里的位置列名，默认是BP
--p-col p_value#指定pfile里的Pvalue的列名，默认是P
--p-value-cutoff 0.05 #指定阈值，默认是0.05
--vcf-ref ~/blink/demo_data/vcf/myData.vcf #指定vcf文件的地址
--pfile EarDia_GWAS_result.txt#指定gwas的输出结果的文件，该文件必须包含至少3列，内容分别是染色体名称、变异位置、Pvalue
--out test#指定输出的文件的前缀
输出文件是：
test.log
test.effective.size.txt.gz
在test.log里有以下的内容：
3093 variant-lines (0 indels) are scanned in ~/soft/blink/demo_data/vcf/myData.vcf; and 2767 variants of 281 individual(s) are retained.
...
INFO 2023-05-10 11:28:32 - Variants: Significance level of p value cutoffs for the overall error rate 0.05 by Standard Bonferroni : 1.616553507921112E-5
...
The p-value cutoff by Bonferroni correction for family-wise error rate 0.05 is 9.34E-04
最终使用的P-value是 9.34E-04。FWER
simpleM
https://simplem.sourceforge.net/
https://yaoulab.group/blog/2022/The-effect_number-of-SNPs/
https://github.com/LTibbs/SimpleM
使用simpleM之前需要先对缺失基因型进行填充或过滤，同时过滤maf<0.05的次等位基因型
填充缺失基因型使用beagle beagle的下载地址
java -jar ~/soft/beagle/beagle.22Jun22.46e.jar gt=test_final.vcf out=test.imputed
beagle输出的文件是test.imputed.vcf.gz
过滤MAF<0.05使用plink
plink --vcf test.imputed.vcf.gz --maf 0.05 --geno 0.1 --recode vcf-iid --out test.filter --allow-extra-chr
一般情况下，进行GWAS分析的vcf文件都已经完成了上面的填充和过滤，所以直接计算numeric矩阵即可。
从vcf获取基因型的矩阵，矩阵里只能有0,1,2.
``
使用SimpleM计算获取最终的Pvalue值。
Rscript ~/soft/SimpleM/SimpleM.R genotype_transpose.txt
SimpleM脚本最常见的错误是：
Error in eigen(CLD) : 'x'里有无穷值或遗漏值
In addition: Warning message:
In cor(dt_My) : 标准差为零
解决办法是：过滤vcf文件的maf<0.05或0.01
下面是我写的使用SimpleM.R自动计算P的阈值的脚本。
使用前先修改软件路径为你的安装路径
#!/bin/bash
##自动分析计算最佳P的阈值的脚本
input_vcf=$1
abbr=$2
##软件路径配置
beagle="java -jar ~/soft/beagle/beagle.22Jul22.46e.jar"
plink="~/soft/LDBlockShadow/LDBlockShow-1.36/bin/plink"
SimpleM="Rscript ~/soft/SimpleM/SimpleM.R"
kggsee="java -Xmx10g -jar ~/soft/GEC/kggsee.jar"
#使用beagle对vcf文件进行填充
$beagle gt=${input_vcf} out=${abbr}.beagle
#输出文件是${abbr}.beagle.vcf.gz
#使用plink过滤输入的数据
$plink --vcf ${abbr}.beagle.vcf.gz --maf 0.05 --geno 0.1 --recode vcf-iid --out ${abbr}.filter
#输出文件是${abbr}.filter.vcf
#把vcf文件转为numeric格式的矩阵只有0,1,2
cat ${abbr}.filter.vcf | grep -v "^##" | cut -f 10- | sed 's/0\/0/0/g' | sed 's/1\/1/2/g' | sed 's/0\/1/1/g' | sed 's/1\/0/1/g'| sed 's/\.\/\./-1/g' | tr "\t" " "|sed '1d' >${abbr}.genotype_transpose.txt
##使用SimpleM.R计算P的阈值
${SimpleM} ${abbr}.genotype_transpose.txt 0.05 >${abbr}.SimpleM.pvalue
将上述脚本保存为getPvalue.bash.
bash getPvalue.bash vcffile abbr #参数1是你的vcf文件，参数2是输出前缀字符串
输出的P值在abbr.SimpleM.pvalue的第三行。
我修改了原版SimpleM.R脚本，从Simple github可以下载。
上述脚本也保存为getPvalue.bash
