Revert "Add exon-only gri filter script"

Author: jylee-bcm

utils/GRI_bed_filter/README.md

@@ -1,23 +1,13 @@
 # Genomic Region of Interest filter
-
 \
 \
 
-##### Data files for generation:
-
-1.  Gencode v46 gene annotation\
-2.  HGMD 2022 database\
-3.  Clinvar 8/6/2024 database\
-4.  SpliceAI v1.3 prediction data\
-
-##### BED file description:
-
-1.  Version 1 --- gene only: BED \<--- protein_coding gene + HGMD mutations (50 flanking bp) + Cinvar mutations (50 flanking bp)\
-    (Note: Promoter upstream 1kb is dropped as requested by Dr.Liu)\
-2.  Version 2 --- exon only: BED \<--- protein_coding gene + HGMD mutations (50 flanking bp) + Cinvar mutations (50 flanking bp) + SpliceAI predicted positions (50 flanking bp)\
-3.  Additionally, initial bed files which also contain a source column can be generated, that indicates where each position entry comes from (gencode, hgmd, clinvar, or spliceAI)
+##### Data files for generation: \
+1. Gencode v46 gene annotation \
+2. HGMD 2022 database \
+3. Clinvar 8/6/2024 database \
 
-##### BED file format:
+##### BED file description: \
+1. Version 1: BED \<--- protein_coding gene + HGMD mutations (50 flanking bp) + Cinvar mutations (50 flanking bp) \
+(Note: Promoter upstream 1kb is dropped as requested by Dr.Liu)
 
-1.  chr \| start \| end\
-2.  The 1st column contains prefix of "chr" for chromosomes, for example, "chr1, chr2, .... chr Y".

utils/GRI_bed_filter/version1_gene-only/gri_v1.R

@@ -22,7 +22,7 @@ library(readr)
 library(data.table)
 
 # Function to generate GRI region based on genome version
-generate_gri_v1 <- function(genome_version = "hg38") {
+generate_gri_geneonly <- function(genome_version = "hg38") {
   # data preparation --------------------------------------------------------
 
   ## hgmd local file name
@@ -34,7 +34,7 @@ generate_gri_v1 <- function(genome_version = "hg38") {
     hg38 <- BSgenome.Hsapiens.UCSC.hg38
     size <- seqlengths(hg38)
     chr_size <-
-      data.frame(seqnames = names(size), chr_end = as.integer(size))[c(1:24), ]
+      data.frame(seqnames = names(size), chr_end = as.integer(size))[c(1:24),]
 
     ##load gencode data
     url <-
@@ -58,7 +58,7 @@ generate_gri_v1 <- function(genome_version = "hg38") {
     hg19 <- BSgenome.Hsapiens.UCSC.hg19
     size <- seqlengths(hg19)
     chr_size <-
-      data.frame(seqnames = names(size), chr_end = as.integer(size))[c(1:24), ]
+      data.frame(seqnames = names(size), chr_end = as.integer(size))[c(1:24),]
 
     ##load gencode data
     url <-
@@ -141,15 +141,15 @@ generate_gri_v1 <- function(genome_version = "hg38") {
     mutate(source = "HGMD_2022")
 
   ## exclude "rejected variants" in HGMD database
-  hgmd_noR <- hgmd[-which(hgmd$CLASS == "R"), ]
+  hgmd_noR <- hgmd[-which(hgmd$CLASS == "R"),]
   hgmd_noR <- hgmd_noR |>
     dplyr::mutate(seqnames = paste0("chr", seqnames))
 
   ## add 50bp flanking region
   hgmd_noR <- merge(hgmd_noR, chr_size, by = "seqnames")
   hgmd_noR_flank <- hgmd_noR |>
     dplyr::mutate(
-      final_start = ifelse(start > 50, start - 50, 1),
+      final_start = ifelse(start > 50, start - 50, 0),
       final_end = ifelse((end + 50) > chr_end, chr_end, end + 50),
       .after = seqnames
     )
@@ -201,7 +201,7 @@ generate_gri_v1 <- function(genome_version = "hg38") {
   ## add flanking region
   clinvar_patho_flank <- clinvar_patho |>
     mutate(
-      final_start = ifelse(start > 50, start - 50, 1),
+      final_start = ifelse(start > 50, start - 50, 0),
       final_end = ifelse((end + 50) > chr_end, chr_end, end + 50),
       .after = seqnames
     )
@@ -233,13 +233,6 @@ generate_gri_v1 <- function(genome_version = "hg38") {
   ## check bed coverage
   coverage_bed(bed_gencode_hgmd_clinvar.dt)
 
-  # convert start position from 1-based into 0-based format of bed file
-  bed_gencode_hgmd_clinvar.dt <- bed_gencode_hgmd_clinvar.dt |>
-    mutate(start = pmax(0, start - 1))
-  
-  bed_gencode_hgmd_clinvar <- bed_gencode_hgmd_clinvar |>
-    mutate(start = pmax(0, start - 1))
-  
   # write the bed file
   file_name <- paste0("gene_only.", genome_version, ".bed")
   write_tsv(bed_gencode_hgmd_clinvar.dt[, c(1:3)],
@@ -255,5 +248,5 @@ generate_gri_v1 <- function(genome_version = "hg38") {
 }
 
 # Apply generate_gri_geneonly function to generate bed files --------------
-generate_gri_v1("hg19")
-generate_gri_v1("hg38")
+generate_gri_geneonly("hg19")
+generate_gri_geneonly("hg38")