繪制曼哈頓圖的時(shí)候，需要指定一個(gè)顯著的閾值線
常規(guī)的閾值是0.05/n或0.01/n具帮。n是基因型的標(biāo)記數(shù)量（snp的數(shù)量）。
但是因?yàn)榇嬖谶B鎖不平衡匪凡，很多時(shí)候按照上面的這個(gè)閾值掘猿，會(huì)篩不出來位點(diǎn)。這種情況稱為假陰線衬衬，所以就需要調(diào)整閾值改橘。
調(diào)整閾值的軟件有：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
參數(shù)講解：
--var-gec 必備參數(shù)
--nt 12 #12個(gè)cpu
--filter-maf-le 0.02 #maf的過濾閾值,默認(rèn)是0.05
--chrom-col chr #指定pfile里的染色體的列名飞主，默認(rèn)是CHR
--pos-col pos #指定pfile里的位置列名，默認(rèn)是BP
--p-col p_value#指定pfile里的Pvalue的列名碾篡，默認(rèn)是P
--p-value-cutoff 0.05 #指定閾值丸冕，默認(rèn)是0.05
--vcf-ref ~/blink/demo_data/vcf/myData.vcf #指定vcf文件的地址
--pfile EarDia_GWAS_result.txt#指定gwas的輸出結(jié)果的文件，該文件必須包含至少3列，內(nèi)容分別是染色體名稱诅迷、變異位置、Pvalue
--out test#指定輸出的文件的前綴
輸出文件是：
test.log
test.effective.size.txt.gz
在test.log里有以下的內(nèi)容：
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之前需要先對(duì)缺失基因型進(jìn)行填充或過濾赋秀，同時(shí)過濾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
一般情況下律想，進(jìn)行GWAS分析的vcf文件都已經(jīng)完成了上面的填充和過濾，所以直接計(jì)算numeric矩陣即可著洼。
從vcf獲取基因型的矩陣，矩陣?yán)镏荒苡?,1,2.
``
使用SimpleM計(jì)算獲取最終的Pvalue值豹悬。
Rscript ~/soft/SimpleM/SimpleM.R genotype_transpose.txt
SimpleM腳本最常見的錯(cuò)誤是：
Error in eigen(CLD) : 'x'里有無窮值或遺漏值
In addition: Warning message:
In cor(dt_My) : 標(biāo)準(zhǔn)差為零
解決辦法是：過濾vcf文件的maf<0.05或0.01
下面是我寫的使用SimpleM.R自動(dòng)計(jì)算P的閾值的腳本液荸。
使用前先修改軟件路徑為你的安裝路徑
#!/bin/bash
##自動(dòng)分析計(jì)算最佳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對(duì)vcf文件進(jìn)行填充
$beagle gt=${input_vcf} out=${abbr}.beagle
#輸出文件是${abbr}.beagle.vcf.gz
#使用plink過濾輸入的數(shù)據(jù)
$plink --vcf ${abbr}.beagle.vcf.gz --maf 0.05 --geno 0.1 --recode vcf-iid --out ${abbr}.filter
#輸出文件是${abbr}.filter.vcf
#把vcf文件轉(zhuǎn)為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計(jì)算P的閾值
${SimpleM} ${abbr}.genotype_transpose.txt 0.05 >${abbr}.SimpleM.pvalue
將上述腳本保存為getPvalue.bash.
bash getPvalue.bash vcffile abbr #參數(shù)1是你的vcf文件娇钱，參數(shù)2是輸出前綴字符串
輸出的P值在abbr.SimpleM.pvalue的第三行。
我修改了原版SimpleM.R腳本忍弛，從Simple github可以下載响迂。
上述腳本也保存為getPvalue.bash
