数据准备
首先要明确SMR分析需要准备哪些文件?为了方便大家去SMR官方网站查找资料时有所参照,我尽量使用官网原文,附上相应的解说:
“We store eQTL summary data in three separate files .esi (SNP information, similar as a PLINK .bim file), .epi (probe information) and .besd (a binary file to store the summary statistics).
Make a BESD file from SMR query output
注意myquery.txt
的行名和行名顺序哦,
smr --qfile myquery.txt --make-besd --out mybesd
有了besd文件后,我们继续来制备另外两种格式的文件:
Update the .esi or .epi files
ps:这里的my_beqtl
文件就是上一步生成的mybesd
文件,不加后缀名
smr --beqtl-summary my_beqtl --update-esi mybigpool.esi
--update-esi reads a .esi file for updating and backup the original .esi file.
代码语言:javascript复制smr --beqtl-summary my_beqtl --update-epi mybigpool.epi
--update-epi reads a .epi file for updating and backup the original .epi file.
要强调的是,这三个格式的文件名要求保持一致,便于后续SMR分析!
然而问题来了,myquery.txt
文件又是从何而来的呢?
靶基因的eqtl文件获取
“eQTLGen - cis-eQTLs
“输入感兴趣的基因,然后记住基因在哪条染色体上,还有Gene ID 比如,
“然后,
“这个文件1G多,下载以后解压缩,一直到文件名的尾巴分别是besd/epi/esi结尾的时候,修改文件名为:
假设你已经下载好了SMR软件,我这里就直接运行:
代码语言:javascript复制PS E: 9-MY_tryMR_PracticeSMRsmr-1.3.1-win-x86_64> .smr-1.3.1-win.exe smr --beqtl-summary ..cis-eQTL-SMR_20191212myeqtl --query 5.0e-8 --snp-chr 1
这个时候就拿到了对应染色体上p值<5.0e-8的所有基因的相关信息:
当然,也可以直接指定基因,这样就更省事了~
代码语言:javascript复制PS E: 9-MY_tryMR_PracticeSMRsmr-1.3.1-win-x86_64> .smr-1.3.1-win.exe smr --beqtl-summary ..cis-eQTL-SMR_20191212myeqtl --query 5.0e-8 --snp-chr 1 --gene ENSG00000119535 --out myquery
这里要注意,gene对应的是ENSG00000
开头的哦
打开myquery文件:
这回总算和第一步衔接上了!
SMR分析前准备
代码语言:javascript复制PS E: 9-MY_tryMR_PracticeSMRsmr-1.3.1-win-x86_64> .smr-1.3.1-win.exe smr --qfile .myquery.txt --make-besd --out myCSF3R
这个时候到对应文件夹下查看就会发现,
哇哦,齐活了!
接下来就可以开始SMR正式分析了——
SMR & HEIDI analysis
这里的结局是primary outcome,别搞错喽!另外,参考基因组记得对应到靶基因对应的染色体上哦~
准备mygwas.ma
这里有几个细节需要关注:The headers are not keywords and will be omitted by the program.
“"A1 "应为效应等位基因,"A2 "应为其他等位基因,"freq "应为 "A1 "的频率。 "n "列不会用于 SMR 或 HEIDI 分析,因此如果没有,可以用 "NA "代替。也就是把列名为n的内容赋值为NA,但不能没有n这一列!
有的GWAS数据是从别人的文献中获取的,这个时候就要手动加freq这一列,之前已经介绍过➡
SMR分析
代码语言:javascript复制PS E: 9-MY_tryMR_PracticeSMRsmr-1.3.1-win-x86_64> .smr-1.3.1-win.exe smr --bfile ..1000G.EUR.QC1000G.EUR.QC.1 --gwas-summary ..my_gwas.ma --beqtl-summary .myCSF3R --out mysmr --thread-num 10