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.

• 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.

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.

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

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.12

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.5        compiler_3.6.1    pillar_1.4.2     
 [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] evaluate_0.14     tibble_2.1.3      lattice_0.20-38  
[13] pkgconfig_2.0.2   rlang_0.4.10      Matrix_1.2-18    
[16] foreach_1.4.4     yaml_2.2.0        xfun_0.8         
[19] stringr_1.4.0     dplyr_0.8.3       knitr_1.23       
[22] vctrs_0.3.1       fs_1.3.1          rprojroot_1.3-2  
[25] grid_3.6.1        tidyselect_1.1.0  glue_1.3.1       
[28] data.table_1.13.2 R6_2.4.0          rmarkdown_1.13   
[31] purrr_0.3.4       magrittr_1.5      backports_1.1.4  
[34] promises_1.0.1    codetools_0.2-16  htmltools_0.3.6  
[37] assertthat_0.2.1  httpuv_1.5.1      logging_0.10-108 
[40] stringi_1.4.3     crayon_1.3.4