## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>"
)

## -----------------------------------------------------------------------------
library(SpliceWiz)

## ----eval = FALSE-------------------------------------------------------------
#  ref_path <- "./Reference"

## ----eval=FALSE---------------------------------------------------------------
#  buildRef(
#      reference_path = ref_path,
#      fasta = "genome.fa", gtf = "transcripts.gtf",
#      genome_type = "hg38"
#  )

## ----eval=FALSE---------------------------------------------------------------
#  getResources(
#      reference_path = ref_path,
#      fasta = "genome.fa",
#      gtf = "transcripts.gtf"
#  )
#  
#  buildRef(
#      reference_path = ref_path,
#      fasta = "", gtf = "",
#      genome_type = "hg38"
#  )

## ----eval=FALSE---------------------------------------------------------------
#  buildRef(
#      reference_path = ref_path,
#      fasta = "genome.fa",
#      gtf = "transcripts.gtf",
#      genome_type = "hg38"
#  )

## ----eval=FALSE---------------------------------------------------------------
#  # Assuming hg38 genome:
#  
#  buildRef(
#      reference_path = ref_path,
#      genome_type = "hg38",
#      overwrite = TRUE
#  )

## ----eval=FALSE---------------------------------------------------------------
#  FTP <- "ftp://ftp.ensembl.org/pub/release-94/"
#  
#  buildRef(
#      reference_path = ref_path,
#      fasta = paste0(FTP, "fasta/homo_sapiens/dna/",
#          "Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz"),
#      gtf = paste0(FTP, "gtf/homo_sapiens/",
#          "Homo_sapiens.GRCh38.94.chr.gtf.gz"),
#      genome_type = "hg38"
#  )

## -----------------------------------------------------------------------------
require(AnnotationHub)

ah <- AnnotationHub()
query(ah, "Ensembl")

## -----------------------------------------------------------------------------
query(ah, c("Homo Sapiens", "release-94"))

## ----eval=FALSE---------------------------------------------------------------
#  buildRef(
#      reference_path = ref_path,
#      fasta = "AH65745",
#      gtf = "AH64631",
#      genome_type = "hg38"
#  )

## ----eval=FALSE---------------------------------------------------------------
#  buildRef(
#      reference_path = ref_path,
#      fasta = "genome.fa", gtf = "transcripts.gtf",
#      genome_type = ""
#  )

## -----------------------------------------------------------------------------
getAvailableGO()

## ----eval=FALSE---------------------------------------------------------------
#  buildRef(
#      reference_path = ref_path,
#      fasta = "genome.fa", gtf = "transcripts.gtf",
#      genome_type = "",
#      ontologySpecies = "Arabidopsis thaliana"
#  )

## -----------------------------------------------------------------------------
STAR_version()

## ----eval = FALSE-------------------------------------------------------------
#  ref_path = "./Reference"
#  
#  # Ensure genome resources are prepared from genome FASTA and GTF file:
#  
#  if(!dir.exists(file.path(ref_path, "resource"))) {
#      getResources(
#          reference_path = ref_path,
#          fasta = "genome.fa",
#          gtf = "transcripts.gtf"
#      )
#  }
#  
#  # Generate a STAR genome reference:
#  STAR_BuildRef(
#      reference_path = ref_path,
#      n_threads = 8
#  )
#  

## ----eval = FALSE-------------------------------------------------------------
#  STAR_BuildRef(
#      reference_path = ref_path,
#      STAR_ref_path = "/path/to/another/directory",
#      n_threads = 8
#  )

## ----eval = FALSE-------------------------------------------------------------
#  # Generate a STAR genome reference:
#  STAR_buildGenome(
#      reference_path = ref_path,
#      STAR_ref_path = "/path/to/hg38"
#      n_threads = 8
#  )

## ----eval = FALSE-------------------------------------------------------------
#  STAR_new_ref <- STAR_loadGenomeGTF(
#      reference_path = ref_path,
#      STAR_ref_path = "/path/to/hg38",
#      STARgenome_output = file.path(tempdir(), "STAR"),
#      n_threads = 8,
#      sjdbOverhang = 100,
#      extraFASTA = "./ercc.fasta"
#  )

## ----eval = FALSE-------------------------------------------------------------
#  STAR_mappability(
#    reference_path = ref_path,
#    STAR_ref_path = file.path(ref_path, "STAR"),
#    map_depth_threshold = 4,
#    n_threads = 8,
#    read_len = 70,
#    read_stride = 10,
#    error_pos = 35
#  )

## ----eval=FALSE---------------------------------------------------------------
#  buildFullRef(
#      reference_path = ref_path,
#      fasta = "genome.fa", gtf = "transcripts.gtf",
#      genome_type = "",
#      use_STAR_mappability = TRUE,
#      n_threads = 8
#  )

## ----eval = FALSE-------------------------------------------------------------
#  require(Rsubread)
#  
#  # (1a) Creates genome resource files
#  
#  ref_path <- file.path(tempdir(), "Reference")
#  
#  getResources(
#      reference_path = ref_path,
#      fasta = chrZ_genome(),
#      gtf = chrZ_gtf()
#  )
#  
#  # (1b) Systematically generate reads based on the SpliceWiz example genome:
#  
#  generateSyntheticReads(
#      reference_path = ref_path
#  )
#  
#  # (2) Align the generated reads using Rsubread:
#  
#  # (2a) Build the Rsubread genome index:
#  
#  subreadIndexPath <- file.path(ref_path, "Rsubread")
#  if(!dir.exists(subreadIndexPath)) dir.create(subreadIndexPath)
#  Rsubread::buildindex(
#      basename = file.path(subreadIndexPath, "reference_index"),
#      reference = chrZ_genome()
#  )
#  
#  # (2b) Align the synthetic reads using Rsubread::subjunc()
#  
#  Rsubread::subjunc(
#      index = file.path(subreadIndexPath, "reference_index"),
#      readfile1 = file.path(ref_path, "Mappability", "Reads.fa"),
#      output_file = file.path(ref_path, "Mappability", "AlignedReads.bam"),
#      useAnnotation = TRUE,
#      annot.ext = chrZ_gtf(),
#      isGTF = TRUE
#  )
#  
#  # (3) Analyse the aligned reads in the BAM file for low-mappability regions:
#  
#  calculateMappability(
#      reference_path = ref_path,
#      aligned_bam = file.path(ref_path, "Mappability", "AlignedReads.bam")
#  )
#  
#  # (4) Build the SpliceWiz reference using the calculated Mappability Exclusions
#  
#  buildRef(ref_path)
#  

## ----eval = FALSE-------------------------------------------------------------
#  buildRef(ref_path, genome_type = "hg38")

## -----------------------------------------------------------------------------
STAR_version()

## ----eval = FALSE-------------------------------------------------------------
#  STAR_alignReads(
#      fastq_1 = "sample1_1.fastq", fastq_2 = "sample1_2.fastq",
#      STAR_ref_path = file.path(ref_path, "STAR"),
#      BAM_output_path = "./bams/sample1",
#      n_threads = 8,
#      trim_adaptor = "AGATCGGAAG"
#  )

## ----eval = FALSE-------------------------------------------------------------
#  Experiment <- data.frame(
#      sample = c("sample_A", "sample_B"),
#      forward = file.path("raw_data", c("sample_A", "sample_B"),
#          c("sample_A_1.fastq", "sample_B_1.fastq")),
#      reverse = file.path("raw_data", c("sample_A", "sample_B"),
#          c("sample_A_2.fastq", "sample_B_2.fastq"))
#  )
#  
#  STAR_alignExperiment(
#      Experiment = Experiment,
#      STAR_ref_path = file.path("Reference_FTP", "STAR"),
#      BAM_output_path = "./bams",
#      n_threads = 8,
#      two_pass = FALSE
#  )

## ----eval = FALSE-------------------------------------------------------------
#  # Assuming sequencing files are named by their respective sample names
#  fastq_files <- findFASTQ(
#      sample_path = "./sequencing_files",
#      paired = TRUE,
#      fastq_suffix = ".fq.gz", level = 0
#  )

## ----eval = FALSE-------------------------------------------------------------
#  STAR_alignExperiment(
#      Experiment = fastq_files,
#      STAR_ref_path = file.path("Reference_FTP", "STAR"),
#      BAM_output_path = "./bams",
#      n_threads = 8,
#      two_pass = FALSE
#  )

## ----eval=FALSE---------------------------------------------------------------
#  bams <- findBAMS("./bams", level = 1)

## ----eval=FALSE---------------------------------------------------------------
#  # assume SpliceWiz reference has been generated in `ref_path` using the
#  # `buildRef()` function.
#  
#  processBAM(
#      bamfiles = bams$path,
#      sample_names = bams$sample,
#      reference_path = ref_path,
#      output_path = "./pb_output",
#      n_threads = 4,
#      useOpenMP = TRUE
#  )

## ----eval=FALSE---------------------------------------------------------------
#  BAM2COV(
#      bamfiles = bams$path,
#      sample_names = bams$sample,
#      output_path = "./cov_output",
#      n_threads = 4,
#      useOpenMP = TRUE
#  )

## -----------------------------------------------------------------------------
se <- SpliceWiz_example_NxtSE()

cov_file <- covfile(se)[1]

cov_negstrand <- getCoverage(cov_file, strand = "-")
bw_file <- file.path(tempdir(), "sample_negstrand.bw")
rtracklayer::export(cov_negstrand, bw_file, "bw")

## ----eval=FALSE---------------------------------------------------------------
#  expr <- findSpliceWizOutput("./pb_output")

## ----eval = FALSE-------------------------------------------------------------
#  collateData(
#      Experiment = expr,
#      reference_path = ref_path,
#      output_path = "./NxtSE_output"
#  )

## ----eval = FALSE-------------------------------------------------------------
#  collateData(
#      Experiment = expr,
#      reference_path = ref_path,
#      output_path = "./NxtSE_output_novelSplicing",
#      novelSplicing = TRUE
#  )

## ----eval = FALSE-------------------------------------------------------------
#  se <- makeSE("./NxtSE_output")

## -----------------------------------------------------------------------------
sessionInfo()