add function seq_seq to get DNA sequence from fa

seq_slice <-	function(

				start,

				end,

				slice_size = NA,

				slice_number = NA,

				gap = 0

){

	if(!is.na(slice_number) & is.na(slice_size)){

		slice_size <- ceiling((end - start + 1 - gap * (slice_number - 1))/slice_number)

	}

	tmp_start <- seq(start,end,by = slice_size + gap)

	tmp_end <- tmp_start + slice_size - 1

	if(last(tmp_end > end)){

		tmp_end[length(tmp_end)] <- end

	}

	data.table(start = tmp_start,end = tmp_end)[,size := end - start + 1][size > 1]

}

 No newline at end of file

seq_slice <-	function(

				start,

				end,

				slice_size = NA,

				slice_number = NA,

				gap = 0

){

	if(!is.na(slice_number) & is.na(slice_size)){

		slice_size <- ceiling((end - start + 1 - gap * (slice_number - 1))/slice_number)

	}

	tmp_start <- seq(start,end,by = slice_size + gap)

	tmp_end <- tmp_start + slice_size - 1

	if(last(tmp_end > end)){

		tmp_end[length(tmp_end)] <- end

	}

	data.table(start = tmp_start,end = tmp_end)[,size := end - start + 1][size > 1]

}

seq_seq <-	function(

				organism = "mouse",

				ref = NA,

				chr,

				start = NA,

				end = NA,

				length = NA

){

	if(is.na(organism) & is.na(ref)){

		stop()

	}else if(!is.na(organism)){

	 	if(organism == "mouse"){

			ref <- "~/Reference/mm39/mm39.fa"

		}else if(organism == "human"){

			ref <- "~/Reference/hg38/hg38.fa"

		}

	}

	chr <-	chr |>

			tolower() %>% 

			gsub("chr","",.) %>% 

			paste0("chr",.)

	if(is.na(start) & is.na(end) & is.na(length)){

		stop()

	}else if(is.na(start) & ! is.na(end) & ! is.na(length)){

		start <- end - length + 1

	}else if(is.na(end) & ! is.na(start) & ! is.na(length)){

		end <- start + length - 1

	}else if(!is.na(start) & !is.na(end & is.na(length))){

		length <- end - start + 1

	}

	tmp_seq <-	paste0("samtools faidx ",ref," ",chr,":",format(start,scientific = F),"-",format(end,scientific = F)) |>

				system(intern = T)

	tmp_str <-	tmp_seq[1] |> strsplit(":")  |> unlist()

	seq <- NULL

	seq$chr <- gsub(">","",tmp_str[1])

	position <- tmp_str[2] |> strsplit("-") |> unlist()

	seq$start <- position[1]

	seq$end <- position[2]

	seq$length <- length

	seq$sequence <- paste0(tmp_seq[2:length(tmp_seq)],collapse = "")

	seq

}

 No newline at end of file