Flanders fails to lift over X chromosome

[shoffm]

Hello.

I would like to follow up on the issue of including the X chromosome #86.

I am running Flanders with two datasets that both include the X chromosome, one that requires liftover as it is on build 37 and one that does not as it is on build 38. The dataset that does not require liftover works well, and X chromosome variants are in the resulting loci file. However, in the dataset that does require liftover the X chromosome variants do not appear in the loci or _dataset_aligned.tsv.gz files even though X chromosome variants pass the significance threshold in this dataset. It appears this is caused by the difference in how X chromosome is coded in Flanders (23, as in --chromosomes 1-23 in the input) vs the chain file (http://hgdownload.cse.ucsc.edu/goldenpath/hg19/liftOver/hg19ToHg38.over.chain.gz) where the chr is coded as X.

I think this should be easily fixed by adding a line after line 27 to s01_alpha_sort_alleles_snpid.R where the chain is imported ch <- import.chain(default_chain_file) to grep and replace "X" with "23" in the chain file. Does this sound appropriate? It would also likely be best to update chr Y to 24.

Thanks,
Sophie

[shoffm]

Hi, I just want to follow up on this as I've done a bit more digging. It appears that this update needs to be implemented in both s01_alpha_sort_alleles_snpid.R and s02_sumstat_munging_and_aligning.R to correctly lift over both the summary stats and bim file. I have updated my own copy of the hg19ToHg38_liftover function as detailed below.

s01_alpha_sort_alleles_snpid.R

hg19ToHg38_liftover <- function(
    dataset_munged,
    default_chain_file = "hg19ToHg38.over.chain"
){
  
  if(!(file.exists(default_chain_file))){
    ### Download chain file and unzip it
    exit_status = system("wget -O - http://hgdownload.cse.ucsc.edu/goldenpath/hg19/liftOver/hg19ToHg38.over.chain.gz | gunzip -c > hg19ToHg38.over.chain")

    # Raise an error if the external command fails
    if (exit_status != 0) {
      cat(paste0("Error: External command failed with exit code: ", exit_status, "\n"))
      quit(status = 1, save = "no")
    }
  }

  ch <- import.chain(default_chain_file)

  dt_for_ranges <- copy(dataset_munged)
  dt_for_ranges[, start := BP]
  dt_for_ranges <- unique(dt_for_ranges, by = c("snp_original", "CHR", "start", "BP"))
  
  dt_for_ranges[, chr_ucsc := as.character(CHR)]
  dt_for_ranges[chr_ucsc == "23", chr_ucsc := "X"]
  dt_for_ranges[chr_ucsc == "24", chr_ucsc := "Y"]
  
  dataset_ranges <- GRanges(
    seqnames = paste0("chr", dt_for_ranges$chr_ucsc),
    ranges = IRanges(start = dt_for_ranges$start, end = dt_for_ranges$BP),
    snp_original = dt_for_ranges$snp_original
  )
  
  rm(dt_for_ranges)
  gc()
  
  dataset_ranges38 <- liftOver(dataset_ranges, ch)
  dataset_ranges38_df <- as.data.table(unlist(dataset_ranges38))
  # Replace UCSC chr names back to your original numeric style
  dataset_ranges38_df[seqnames == "chrX", seqnames := "chr23"]
  dataset_ranges38_df[seqnames == "chrY", seqnames := "chr24"]
  
  setnames(dataset_ranges38_df, "end", "BP")
  dataset_ranges38_df <- dataset_ranges38_df[, .(BP, snp_original)]
  
#  dataset_lifted <- merge(dataset_munged[, !"BP"], dataset_ranges38_df, by = "snp_original", all = FALSE)

# New merge which keeps original file order - CRUCIAL TO MANTAIN ALIGNMENT WITH .BED!!
  dataset_lifted <- dataset_ranges38_df[
    dataset_munged[, !"BP"],
    on = "snp_original",
    nomatch = 0
  ]

  return(dataset_lifted)
}

s02_sumstat_munging_and_aligning.R

hg19ToHg38_liftover <- function(
    dataset_munged,
    default_chain_file = "hg19ToHg38.over.chain"
) {
  
  chain_file <- default_chain_file
  
  if (!file.exists(chain_file)) {
    message("Downloading chain file...")
    download.file("http://hgdownload.cse.ucsc.edu/goldenpath/hg19/liftOver/hg19ToHg38.over.chain.gz",
                  destfile = paste0(chain_file, ".gz"))
    R.utils::gunzip(paste0(chain_file, ".gz"))
  }
  
  ch <- import.chain(chain_file)
  
  dt_for_ranges <- copy(dataset_munged)
  dt_for_ranges[, start := BP]
  dt_for_ranges <- unique(dt_for_ranges, by = c("snp_original", "CHR", "start", "BP"))
  
  # make X and Y chr compatible with chain file for liftover
  dt_for_ranges[, chr_ucsc := as.character(CHR)]
  dt_for_ranges[chr_ucsc == "23", chr_ucsc := "X"]
  dt_for_ranges[chr_ucsc == "24", chr_ucsc := "Y"]
  
  dataset_ranges <- GRanges(
    seqnames = paste0("chr", dt_for_ranges$chr_ucsc),
    ranges = IRanges(start = dt_for_ranges$start, end = dt_for_ranges$BP),
    snp_original = dt_for_ranges$snp_original
  )
  
  rm(dt_for_ranges)
  gc()
  
  dataset_ranges38 <- liftOver(dataset_ranges, ch)
  dataset_ranges38_df <- as.data.table(unlist(dataset_ranges38))
  # Replace UCSC chr names back to your original numeric style
  dataset_ranges38_df[seqnames == "chrX", seqnames := "chr23"]
  dataset_ranges38_df[seqnames == "chrY", seqnames := "chr24"]
  
  setnames(dataset_ranges38_df, "end", "BP")
  dataset_ranges38_df <- dataset_ranges38_df[, .(BP, snp_original)]
  
  dataset_lifted <- merge(dataset_munged[, !"BP"], dataset_ranges38_df, by = "snp_original", all = FALSE)
  
  return(dataset_lifted)
}