SNP nomenclature

I found this paper particularly useful for clarifying SNP nomenclature:

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6099125/

File formats for SNP information

Different file formats have different names for describing alleles and coding alleles.

• .bim format: they use A1 and A2 alleles, they claim A1 is usually the minor allele, see here: https://www.cog-genomics.org/plink/1.9/formats#bim

• .pvar format: they use Alt and Ref alleles, see here: https://www.cog-genomics.org/plink/2.0/formats#pvar

• .traw format: see here: https://www.cog-genomics.org/plink/2.0/formats#traw the Ref allele is counted (coded)

• Fusion summary stats files: They use A1, A2. A1 is the effect allele.

• Fusion weights files: They use .bim format for SNP information, so 5th column is Alt, 6th column is Ref. Effect size in wgt.matrix variable is defined using Alt as effect ellele (5th column in snps variable).

• SMR files: They use A1, A2. A1 is the effect allele. They also call A1 as Ref allele, A2 as Alt allele.

• ctwas format: ctwas package uses R package pgenlibr to read data. The Alt allele in .bim or .pvar files are counted, so the effect allele is Alt allele.

• plink association tests: see here. For the additive effects of SNPs, the direction of the regression coefficient represents the effect of each extra minor allele (i.e. a positive regression coefficient means that the minor allele increases risk/phenotype mean). If the –beta command is added along with –logistic, then the regression coefficients rather than the odds ratios will be returned. NOTE:Elsewhere in plink, the term reference allele is sometimes used to refer to A1, i.e. the –reference-allele command can be used to specify which allele is A1. Note that in association testing, the odds ratios, etc are typically calculated with A2 as the actual reference allele (i.e. a positive OR means A1 increases risk relative to A2).

• GTEx summary statistics: The normalized effect size (NES) of the eQTLs is defined as the slope of the linear regression, and is computed as the effect of the alternative allele (ALT) relative to the reference allele (REF) in the human genome reference GRCh38/hg38 (i.e., the eQTL effect allele is the ALT allele). See here.

• GEMMA: “chr rs ps n_miss allel1 allel0 af beta se l_remle p_wald”. The 11 columns are: chromosome, SNP ID, basepair position, number of missing values for a given SNP, the effect (coded) allele, the other allele, frequency of the effect allele, effect size, standard error, lambda and p-value ([http://www.xzlab.org/software.html]).

Consistency between GTEx v7 and UKB data

• We check if SNP information provided by FUSION GTEx v7 trained and UKB SNP information are consistent.

UKBiobank SNP information were generated by plink2, we only selected variants with MAF > 0.05 variants and missing genotype rate < 0.05.

eQTL information were obtained from FUSION pre-trained models. We used all SNPs with weights != 0 in lasso trained models.

library(ctwas)
library(plyr)
ld_pgenfn = "~/causalTWAS/ukbiobank/ukb_pgen_s40.22/ukb-s40.22.2_pgenfs.txt"
eqtlfn = "~/causalTWAS/fusion_weights/Adipose_Subcutaneous.lasso.eqtl.txt"

eqtl = read.table(eqtlfn, header = T)
colnames(eqtl)[1:6] <- c("chrom", "id", "cm", "pos", "alt", "ref")
ld_pgenfs <- read.table(ld_pgenfn, header = F, stringsAsFactors = F)[,1]
ld_pvarfs <- sapply(ld_pgenfs, prep_pvar, outputdir = "~/temp/")

i <- 1 # chrom 1 as an example
pvar <- read_pvar(ld_pvarfs[i])

cat("number of eQTLs on chromosome ", i,  "(by lasso): ", nrow(eqtl[eqtl$chrom == i,]))

number of eQTLs on chromosome  1 (by lasso):  4488

a <- join(pvar, eqtl, type = "inner", by = c("id", "pos"))
cat("number of eQTLs on chromosome ", i,  "(by lasso, shared with UKB): ", nrow(a))

number of eQTLs on chromosome  1 (by lasso, shared with UKB):  3885

b <- join(pvar, eqtl, type = "inner")

Joining by: chrom, id, pos, alt, ref

cat("number of eQTLs on chromosome ", i,  "(by lasso, shared with UKB) that match: ", nrow(b))

number of eQTLs on chromosome  1 (by lasso, shared with UKB) that match:  1117

cat("number of eQTLs on chromosome ", i,  "(by lasso, shared with UKB) that match after reverse alt and ref: ",
    which(a[!(a$id %in% b$id), 5] != a[!(a$id %in% b$id), 8]))

number of eQTLs on chromosome  1 (by lasso, shared with UKB) that match after reverse alt and ref:  2645

cat("number of eQTLs on chromosome ", i,  "(by lasso, shared with UKB) that match after flip: ",
sum(a[!(a$id %in% b$id), 5] != a[!(a$id %in% b$id), 8]))

number of eQTLs on chromosome  1 (by lasso, shared with UKB) that match after flip:  1

cat("number of eQTLs on chromosome ", i,  "(by lasso, shared with UKB) that are strand ambiguous: ",
nrow(a[(a$alt=="A" & a$ref=="T") | (a$alt=="T" & a$ref=="A") | (a$alt=="C" & a$ref=="G") | (a$alt=="G" & a$ref=="C"),]))

number of eQTLs on chromosome  1 (by lasso, shared with UKB) that are strand ambiguous:  0

Procedures to harmonize SNPs

See here about some instructions for merging two datasets: https://zzz.bwh.harvard.edu/plink/dataman.shtml#flip

It is assumed that reported two alleles are on the same strand, but it can be either positive strand or negative strand.

• FUSION filters out all strand ambiguous SNPs (i.e. AT, TA, CG, GC alleles) and flipped strand unambiguous SNPs if there is a match after flipping. It also switch alt ref SNPs if there is a match after switching.

• SMR switch alt ref SNPs if there is a match after switching. No flipping.

• ctwas currently (v.0.1.11) does not harmonize data and assumes the alleles were harmonized before running ctwas, so they do not need to be flipped or switched.

Session information

sessionInfo()

R version 3.6.1 (2019-07-05)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Scientific Linux 7.4 (Nitrogen)

Matrix products: default
BLAS/LAPACK: /software/openblas-0.2.19-el7-x86_64/lib/libopenblas_haswellp-r0.2.19.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] plyr_1.8.4   ctwas_0.1.18

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.5        pillar_1.5.1      compiler_3.6.1   
 [4] later_0.8.0       git2r_0.26.1      workflowr_1.6.2  
 [7] iterators_1.0.10  tools_3.6.1       digest_0.6.20    
[10] debugme_1.1.0     evaluate_0.14     lifecycle_1.0.0  
[13] tibble_3.1.0      lattice_0.20-38   pkgconfig_2.0.2  
[16] rlang_0.4.10      Matrix_1.2-18     foreach_1.4.4    
[19] DBI_1.1.0         yaml_2.2.0        xfun_0.8         
[22] stringr_1.4.0     dplyr_1.0.5       knitr_1.23       
[25] pgenlibr_0.2      generics_0.0.2    fs_1.3.1         
[28] vctrs_0.3.7       tidyselect_1.1.0  rprojroot_1.3-2  
[31] grid_3.6.1        glue_1.4.2        data.table_1.13.2
[34] R6_2.4.0          fansi_0.4.0       rmarkdown_1.13   
[37] purrr_0.3.4       magrittr_1.5      whisker_0.3-2    
[40] ellipsis_0.2.0.1  backports_1.1.4   promises_1.0.1   
[43] codetools_0.2-16  htmltools_0.3.6   httpuv_1.5.1     
[46] logging_0.10-108  utf8_1.1.4        stringi_1.4.3    
[49] crayon_1.3.4