## ----setup, eval=TRUE, echo=FALSE, include=FALSE------------------------------
library("knitr") # used for a kable table later in the document

## ----eval=TRUE----------------------------------------------------------------
# Gene annotations (GTF) file. The format of the GTF file is a tab-delimited 9-column file.
# data.dir = "/inst/extdata/morex3_subset_1H_2H.gtf.gz"

gtf_file <- system.file(
  "extdata/morex3_subset_1H_2H.gtf.gz",
  package = "scShardSplitRef"
)
decompress_gtf <- gzfile(gtf_file, open = "rt")
gtf_df <- read.table(file = decompress_gtf, header = FALSE, sep = "\t")
close(decompress_gtf)
dim(gtf_df)
gtf_df[1:10, ]

## ----eval=TRUE----------------------------------------------------------------
# The format of the centromere coordinates file is a tab-delimited 3-column file, where
#   the 1st column contains chromosome/contig names,
#   the 2nd column contains the centromere coordinate (or 0 if this information is not available for your species),
#   the 3rd column contains the chromosome/contig end coordinates.
# Manually formatted centromere file (please check the expected format below):
# data.dir = "/inst/extdata/morex3_subset_1H_2H.bed"

bed_file <- system.file(
  "extdata/morex3_subset_1H_2H.bed",
  package = "scShardSplitRef"
)
centr_bed_df <- read.table(file = bed_file, header = FALSE, sep = "\t")
dim(centr_bed_df)
centr_bed_df

## ----eval=TRUE----------------------------------------------------------------
# Load the library
library(scShardSplitRef)

# Load MOREX V3 (subset to 1H, 2H)
# Note1: all files need to be decompressed for use in scShardSplitRef and Cell Ranger ARC
# Note2: use correctly formatted BED file with centromere coordinates
# Task: Check formatting of input files. Build BED-like split intervals using genomic regions (BED) and gene annotations (GTF) input files

# Read GTF file
gtf_path <- system.file(
  "extdata/morex3_subset_1H_2H.gtf.gz",
  package = "scShardSplitRef"
)
tmp_gtf_path <- tempfile(fileext = ".gtf")
con_in <- gzcon(file(gtf_path, open = "rb"))
con_out <- file(tmp_gtf_path, open = "wb")
writeBin(readBin(con_in, raw(), n = 1e8), con_out)
close(con_in)
close(con_out)

# Read BED file
bed_path <- system.file(
  "extdata/morex3_subset_1H_2H.bed",
  package = "scShardSplitRef"
)

# Define output path
output_file <- file.path(tempdir(), "split")
dir.create(output_file, recursive = TRUE, showWarnings = FALSE)
output_file_name_dir <- file.path(
  output_file,
  "morex3_subset_1H_2H_Split_Regions.bed"
)

# Determine split regions
get_split_reg <- determine_split_regions(
  bed = bed_path,
  gtf = tmp_gtf_path,
  output_bed = output_file_name_dir,
  limit = 2L^29L,
  shift_by = 20L,
  clearance = 20L
)

get_split_reg

## ----eval=TRUE----------------------------------------------------------------
bed_f <- system.file(
  "extdata/morex3_subset_1H_2H_Split_Regions.bed",
  package = "scShardSplitRef"
)
split_region_bed_df <- read.table(file = bed_f, header = FALSE, sep = "\t")
dim(split_region_bed_df)

split_region_bed_df

## ----eval=TRUE----------------------------------------------------------------
# Process GTF file and write file to tempdir()
print(output_file_name_dir)
print(tmp_gtf_path)
process_gtf(
  split_regions_bed = output_file_name_dir,
  gtf = tmp_gtf_path,
  genome_name = "subset_barley",
  genome_version = "3",
  out_path = tempdir()
)

## ----eval=TRUE----------------------------------------------------------------
gtf_modified <- list.files(
  tempdir(),
  pattern = "B1_FINAL_MODIFIED_GTF_subset.*\\.gtf$",
  full.names = TRUE
)

processed_gtf <- read.table(file = gtf_modified, header = FALSE, sep = "\t")
print(gtf_modified)
dim(processed_gtf)

sep <- as.data.frame(matrix("...", nrow = 1, ncol = 5))
colnames(sep) <- colnames(processed_gtf)[1:5]

kable(rbind(
  processed_gtf[200:203, 1:5],
  sep,
  processed_gtf[(nrow(processed_gtf) - 2):nrow(processed_gtf), 1:5]
))

## ----citation, eval=TRUE------------------------------------------------------
citation("scShardSplitRef")