全基因组关联分析（GWAS）软件：mrMLM

mrMLM 官网介绍，它的优势是先用一个比 Bonferroni 略宽松的标准筛选SNP，然后用多位点遗传模型分析，可以找到更多的与目标性状关联的SNP。
不得不吐槽一下，使用体验太差，输入文件不是常用格式（vcf、ped等），也没有解释说明（说明书里的信息基本等于零），逼的我只好看程序代码如何读取数据的。。。；后来发现，R包的文件夹里有个原始数据的示例文件，终于有了点头绪。。。
ps：后来终于在bioRxiv上的关于mrMLM 4.0的文章的附件中找到一份比较详细的说明.

前言

输入文件有三种格式，numeric 格式、charactor格式、hapmap格式。
numeric 格式和charactor格式也跟通常理解的有所不同，并且必须经过imputation，不能有缺失基因型，因为源代码里没有识别缺失基因型数据

numeric 格式

rs#     chrom   pos     genotype for code 1    line1    line2    line3    ...
SNP1    1    100    C    1    0    -1

前四列分别是SNP，染色体，位置，1代表的基因型
注意，前四列列名不能变；1/-1表示纯合基因型，0表示杂合基因型和缺失型

charactor 格式

rs#     chrom   pos    line1    line2    line3    ...
SNP1    1    100    A    N    T

前列三列列名不能变；AGCT之外的字母表示杂合基因型及缺失

hapmap 格式

没有仔细看，应该是hapmap的正常格式，前几列信息的列名也不能变

构造输入文件

Genotype（numeric格式）

#csv格式
 bcftools query -f '%CHROM\_%POS,%CHROM,%POS,%REF[,%GT]\n' snp_hardfiltered_biallelic_maf5_ms50_imputated.vcf >num_format.csv
# 0|0 基因型改为1（即参考基因型）、1|1基因型改为-1、0|1或1|0基因型改为0
sed -i -r -e "s/0\|0/1/g" -e "s/1\|1/\-1/g" -e "s/(0\|1)|(1\|0)/0/g" num_format.csv
grep "#CHROM" snp_hardfiltered_biallelic_maf5_ms50_imputated.vcf|cut -f1-4,10-|xargs -I [] sed -i '1i []' num_format.csv
vi num_format.txt #手动改前四列列名
# mrMLM 内置函数 DoData():gen <- matrix(as.numeric(gen), nrow = nrow(gen))，在做矩阵转换的时候，如果chrom命名是字符类型，会强制转换为NA；因此，染色体命名只能是纯数字形式。（坑！说明文件里没有说明，也没有错误提示，只有个as.numeric（）的强制转换提醒，只能回溯代码一行行检查）
# R通过as.datatype(object)函数将对象类型进行强制转换，转换原则也遵循高等级对象可转换为低等级，而低等级对象不能转换为高等级。即数字可以转换为字符，字符不能转换为数字；因子转换为数字会按序以1,2,3...代替。

phenotype

表型文件phenotype：taxa、trait1、trait2…；第一列为材料号（表头为，之后为性状，列名不作要求,可以有多列。

<Phenotype>	trait1
B46	42
B57	72.5
…	…

计算structure 矩阵

用 admixture 计算群体结构矩阵，输入格式为plink，bed

#转为plink格式
plink --vcf ~/data/SLAFseq-rawdata/snp_hardfiltered_biallelic_maf5_ms50_imputated.vcf --make-bed --out SLAF --allow-extra-chr 
for K in `seq 1 5`;do admixture --cv 185SLAF.bed $K|tee log${K}.out;done

 #交叉验证，计算最小K值 
grep -h CV *.out
CV error (K=1): 0.78357
CV error (K=2): 0.73175
CV error (K=3): 0.69746
CV error (K=4): 0.68539
CV error (K=5): 0.67758
CV error (K=6): 0.66330
CV error (K=7): 0.66792

admixture SLAF.bed 6
cut -d ' ' -f1 SLAF.fam |paste -d ' ' - admixture/SLAF.6.Q |sed 's/ /,/g' >PopStr.csv

#加表头
<PopStr>,,,,,,
<ID>,Q1,Q2,Q3,Q4,Q5,Q6
S1,0.491380,0.401732,0.000010,0.000010,0.000010,0.106858
S10,0.000010,0.000010,0.000010,0.999950,0.000010,0.000010
S100,0.575853,0.171079,0.000010,0.000010,0.000010,0.253038
S101,0.081143,0.804484,0.114343,0.000010,0.000010,0.000010
S102,0.353952,0.254454,0.000010,0.000010,0.047595,0.343979

运行

R 语言

library('mrMLM')

# 建议直接在 mrMLM() 直接用genotype和phenotype的路径，mrMLM可以自动识别格式，如果用 readtable() 读取genotype的时候，一定要加入参数comment.char=""，因为第一列名为 rs#，“#”是默认的注释字符。

mrMLM(fileGen="num_format.txt",filePhe="SPAD.txt",fileKin=NULL,filePS="PopStr.csv",fileCov=NULL,Genformat="Num",method="mrMLM",trait=1,SearchRadius=100,CriLOD=3,DrawPlot=TRUE,Plotformat="jpeg",Resolution="Low",dir="./",PopStrType="Q")

kinship 矩阵(fileKin)和群体结构矩阵（fileSP）都可以设置为空（NULL），但是软件只会自动计算kinship矩阵，而不会计算群体结构矩阵（即默认没有群体结构）；群体结构矩阵可以是Q矩阵（fastStructure）或者PCA矩阵(admixture)等