Skip to content

Commit

Permalink
spell check
Browse files Browse the repository at this point in the history
  • Loading branch information
HendrikSchultheis committed Dec 20, 2018
1 parent d60faa7 commit 6507643
Showing 1 changed file with 38 additions and 38 deletions.
76 changes: 38 additions & 38 deletions bin/reduce_sequence.R
Original file line number Diff line number Diff line change
Expand Up @@ -3,15 +3,15 @@ library("optparse")

option_list <- list(
make_option(opt_str = c("-i", "--input"), default = NULL, help = "Input bed-file. Last column must be sequences.", metavar = "character"),
make_option(opt_str = c("-k", "--kmer"), default = 10, help = "Kmer length. Default = %default", metavar = "integer"),
make_option(opt_str = c("-k", "--kmer"), default = 10, help = "k-mer length. Default = %default", metavar = "integer"),
make_option(opt_str = c("-m", "--motif"), default = 10, help = "Estimated motif length. Default = %default", metavar = "integer"),
make_option(opt_str = c("-o", "--output"), default = "reduced.bed", help = "Output file. Default = %default", metavar = "character"),
make_option(opt_str = c("-t", "--threads"), default = 1, help = "Number of threads to use. Use 0 for all available cores. Default = %default", metavar = "integer"),
make_option(opt_str = c("-c", "--clean"), default = TRUE, help = "Delete all temporary files. Default = %default", metavar = "logical"),
make_option(opt_str = c("-s", "--min_seq_length"), default = NULL, help = "Remove sequences below this length. Defaults to the maximum value of motif and kmer and can not be lower.", metavar = "integer", type = "integer"),
make_option(opt_str = c("-n", "--minoverlap_kmer"), default = NULL, help = "Minimum required overlap between kmer. Used to create reduced sequence ranges out of merged kmer. Can not be greater than kmer length. Default = kmer - 1", metavar = "integer", type = "integer"),
make_option(opt_str = c("-v", "--minoverlap_motif"), default = NULL, help = "Minimum required overlap between motif and kmer to consider kmer significant. Used for kmer cutoff calculation. Can not be greater than motif and kmer length. Default = ceiling(motif / 2)", metavar = "integer", type = "integer"),
make_option(opt_str = c("-f", "--motif_occurence"), default = 1, help = "Define how many motifs are expected per sequence. This value is used during kmer cutoff calculation. Default = %default meaning that there should be approximately one motif per sequence.", metavar = "double")
make_option(opt_str = c("-s", "--min_seq_length"), default = NULL, help = "Remove sequences below this length. Defaults to the maximum value of motif and k-mer and can not be lower.", metavar = "integer", type = "integer"),
make_option(opt_str = c("-n", "--minoverlap_kmer"), default = NULL, help = "Minimum required overlap between k-mer. Used to create reduced sequence ranges out of merged k-mer. Can not be greater than k-mer length. Default = kmer - 1", metavar = "integer", type = "integer"),
make_option(opt_str = c("-v", "--minoverlap_motif"), default = NULL, help = "Minimum required overlap between motif and k-mer to consider k-mer significant. Used for k-mer cutoff calculation. Can not be greater than motif and k-mer length. Default = ceiling(motif / 2)", metavar = "integer", type = "integer"),
make_option(opt_str = c("-f", "--motif_occurrence"), default = 1, help = "Define how many motifs are expected per sequence. This value is used during k-mer cutoff calculation. Default = %default meaning that there should be approximately one motif per sequence.", metavar = "double")
)

opt_parser <- OptionParser(option_list = option_list,
Expand All @@ -23,19 +23,19 @@ opt <- parse_args(opt_parser)
#' Reduces each sequence to its most frequent region.
#'
#' @param input Input bed-file. Last column must be sequences.
#' @param kmer Kmer length. Default = 10
#' @param kmer k-mer length. Default = 10
#' @param motif Estimated motif length. Default = 10
#' @param output Output file. Default = reduced.bed
#' @param threads Number of threads to use. Default = 1. Use 0 for all cores.
#' @param clean Delete all temporary files.
#' @param minoverlap_kmer Minimum required overlap between kmer. Used to create reduced sequence ranges out of merged kmer. Can not be greater than kmer length . Default = kmer - 1
#' @param minoverlap_motif Minimum required overlap between motif and kmer to consider kmer significant. Used for kmer cutoff calculation. Can not be greater than motif and kmer length. Default = ceiling(motif / 2)
#' @param min_seq_length Remove sequences below this length. Defaults to the maximum value of motif and kmer and can not be lower.
#' @param motif_occurence Define how many motifs are expected per sequence. This value is used during kmer cutoff calculation. Default = 1 meaning that there should be approximately one motif per sequence.
#' @param minoverlap_kmer Minimum required overlap between k-mer. Used to create reduced sequence ranges out of merged k-mer. Can not be greater than k-mer length . Default = kmer - 1
#' @param minoverlap_motif Minimum required overlap between motif and k-mer to consider k-mer significant. Used for k-mer cutoff calculation. Can not be greater than motif and k-mer length. Default = ceiling(motif / 2)
#' @param min_seq_length Remove sequences below this length. Defaults to the maximum value of motif and k-mer and can not be lower.
#' @param motif_occurrence Define how many motifs are expected per sequence. This value is used during k-mer cutoff calculation. Default = 1 meaning that there should be approximately one motif per sequence.
#'
#' @details If there is a header supplied other then the default data.table nameing scheme ('V1', 'V2', etc.) it will be kept.
#' @details If there is a header supplied other then the default data.table naming scheme ('V1', 'V2', etc.) it will be kept.
#'
reduce_sequence <- function(input, kmer = 10, motif = 10, output = "reduced.bed", threads = NULL, clean = TRUE, minoverlap_kmer = kmer - 1, minoverlap_motif = ceiling(motif / 2), min_seq_length = max(c(motif, kmer)), motif_occurence = 1) {
reduce_sequence <- function(input, kmer = 10, motif = 10, output = "reduced.bed", threads = NULL, clean = TRUE, minoverlap_kmer = kmer - 1, minoverlap_motif = ceiling(motif / 2), min_seq_length = max(c(motif, kmer)), motif_occurrence = 1) {
if (system("which jellyfish", ignore.stdout = TRUE) != 0) {
stop("Required program jellyfish not found! Please check whether it is installed.")
}
Expand Down Expand Up @@ -76,7 +76,7 @@ reduce_sequence <- function(input, kmer = 10, motif = 10, output = "reduced.bed"

# remove sequences below minimum length
if (min_seq_length < max(c(kmer, motif))) {
stop("Minimum sequence length must be greater or equal to ", max(c(motif, kmer)), " (maximum value of kmer and motif).")
stop("Minimum sequence length must be greater or equal to ", max(c(motif, kmer)), " (maximum value of k-mer and motif).")
}

total_rows <- nrow(bed_table)
Expand All @@ -91,7 +91,7 @@ reduce_sequence <- function(input, kmer = 10, motif = 10, output = "reduced.bed"
fasta_file <- paste0(basename(input), ".fasta")
seqinr::write.fasta(sequences = as.list(bed_table[[ncol(bed_table)]]), names = bed_table[[4]], as.string = TRUE, file.out = fasta_file)

message("Counting kmer...")
message("Counting k-mer...")
# count k-mer
hashsize <- 4 ^ kmer
count_output_binary <- "mer_count_binary.jf"
Expand All @@ -105,20 +105,20 @@ reduce_sequence <- function(input, kmer = 10, motif = 10, output = "reduced.bed"

system(command = jellyfish_dump_call, wait = TRUE)

message("Reduce kmer.")
message("Reduce k-mer.")
# load mer table
# columns: kmer, count
kmer_counts <- data.table::fread(input = mer_count_table, header = FALSE)
# order kmer descending
# order k-mer descending
data.table::setorder(kmer_counts, -V2)

# compute number of hits to keep
keep_hits <- significant_kmer(bed_table, kmer = kmer, motif = motif, minoverlap = minoverlap_motif, motif_occurence = motif_occurence)
keep_hits <- significant_kmer(bed_table, kmer = kmer, motif = motif, minoverlap = minoverlap_motif, motif_occurrence = motif_occurrence)

# reduce kmer
# reduce k-mer
reduced_kmer <- reduce_kmer(kmer = kmer_counts, significant = keep_hits)

message("Find kmer in sequences.")
message("Find k-mer in sequences.")
# find k-mer in sequences
# columns: name, start, end, width
ranges_table <- find_kmer_regions(bed = bed_table, kmer_counts = reduced_kmer, minoverlap = minoverlap_kmer, threads = threads)
Expand Down Expand Up @@ -155,21 +155,21 @@ reduce_sequence <- function(input, kmer = 10, motif = 10, output = "reduced.bed"
data.table::fwrite(merged, file = output, sep = "\t", col.names = keep_col_names)
}

#' Predict how many interesting kmer are possible for the given data.
#' Predict how many interesting k-mer are possible for the given data.
#'
#' @param bed Bed table with sequences in last column
#' @param kmer Length of kmer
#' @param kmer Length of k-mer
#' @param motif Length of motif
#' @param minoverlap Minimum number of bases overlapping between kmer and motif. Must be <= motif & <= kmer. Defaults to ceiling(motif / 2).
#' @param motif_occurence Define how many motifs are expected per sequence. Default = 1
#' @param minoverlap Minimum number of bases overlapping between k-mer and motif. Must be <= motif & <= kmer. Defaults to ceiling(motif / 2).
#' @param motif_occurrence Define how many motifs are expected per sequence. Default = 1
#'
#' @return Number of interesting kmer.
significant_kmer <- function(bed, kmer, motif, minoverlap = ceiling(motif / 2), motif_occurence = 1) {
#' @return Number of interesting k-mer.
significant_kmer <- function(bed, kmer, motif, minoverlap = ceiling(motif / 2), motif_occurrence = 1) {
if (minoverlap > kmer || minoverlap > motif) {
stop("Kmer & motif must be greater or equal to minoverlap!")
}
if (motif_occurence <= 0) {
stop("Motif_occurence must be a numeric value above 0!")
if (motif_occurrence <= 0) {
stop("Motif_occurrence must be a numeric value above 0!")
}

# minimum sequence length to get all interesting overlaps
Expand All @@ -180,15 +180,15 @@ significant_kmer <- function(bed, kmer, motif, minoverlap = ceiling(motif / 2),
# reduce to max interesting length
seq_lengths <- ifelse(seq_lengths > min_seq_length, min_seq_length, seq_lengths)

# calculate max possible kmer
# calculate max possible k-mer
topx <- sum(seq_lengths - kmer + 1)

return(topx * motif_occurence)
return(topx * motif_occurrence)
}

#' Orders kmer table after count descending and keeps all kmer with a cumulative sum below the given significance threshold.
#' Orders k-mer table after count descending and keeps all k-mer with a cumulative sum below the given significance threshold.
#'
#' @param kmer Kmer data.table columns: kmer, count
#' @param kmer K-mer data.table columns: kmer, count
#' @param significant Value from significant_kmer function.
#'
#' @return reduced data.table
Expand All @@ -204,16 +204,16 @@ reduce_kmer <- function(kmer, significant) {
#' create list of significant ranges (one for each bed entry)
#'
#' @param bed Data.table of bed with sequence in last column
#' @param kmer_counts Data.table of counted kmer. Column1 = kmer, column2 = count.
#' @param minoverlap Minimum overlapping nucleotides between kmers to be merged. Positive integer. Must be smaller than kmer length.
#' @param kmer_counts Data.table of counted k-mer. Column1 = kmer, column2 = count.
#' @param minoverlap Minimum overlapping nucleotides between k-mers to be merged. Positive integer. Must be smaller than k-mer length.
#' @param threads Number of threads.
#'
#' @return Data.table with relative positions and width (start, end, width).
#'
#' TODO Include number of motifs per sequence (aka motif_occurence). Attempt to keep best 2 regions for occurence = 2? Probably high impact on performance.
#' TODO Include number of motifs per sequence (aka motif_occurrence). Attempt to keep best 2 regions for occurrence = 2? Probably high impact on performance.
find_kmer_regions <- function(bed, kmer_counts, minoverlap = 1 , threads = NULL) {
if (nchar(kmer_counts[1, 1]) <= minoverlap) {
stop("Minoverlap must be smaller than kmer length!")
stop("Minoverlap must be smaller than k-mer length!")
}

names(kmer_counts)[1:2] <- c("kmer", "count")
Expand All @@ -232,9 +232,9 @@ find_kmer_regions <- function(bed, kmer_counts, minoverlap = 1 , threads = NULL)
return(data.table::data.table(start = NA, end = NA, width = NA, name = name))
}

# add kmer sequences
# add k-mer sequences
ranges[, sub_seq := stringr::str_sub(seq, start, end)]
# add kmer count
# add k-mer count
ranges[, count := kmer_counts[ranges[["sub_seq"]], "count", on = "kmer"]]

#### reduce ranges
Expand All @@ -253,7 +253,7 @@ find_kmer_regions <- function(bed, kmer_counts, minoverlap = 1 , threads = NULL)
which(member == x)
})

# calculate component score (= sum of kmer count)
# calculate component score (= sum of k-mer count)
score <- vapply(node_membership, FUN.VALUE = numeric(1), function(x) {
sum(kmer_counts[x, "count"])
})
Expand Down

0 comments on commit 6507643

Please sign in to comment.