基于MatrixEQTL进行eQTL mapping analysis——以GTEx数据库为例
MatrixEQTL是一个快速有效的cis/trans eQTL mapping分析工具，被The Genotype-Tissue Expression (GTEx)数据库用于eQTL mapping分析。使用MatrixEQTL进行eQTL mapping需要五个输入文件，即SNP.txt（snpmatrix文件，行名为样本名称，列名为snpid）， GE.txt（基因表达量文件，行名为样本名称，列名为基因或者转录本名称）， Covariates.txt（协变量文件，行名为样本名称，列名为协变量名称），geneloc.txt（基因位置文件，各列的含义分别是基因名称、基因所在染色体编号、基因起始位置和基因结束位置），snpsloc.txt（snp的物理位置文件，各列含义分别为snpid、所在染色体名称、物理位置）。
SNP.txt文件的准备
需要准备的数据：.vcf.gz或者.vcf文件。
第一步：获取重复snpid
由于第七版的GTEx基因型数据存在snpid重复的情况，为了准确计算maf，所以需要去重。
具体方法一：
① vcf文件转化为ped和map文件
plink --vcf file.vcf --recode --out file
或者
plink --vcf file.vcf --recode --out file
输出：file.ped file.map
②基于得到的file.map文件，得到重复的snpid，存储到plink.dupvar文件
cut -f 2 file.map | sort | uniq -d > plink.dupvar
具体方法二：
plink --file file.vcf --list-duplicate-vars ids-only suppress-first
输入：file.vcf
输出：plink.dupvar（重复的snpid文件）
注：方法一和方法二的区别是方法二更准确，其对于每一个重复的项保留了一项（suppress-first），而方法一删除了全部重复项。
第二步：过滤snpid和样本id
plink --vcf .vcf.gz --recode --out .name --keep .sampleid --missing --freq --maf 0.02 --exclude plink.dupvar
输入：
.vcf.gz，原始vcf文件；
.sampleid，用户想保留的sampleid；
plink.dupvar，用户想过滤掉的snpid。
输出：
.name.frq， 所有snpid的等位基因（A1，A2)以及MAF；
.name.lmiss， 所选样本的snpid确实情况；
.name.imiss， 每个样本基因型缺失情况；
.name.ped，更新后的ped文件；
.name.map，更新后的map文件。
第三步：计算snp矩阵
vcftools --gzvcf .vcf.gz --012 --out .name --keep .sampleid --maf 0.02 --exclude plink.dupvar
输入：
.vcf.gz，原始vcf文件；
.sampid，用户想保留的样品id；
plink.dupvar，用户想过滤掉的snpid。
输出：
name.012，snpmatrix（样本数*snpid数），纯数字矩阵，只有行索引（是样本id在原始.vcf.gz文件中的索引号），没有列索引。与MatrixEQTL要求的输入文件SNP.txt正好相反，所以后面需要转置；
name.012.pos，snpid 所在的染色体及物理位置；
name.012.indv，样本id。
第四步：snp矩阵转置
因为SNP.tx有格式要求，所以行名必须是snpid，列名必须是样本名称。因此需要删除name.012的第一列（即索引号），然后转置。
① 快速删除第一列Index
awk '{ $1=null;print }' name.012 > name_del0.012
②添加sampleid
paste name.012.indv name_del0.012 > name_del0.012_indv.txt
③ 快速添加行snpid
首先从.map文件中获得snpid
cat file | awk '{print $2}' > snpid.txt
echo "ID" > snpid1.txt ; cat snpid.txt >> snpid1.txt
注：snpid.txt即为去除重复后的snpid，从name.map中获得。
④合并多行为一行
可以用很多方法实现该目的，在此以python脚本sipid_2_oneline.py为例，输入文件为上述snpid1.txt，输出文件为一行snpid，以tab符分隔：
from sys import argv
script,input,output = argv
with open(output,"w") as f:
with open(input) as f1:
for line in f1.readlines():
f.write(line.strip()+"\t")
f.write("\n")
python sipid_2_oneline.py input output
⑤添加snpid并转置
cat snpid1.txt name_del0.012_indv.txt> genotype.txt
awk '{i=1;while(i <= NF){col[i]=col[i] $i " ";i=i+1}} END {i=1;while(i<=NF){print col[i];i=i+1}}' genotype.txt | sed 's/[ \t]*$//g' > genotype_transpose.txt
snp_loc.txt准备和gene_loc.txt准备
1）snp位置（snp position)
从.map中获得，只需调换位置即可。
2）基因或者转录本位置（gene position）
从ensembl网站下载。
GE.txt文件准备和Covariate.txt文件准备
根据实际情况准备GE.txt，协变量文件非必需。
运行MatrixEQTL
将上述五个（或者四个，不需要covariate.txt文件）准备好后，直接使用MatrixEQTL的示例代码进行eQTL分析。
参考资料：
R语言实现eQTL分析
