SMR 教程

1.数据获取

1.1eqtl数据

qtl数据获取主要从GTEx数据库，可以得到各种组织的数据。GTEx数据库的速度很慢建议从SMR官网的数据库下载（华人做的）SMR | Yang Lab (westlake.edu.cn)里面有各种qtl的数据集，下得很快！

eqtl还有血液的数据，样本量很大eQTLGen - cis-eQTLs

直接下载besd二进制格式即可，后期数据可以用smr软件自提

1.2GWAS数据

GWAS数据获取的方法很多可以直接去IEU OpenGWAS project (mrcieu.ac.uk)下载，也可以通过R包TwoSampleMR直接从网站上拽。available_outcomes()函数可以看到openGWAS里面所有的gwas数据

当然还有很多别的获取手段。也可以从finngen的GWAS获取

1.3参考基因组

用R包下载Download 1000G — download_1000G • bigsnpr (privefl.github.io)

也可以直接链接：http://fileserve.mrcieu.ac.uk/ld/1kg.v3.tgz

也可以自行从官网下载

不难发现，所需数据全是公共数据，而且获取比较轻松。

2.SMR

2.1软件下载

SMR | Yang Lab (westlake.edu.cn)下载即可，Windows和Linux均可

1	smr --bld mybld --gwas-summary mygwas.ma --beqtl-summary myeqtl --out mysmr --thread-num 10

smr输入有

--bld 指的是1000人参考基因组的目录，

--gwas-summary输入的是GWAS数据，

--beqtl-summary 输入的是eqtl数据，

--query 5e-8

--gene ENSG00000113161

--out 输出文集名

--make-besd

*mygwas.ma*

SNP    A1  A2  freq    b   se  p   n
rs1001    A   G   0.8493  0.0024  0.0055  0.6653  129850
rs1002    C   G   0.03606 0.0034  0.0115  0.7659  129799
rs1003    A   C   0.5128  0.045   0.038   0.2319  129830

exp<-fread('GCST004744-EFO_0000571-Build37.f.tsv.gz',data.table = F,header=T)
tail(exp)
exp$n<-66756
##ebi常用格式
LUAD<-exp%>%
  dplyr::filter(str_detect(variant_id, "^rs"))%>%
  dplyr::rename(
    SNP=variant_id,
    A1=effect_allele, 
    A2=other_allele,
    freq=effect_allele_frequency, 
    b=beta,
    se=standard_error,
    p=p_value,
    N=n 
  )%>%
  dplyr::select(SNP, A1, A2, freq, b, se, p, N) %>%
  dplyr::filter(!str_detect(SNP, ","))
##finngen常用格式
LUAD<-exp%>%
  dplyr::filter(str_detect(rsids, "^rs"))%>%
  dplyr::rename(
    SNP=rsids,
    A1=alt, 
    A2=ref,
    freq=af_alt, 
    b=beta,
    se=sebeta,
    p=pval,
    N=n 
  )%>%
  dplyr::select(SNP, A1, A2, freq, b, se, p, N) %>%
  dplyr::filter(!str_detect(SNP, ","))
head(LUAD,804)
write.table(LUAD,'LUAD.ma', sep = "\t", row.names = FALSE, col.names = TRUE, quote = FALSE)

SMR 具体eqtlgen转换要求

#超算平台
cd "/public/home/fanfangzhou/apps/SMR/smr_linux" 

##这里获取的是含有eQTLGen的HMGCR这基因的besd格式的
"/public/home/fanfangzhou/apps/SMR/smr_linux" 
	--beqtl-summary "/public/home/fanfangzhou/R/GTEx8/SMR/eQTLGen/cis-eQTLs-full_eQTLGen_AF_incl_nr_formatted_20191212.new.txt_besd-dense" 
	--query 1 
	--gene ENSG00000113161 
	--out "/public/home/fanfangzhou/R/2024science/6. Drugs mendelian/rotator cuff syndrome/rotator cuff injury/SMR/" 
	--make-besd

##这里获取了HMGCR的besd格式文件，就可以在这个区间进行SMR分析
"/public/home/fanfangzhou/apps/SMR/smr_linux"  
	--bfile "/public/home/fanfangzhou/R/Rcodes/486metabolites/1kg.v3/EUR" 
	--gwas-summary SMR_rotatorcuffR9.ma 
	--beqtl-summary HMGCR 
	--out rotatorblood  
	--thread-num 10

2.3运行smr

命令很简单，找着官网的改一改就行

1
2

##超算平台
"/public/home/fanfangzhou/apps/SMR/smr_linux"  --bfile "/public/home/fanfangzhou/R/Rcodes/486metabolites/1kg.v3/EUR" --gwas-summary LUAD.ma --beqtl-summary "/public/home/fanfangzhou/R/GTEx8/SMR/Whole_Blood/Whole_Blood" --out Lung_SMR_wholeblood --thread-num 10 --diff-freq-prop 0.92

–diff-freq-prop 参数是等位基因频率的设定阈值，如果小的话程序是跑不通的。因为我输入的eqtl和gwas数据没经过过滤，被筛掉的snp的比例非常大，所以需要把这个阈值设高点防止程序终止。

3.共定位

细心的朋友会发现，刚才我们并没有下载eqtl的txt列表格式的文件，理论上和gwas一样，eqtl也应该有一个非常大数据量的列表文件，那么如何得到呢，依然是靠smr软件

这里可以用eqtlgen的也可以用V8的，***–make-besd加这个就是输出besd格式的（三种文件那种），或者不加就是默认*txt

1	smr-1.3.1-win.exe --beqtl-summary Artery_Aorta --query 0.05 --out Aorta_eqtl

smr可以从二进制的besd文件中导出txt格式的eqtl文件，得到这里是Aorta_eqtl.txt文件

4.SMR locus plot

SMR | Yang Lab (westlake.edu.cn)

SMR locus plot

Here we provide an R script to plot SMR results as presented in Zhu et al. (2016 Nature Genetics). The data file for plot can be generated by the command below. The R script is avaliable at Download.

# SMR command line to generate a data file for plot

1	smr --bfile mydata --gwas-summary mygwas.ma --beqtl-summary myeqtl --out myplot --plot --probe ILMN_123 --probe-wind 500 --gene-list glist-hg19

5.Multi-SNP-based SMR test

Below shows an option to combine the information from all the SNPs in a region that pass a p-value threshold (the default value is 5.0e-8 which can be modified by the flag –peqtl-smr) to conduct a multi-SNP-based SMR analysis (Wu et al. 2018 Nature Communications).

The SNPs are pruned for LD using a weighted vertex coverage algorithm with a LD r2 threshold (the default value is 0.9 which can be modified by the flag –ld-pruning) and eQTL p-value as the weight.

1	smr --bfile mydata --gwas-summary mygwas.ma --beqtl-summary myeqtl --out mymulti --smr-multi

–smr-multi turns on set-based SMR test in the cis-region.

1	smr --bfile mydata --gwas-summary mygwas.ma --beqtl-summary myeqtl --out mymulti --smr-multi --set-wind 500

–set-wind defines a window width (Kb) centred around the top associated cis-eQTL to select SNPs in the cis-region. The default value is -9 resulting in selecting SNPs in the whole cis-region if this option is not specified.

1	smr --bfile mydata --gwas-summary mygwas.ma --beqtl-summary myeqtl --out mymulti --smr-multi --ld-multi-snp 0.1

–ld-multi-snp LD r-squared threshold used to prune SNPs (eQTLs) in the Multi-SNP based SMR test. The default value is 0.1.

最初提出的定义独立顺式-eQTL的阈值为r2< 0.9。然而，该阈值可能包括具有高 LD 的 SNP（定义为 r 2≥ 0.1）。因此，我们采用了 r 2< 0.1 作为定义独立 SNP 的阈值

参考文献：基于转录组全汇总数据的孟德尔随机化分析揭示了 38 个与严重 COVID-19 相关的新基因 - PMC (nih.gov)

GTEx V8 基因注释

GTEx Portal网址里面的：https://storage.googleapis.com/adult-gtex/bulk-gex/v8/rna-seq/GTEx_Analysis_2017-06-05_v8_RNASeQCv1.1.9_gene_reads.gct.gz

GTEx<-read.table("GTEx_Analysis_2017-06-05_v8_RNASeQCv1.1.9_gene_reads.gct", skip = 2, header = TRUE, sep = "\t") save(GTEx,file = 'GTEx.Rdata') 
GTEx[1:4,1:4] ##行是基因 列是样本 
# Name    Description       GTEX.1117F.0226.SM.5GZZ7 GTEX.1117F.0426.SM.5EGHI 
# 1 ENSG00000223972.5 DDX11L1 0 0 
# 2 ENSG00000227232.5 WASH7P 187 109 
# 3 ENSG00000278267.1 MIR6859-1 0 0 
# 4 ENSG00000243485.5 MIR1302-2HG 1 0 
colnames(GTEx)

5.常用数据

V8 release of the GTEx eQTL/sQTL summary data

V8 release of the GTEx eQTL/sQTL summary data (n = 73 - 670; [GTEx Consortium 2020 Science](https://www.science.org/doi/10.1126/science.aaz1776?url_ver=Z39.88-2003&rfr_id=ori:rid:crossref.org&rfr_dat=cr_pub 0pubmed)) in SMR binary (BESD) format:

GTEx_V8_cis_eqtl_summary(hg19) (48 GB)

GTEx_V8_cis_eqtl_summary_lite (hg19) Lite version of the GTEx V8 eqtl data (only SNPs with P < 1e-5 are included; 506MB)

This is a set of cis-eQTL summary data across 49 human tissues from the GTEx project. Only SNPs within 1Mb of the transcription start site are available. The forth column of the *.epi file is the middle position of the probe sequence rather than the transcription start site.

cis-eQTLs

This page contains the cis-eQTL results. The statistically significant cis-eQTLs and SMR-prioritised genes for several traits are browsable, the other files can be downloaded.

Downloads

Significant cis-eQTLs
Full cis-eQTL summary statistics
README
Initial release Significant cis-eQTLs (deprecated)
Initial release Full cis-eQTL summary statistics (deprecated)

6.SMR使用技巧

6.1.批量使用49个组织的SMR

首先在R里面创建我们需要的sh格式的文件，当然我们可以进行txt简单创建

original_string1 <- '
echo Whole_Blood
"/public/home/fanfangzhou/apps/SMR/smr_linux" --bfile "/public/home/fanfangzhou/R/Rcodes/486metabolites/1kg.v3/EUR" --gwas-summary LUAD.ma --beqtl-summary "/public/home/fanfangzhou/R/GTEx8/SMR/Whole_Blood/Whole_Blood" --out Lung_SMR_Whole_Blood --thread-num 10 --diff-freq-prop 0.92
'
original_string2 <- '
echo Whole_Blood
"/public/home/fanfangzhou/apps/SMR/smr_linux" --bfile "/public/home/fanfangzhou/R/Rcodes/486metabolites/1kg.v3/EUR" --gwas-summary LUAD.ma --beqtl-summary "/public/home/fanfangzhou/R/GTEx8/SMR/Whole_Blood/Whole_Blood" --out Lung_SMR_Whole_Blood_test  --smr-multi  --ld-multi-snp 0.1 
'


# 待替换的字符串
replace_values <- original_string1  # 替换值的向量

# 用替换值逐个替换原始字符串中的"Whole_Blood"
replaced_strings <- sapply(replace_values, function(value) {
  gsub("Whole_Blood", value, original_string1)
})

writeLines(replaced_strings, "script.sh")

command <- '
#!/bin/bash

ls

# 输出 "Hello, world!" 消息
echo "Hello, world!"
'
writeLines(command, "script.sh")

然后我们在linux系统里面使用sh格式运行

#使用 chmod 命令将该文件设置为可执行权限：
chmod +x my_script.sh

#执行该脚本文件：
./my_script.sh

6.2.SMR 执行mult

"/public/home/fanfangzhou/apps/SMR/smr_linux" --bfile "/public/home/fanfangzhou/R/Rcodes/486metabolites/1kg.v3/EUR" --gwas-summary ON.ma --beqtl-summary "/public/home/fanfangzhou/R/GTEx8/SMR/eQTLGen/cis-eQTLs-full_eQTLGen_AF_incl_nr_formatted_20191212.new.txt_besd-dense"  --out osteonecrosis_blood  --smr-multi  --ld-multi-snp 0.1