#Extracts BAM reads overlapping Cyp2e1 from large BAM file
#Smaller BAM subset allows for more rapid visualization and analysis
library(Rsamtools)
library(GenomicRanges)

#Cyp2e1 in mm39 range
chromosome_to_subset <- 7
range_start_to_subset <- 140341652
range_end_to_subset <- 140356894

bam_directory_path <- "D:/BAM/"
bam_output_directory_path <- "D:/BAM_subsets/"
prefix_name <- "Cyp2e1"
bam_paths <- list.files(bam_directory_path, pattern = ".bam$")

output_range <- GRanges(
  seqnames = chromosome_to_subset,
  ranges = IRanges(start = range_start_to_subset, end = range_end_to_subset),
  strand = c("+", "-")
)

for(bam_num in 1:length(bam_paths)){ #bam_num <- 1

filterBam(file = paste0(bam_directory_path, bam_paths[bam_num])
          ,destination = paste0(bam_output_directory_path,prefix_name,bam_paths[bam_num])
          ,param = ScanBamParam(which = output_range))
}