positives_tmp <- rbind(selected_fragments$positive, selected_fragments_identifier$positive)
rejected_tmp <- rbind(selected_fragments$rejected, selected_fragments_identifier$rejected)
selected_fragments <- list("positive"=positives_tmp, "rejected"=rejected_tmp)
}
}
selected_fragments$positive
print("Normalized and Selected LS-GKM output...")
print("Printing results to files...")
create_output(selected_fragments,
paste(path_to_output_result, paste("BOriS_", strsplit(strsplit(path_to_data, "/")[[1]]
[length(strsplit(path_to_data, "/")[[1]])], "\\.fa")[[1]][1], "_fa.fasta", sep=""), sep = "/"),
paste(path_to_output_result, paste("BOriS_", strsplit(strsplit(path_to_data, "/")[[1]]
[length(strsplit(path_to_data, "/")[[1]])], "\\.fa")[[1]][1], "_df.csv", sep=""), sep = "/"),
paste(path_to_output_result, paste("BOriS_", strsplit(strsplit(path_to_data, "/")[[1]]
[length(strsplit(path_to_data, "/")[[1]])], "\\.fa")[[1]][1], "_ambiguous_df.csv", sep=""), sep = "/")
)
call_BOriS_with_ftp_file <- function(ftpurl, targetfolder, source, target){
library(RCurl)
library(R.utils)
filenames = strsplit(getURL(ftpurl, dirlistonly = TRUE), "\n")[[1]]
correctname <- filenames[grep("_genomic.fna.gz", filenames)][-grep("_from_", filenames[grep("_genomic.fna.gz", filenames)])]
download.file(paste(ftpurl, correctname, sep = ""), paste(targetfolder, correctname, sep = ""), quiet=T)
gunzip(paste(targetfolder, correctname, sep = ""))
newfilename = sub(".gz", "", paste(targetfolder, correctname, sep = ""))
call_BOriS(path_to_data = newfilename,
path_to_seeds = paste(source, "seeds.fa", sep = ""),
path_to_lsgkmpredict = "gkmpredict",
path_to_model = paste(source, "ls-gkm_train_model_1250_L10_K6_D4_T4.model.txt", sep = ""),
path_to_output_result = target,
calibratR_model_path = paste(source, "ls-calibratR_model_1250_train_selected.RDS", sep = ""),
1250, FALSE, single_chromosome = TRUE)
}
# R-Function ----------------------------
#
# Glossar:
#
# Seed - Motive mostly found in ori-areas, basepairs of length 9
# Seedsequencelist - List of seeds found in data of C.weigle
# Fragment - part of the genome with length 1250 and seedmotive in the middle. ( XXXXXXXXXXX_Seed_XXXXXXXXXXX )
#            The area around the seed most likely showed to be the place the ori can be located.
#
# Input:
# The FASTA-File input in either a file fillled with one or more sequences.
# These sequences are genoms either of the type linear or circular.
#
# Output:
# The output will be a FASTA-File with oriC-sequences, their positions and the probability of beeing an ori
#
# Example:
# >Name_of_sequence    position_begin_seed:posisiton_end_seed    probability
#  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAATTTTTTTTTTTTTTTTTTTTTGGGGGGGGGGGGGGGGGGGGGGGGCCCCCCCCCCCCCCCCC
#
# >NC_0009.3.3    315000:316250    80%
#  ATCGGATCAAAGTAGCAAATACGGGACCCGTAAAAAAAAAA....
#
#
#
# Steps:
# 1. Input Check -
#    First the Input-File will be checked if the requirements are fullfilled - FASTA-File & unique ID's
# 2. Preparing -
#    For each sequence "s" in the FASTA-File the following steps are preformed.
#       2.1. Every sequence s will be searched though for seedsequencesmotives.
#       2.2. For each located seed in the sequence s there will be made a fragment "f" of length 1250.
#       2.3. linear / circular
# 3. Each fragment f will be predicted by the LS-GKM Model (t=4, l=10, k=6, d=4) by Dongwon Lee (https://github.com/Dongwon-Lee/lsgkm)
# 4. The predictions for the fragments will be cuttoff or calibrated
# 5. All rated fragments f of all sequences s will be returned as a FASTA-File
#
# Input Check ===========
# Preparing ==============
# Make sure all Sequences have individuall ID's
rename_FASTA_ID<- function(data){
names <- names(data)
counter<-1
for(i in names){
i <-paste(counter,"/begin:end/percent/Description:",i,sep="")#paste(paste(data[[counter]]))#
names[counter]<-i
counter<-counter+1
}
names(data)<-names
return(data)
}
# Predict ========
# Input:
# FASTA-File with genome data & seeds to be searched for in the data
# Output:
# Location of the searched seeds in the data.
#
# Example output:
#
# Dataframe of Seeds found in sequences, with:
#                                             seedpattern
#                                             sequencename = name of FASTA-Input-File (genomedata),
#                                             begin = position of the seed found in sequence in,
#                                             end = position of the seed found in sequence end,
#                                             sequencelength = length of the sequence of FASTA-Input-File (genomedata),
#                                             fragment = the fragment around the seed (here still empty, field is a placeholder),
#                                             non_calibrated = the endresult of prediction for the fragment (here still empty, field is a placeholder),
#                                             calibrated = the calibrated endresult of prediction for the fragment (here still empty, field is a placeholder)
#
#Example:
#
#        seed                                       sequencename   begin     end sequencelength fragment non_calibrated calibrated
#   TGTTAATAA  74/196:204/percent/Description:271_negative_400 0     196     204            400    FALSE          FALSE      FALSE
#   TGTTAATAA 158/139:147/percent/Description:191_negative_400 0     139     147            400    FALSE          FALSE      FALSE
#   TGTTAATAA 274/163:171/percent/Description:285_negative_400 0     163     171            400    FALSE          FALSE      FALSE
#   TGTTAATAA   315/196:204/percent/Description:4_negative_400 0     196     204            400    FALSE          FALSE      FALSE
#
#
search_for_seeds<- function(data, seeds, fragmentlength){
X <- c()
for(s in 1:length(seeds)){
# for each seed in our seedlist the given sequences will be looked through
pattern <- names(seeds)[s]
genome_as_string <- toString(data)
sequencelength <- width(data)
# if the seed is part of the pattern
if(str_detect(genome_as_string,pattern)){
# we locate him and save the positions in the result vector
# as well as we're writing the position into the name og the sequence
matching <- str_locate_all(genome_as_string,pattern)
matching_matrix<- as.matrix(as.matrix(matching)[[1]])
for(i in 1:length(matching_matrix[,1])){
begin <-matching_matrix[[i,1]] - ((fragmentlength-9)/2+0.5)
if (begin < 0){
begin <- sequencelength + begin
}
end <-matching_matrix[[i,2]] + ((fragmentlength-9)/2+0.5)
if (end > sequencelength){
end <- end - sequencelength
}
seqname <- strsplit(strsplit(names(data)[1], "Description:")[[1]][2], " ")[[1]][1]
X<- rbind(X, c(pattern, seqname, begin, end, sequencelength))
}
}
}
# changing the list into a data.frame because R is stupid
seed <- as.vector(X[,1])
sequencename <- as.vector(X[,2])
begin <- as.vector(X[,3])
end <-as.vector(X[,4])
sequencelength <- as.vector(X[,5])
identifier <- as.vector(X[,2])
fragment <- vector(mode="logical",length=length(seed))
raw_output <- vector(mode="logical",length=length(seed))
seed_dataframe <- data.frame(seed, identifier, begin, end, sequencelength, fragment, raw_output)
return(seed_dataframe)
}
# Input:
# FASTA-file with genome data & seeds to be searched for in the data
# Will create a seed_dataframe_test generated by the search_for_seeds function
# If the genomedata in the FASTA-file is linear = TRUE (default) otherwise if it is circular, then linear = FALSE
#
# Output:
# fragmentdataframe - List of fragments of length n = fragmentsize  (default = 1250)
#
# 1-621 XXXXXXXXX 630-1250
#
#
get_fragments_of_size_n <- function(data, seeds, path_to_fragmentfile, n = 1250, linear = TRUE){
genome_sequence <- toString(data)
genome_length <- width(data)
# call search_for_seeds to generate the seed_dataframe
fragment_dataframe <- search_for_seeds(data, seeds, n)
#for each row in the data.frame (looks like this)
#          seed                                       sequencename   begin     end sequencelength fragment non_calibrated calibrated
# 1   TGTTAATAA  74/196:204/percent/Description:271_negative_400 0     196     204            400    FALSE          FALSE      FALSE
#
# we will calculate the fragment of size fragmentsize  around the seed.
# to center the seed we need a fixed number of entries left and right of the seed
#
# find out the seedlength, to calculate the left- and rightside dependent on the seed
seedlength <- length(seeds[[1]])
# for each seed in the dataframe we'll generate a fragment
for(i in 1:nrow(fragment_dataframe)){
fragment_dataframe_row <- as.matrix(as.matrix(fragment_dataframe)[i,])
begin <-  as.numeric(fragment_dataframe_row[3])
end <-  as.numeric(fragment_dataframe_row[4])
# CASE 1: GENOME IS LINEAR
if(linear){
# linear = TRUE; (the geonme can be interpretated as a line)
# (Position means position on the genome)
#
#
# CASE_1.1: Seed is in the middle of the linear genome
# IF the starting position of the seed minus the needed leftside of the fragment is not surpassing the beginnboarder of the geonme (here pos = 0)
# &&
# IF the end point of the seed plus the rightside of the fragment is not supassing the endboarder of the geonme the searched fragment is in the middle
# and can be cutted out. If not so, the fragment can't be created because the genome isn't circular!
#
if((begin >= 0) && (end <= genome_length)){
fragment <-substr(genome_sequence, begin, end)
#Check if fragment contains other characters than A,C,G,T; if yes -> don't accept
if(!(grepl("[^ACGT]",fragment))){fragment_dataframe[i,6]<- fragment}
}
# CASE 2: GENOME IS CIRCULAR
}else{
# linear = FALSE; the geonme can be interpretated as a circle)
# (Position means position on the genome)
#
#
# CASE_2.1:
# IF the starting position of the seed minus the needed leftside of the fragment is surpassing the beginnboarder of the geonme (here pos = 0)
# There's an overlay to be calculated;
# If we're surpassing the boarder with the leftside, we need to get the "rest" at the oposite side of the genome.
# The overlay is the length of the needed fragment (default fragmentsize =1250) minus the part we already generated on the leftside
#
#
# CASE_2.2:
# IF the end point of the seed plus the rightside of the fragment is not supassing the endboarder of the geonme the searched fragment is in the middle
# and can be cutted out.
#
#
# CASE_2.3: Seed is in the middle of the linear genome
# IF the starting position of the seed minus the needed leftside of the fragment is not surpassing the beginnboarder of the geonme (here pos = 0)
# &&
# IF the end point of the seed plus the rightside of the fragment is not supassing the endboarder of the geonme the searched fragment is in the middle
# and can be cutted out. If not so, the fragment can't be created because the genome isn't circular!
#
if((begin<0)&&(end<=genome_length)){ #Seed is over the left side of the genome
right_fragment <- substr(genome_sequence,end+1,end)
seed <- substr(genome_sequence,begin,end)
left_fragment <- substr(genome_sequence,1,begin-1)
overlay <- abs(width(left_fragment))
overlay_left_fragment <- substr(genome_sequence,width(genome_sequence)-(overlay-1),width(genome_sequence))
fragment <-paste(overlay_left_fragment,left_fragment,seed,right_fragment,sep="")
if(!(grepl("[^ACGT]",fragment))){fragment_dataframe[i,6]<- fragment}
}else if((begin>=0)&&(end>genome_length)){ #Seed is over the right side of the genome
left_fragment <- substr(genome_sequence,begin,begin-1)
seed <- substr(genome_sequence,begin,end)
right_fragment <- substr(genome_sequence,end+1,width(genome_sequence))
overlay <- abs(width(right_fragment))
overlay_right_fragment <- substr(genome_sequence,1,(overlay))
fragment <-paste(left_fragment,seed,right_fragment,overlay_right_fragment,sep="")
if(!(grepl("[^ACGT]",fragment))){fragment_dataframe[i,6]<- fragment}
}else if((begin>=0)&&(end<=genome_length)){
fragment <-substr(genome_sequence, begin, end)
if(!(grepl("[^ACGT]",fragment))){fragment_dataframe[i,6]<- fragment}
}
}
}
remove_rows = which(fragment_dataframe[[6]] == "FALSE")
if (length(remove_rows) != 0){
fragment_dataframe = fragment_dataframe[-remove_rows, ]
}
# save the fragmentdataframe to a FASTA-file; needed to call the LS-GKM algorithm
return(fragment_dataframe)
}
write_fastafile_before_lsgkm <- function(fragment_dataframe, path_to_fragmentfile){
for (i in 1:nrow(fragment_dataframe)){
fragment_dataframe_row<- as.matrix(as.matrix(fragment_dataframe)[i,])
name <- unname(fragment_dataframe_row["identifier", ])
fragment <-  unname(fragment_dataframe_row["fragment", ])
write(paste(">",i,"/", name, "\n", fragment, sep=""), file=path_to_fragmentfile, append = T)
}
}
# Input:
# List of fragments from get_fragments_of_size_n()
# saves the list of fragments into a FASTAfile and Will predict the propabilitly of the fragment being an ori
# by using lsgkm (N1250_L10_K6_D4_T4).
# Reads out the predvalues pred by lsgkm method writes them into a File
# Output:
# fragment_dataframe - List of framents and their predictionvalue
#
#
predict_fragments <- function(fragment_dataframe_input, path_to_lsgkmpredict, path_to_model, path_to_fragmentdataframe, path_to_output_result,
fragmentsize = 1250, linear = TRUE, cores = 4){
lsgkm_command =  paste(path_to_lsgkmpredict, path_to_fragmentdataframe, path_to_model, path_to_output_result, "-T", cores, sep = " ")
system(lsgkm_command)
# write the newest data of prediction by lsgkm into the fragment_dataframe
fragmentpredictionlist <- read.csv(path_to_output_result, sep ="\t", header = F)
fragment_dataframe_input[["raw_output"]] <- fragmentpredictionlist[, 2]
return(fragment_dataframe_input)
}
# Scaling =======
# Input:
# Dataframe of fragments
# normalizes the predicitionvalues from lsgkm by:
#     - adding the absolute value of the minimum to all elements.
#     - dividing all values by the absolute value of the maxmimum.
#     => scales the values into [0,1]
# Output:
# fragment_dataframe - column "calibrated" filled with calibrated values between [0,1] of lsgkm-model-prediction
#
#
normalize<-function(fragment_dataframe){
abs_min <-abs(min(fragment_dataframe[["raw_output"]], na.rm = T))
for(i in 1:nrow(fragment_dataframe)){
fragment_dataframe[i, "normalized"]<- (fragment_dataframe[i,  "raw_output"]+abs_min)
}
abs_max <- abs(max(fragment_dataframe[["normalized"]], na.rm = T))
for(i in 1:nrow(fragment_dataframe)){
fragment_dataframe[i, "normalized"]<- fragment_dataframe[i, "normalized"]/abs_max
}
return(fragment_dataframe)
}
### Cutoff ========
# Input:
# Dataframe of fragments
# decimates the predicitionvalues from lsgkm by "throwing away" all values less than the cut_value
# Output:
# fragment_dataframe - only columns where the value in "calibrated" is equal or bigger than the cutvalue
# #TODO correct this documentation
#
apply_cutoffs <- function(fragment_dataframe, lower.cutoff = 0.41, upper.cutoff = 0.99){
rejected_dataframes <- fragment_dataframe[which(fragment_dataframe[["normalized"]] > lower.cutoff & fragment_dataframe[["normalized"]] < upper.cutoff), ]
positive_dataframes <- fragment_dataframe[which(fragment_dataframe[["normalized"]] >= upper.cutoff), ]
return(list("positive" = positive_dataframes, "rejected" = rejected_dataframes) )
}
### CalibratR ========
#TODO Document
calibrate_results <- function(fragment_dataframe, calibratR_model_path){
library(CalibratR)
library(R.filesets)
calib_model <- loadRDS(file = calibratR_model_path)
naindexes = which(is.na(fragment_dataframe[[7]]))
fragment_dataframe[[7]][naindexes] = 0
calib_preds <- predict_calibratR(calib_model$calibration_models, fragment_dataframe[[7]])
fragment_dataframe[[7]][naindexes] = NaN
fragment_dataframe[["Calibrated"]] = calib_preds$hist_scaled
fragment_dataframe[["Calibrated"]] = NaN
return(fragment_dataframe)
}
### Output =========
# Input:
# (selected) Dataframe of fragments and url to create the output
# Output:
# fragment_dataframe
#
#
create_output <- function(selected_fragment_dataframe, target_fasta, target_dataframe, target_ambiguous_dataframe, calibrated = F){
write.csv(selected_fragment_dataframe$positive, target_dataframe)
write.csv(selected_fragment_dataframe$rejected, target_ambiguous_dataframe)
print(selected_fragment_dataframe$positive)
if (calibrated){
NA#TODO implement, if CalibratR is implemented
}else{
tmp <- apply(selected_fragment_dataframe$positive, 1, function(x) write(paste(">", x[["identifier"]], " ", x[["begin"]], ":", x[["end"]], " ", x[["raw_output"]], " ",
x[["normalized"]], "\n", x[["fragment"]], sep=""), file = target_fasta, append = T))
}
}
isValidEmail <- function(x) {
grepl("\\<[A-Z0-9._%+-]+@[A-Z0-9.-]+\\.[A-Z]{2,}\\>", as.character(x), ignore.case=TRUE)
}
writeEmail<-function(from,to,subject,message,smtp){
require('mailR')
library(mailR)
#TODO write something that sends the output of BOriS to the user. use send.mail or smth
}
### Programmcall =========
#
# Input:
# The FASTA-File paths to data/seeds/lsgkmpredict/model/fragmentdataframe/result,
# The wanted fragmentsize = 1250 and if the genome is linear = TRUE
#
# Output:
# The output will be a FASTA-File with the ori-sequences, their positions and the probability of beeing an ori
#
# Example output:
# >Name_of_sequence    position_begin_seed:posisiton_end_seed    probability (non_calibrated) probability (calibrated)
#  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAATTTTTTTTTTTTTTTTTTTTTGGGGGGGGGGGGGGGGGGGGGGGGCCCCCCCCCCCCCCCCC
#
#
call_BOriS <- function(path_to_data, path_to_seeds, path_to_lsgkmpredict, path_to_model, path_to_output_result, calibratR_model_path,
fragmentsize = 1250, linear = TRUE, cores = 4, single_chromosome = TRUE){#TODO implement single_chromosomes
library(stringr)
library(Biostrings)
print("Started BOriS...")
data <- readDNAStringSet(path_to_data)
data <- rename_FASTA_ID(data) #TODO do I need to cut off the description?
seeds<- readDNAStringSet(path_to_seeds)
print("Creating fragments")
fragment_path = paste(path_to_output_result, "/fragments_", as.numeric(Sys.time())*1000, ".fa", sep = "")
# get the fragments to predict their possibilies to be an Ori by LS-GKM -- loop is necessary due to the fact that multiple sequences might be in data
for (counter in c(1:length(data))){
if (counter == 1){
fragment_dataframe <- get_fragments_of_size_n(data[counter], seeds, fragment_path, n = fragmentsize, linear = linear)
}else{
tmp <- get_fragments_of_size_n(data[counter], seeds, fragment_path, n = fragmentsize, linear = linear)
fragment_dataframe <- rbind(fragment_dataframe, tmp)
}
}
write_fastafile_before_lsgkm(fragment_dataframe, fragment_path)
print("Calling LS-GKM...")
tmp_path = paste(path_to_output_result, "/tmp_", as.numeric(Sys.time())*1000, ".fa", sep = "")
fragment_dataframe <- predict_fragments(fragment_dataframe, path_to_lsgkmpredict, path_to_model, fragment_path, tmp_path,
fragmentsize = fragmentsize, linear = linear, cores = cores)
invisible(file.remove(tmp_path))
invisible(file.remove(fragment_path))
print("Done calling LS-GKM...")
fragment_dataframe["orientation"] = "+"
rv_data <- reverseComplement(data)
print("Creating reverse complement fragments")
fragment_path = paste(path_to_output_result, "/rev_fragments_", as.numeric(Sys.time())*1000, ".fa", sep = "")
# get the fragments to predict their possibilies to be an Ori by LS-GKM
for (counter in c(1:length(rv_data))){
if (counter == 1){
rev_fragment_dataframe <- get_fragments_of_size_n(rv_data[counter], seeds, fragment_path, n = fragmentsize, linear = linear)
}else{
tmp <- get_fragments_of_size_n(rv_data[counter], seeds, fragment_path, n = fragmentsize, linear = linear)
rev_fragment_dataframe <- rbind(rev_fragment_dataframe, tmp)
}
}
write_fastafile_before_lsgkm(rev_fragment_dataframe, fragment_path)
print("Calling LS-GKM...")
tmp_path = paste(path_to_output_result, "/tmp_", as.numeric(Sys.time())*1000, ".fa", sep = "")
rev_fragment_dataframe <- predict_fragments(rev_fragment_dataframe, path_to_lsgkmpredict, path_to_model, fragment_path, tmp_path,
fragmentsize = fragmentsize, linear = linear, cores = cores)
invisible(file.remove(tmp_path))
invisible(file.remove(fragment_path))
print("Done calling LS-GKM...")
rev_fragment_dataframe["orientation"] = "-"
full_dataframe <- rbind(fragment_dataframe, rev_fragment_dataframe)
if (single_chromosome){
full_dataframe <- normalize(full_dataframe)
selected_fragments <- apply_cutoffs(full_dataframe, 0.41, 0.99)
}else{
chromosome_ids <- unique(as.character(full_dataframe[["identifier"]]))
for (counter in c(1:length(chromosome_ids))){
identifier_fragments <- full_dataframe[which(full_dataframe[["identifier"]] == chromosome_ids[counter]), ]
identifier_fragments <- normalize(identifier_fragments)
if (counter == 1){
selected_fragments <- apply_cutoffs(identifier_fragments, 0.41, 0.99)
}else{
selected_fragments_identifier <- apply_cutoffs(identifier_fragments, 0.41, 0.99)
positives_tmp <- rbind(selected_fragments$positive, selected_fragments_identifier$positive)
rejected_tmp <- rbind(selected_fragments$rejected, selected_fragments_identifier$rejected)
selected_fragments <- list("positive"=positives_tmp, "rejected"=rejected_tmp)
}
}
}
print("Normalized and Selected LS-GKM output...")
#print("Calibrating results...")
#if (calibratR_model_path == ""){ # TODO think about whether to use this
#  full_dataframe <- calibrate_results(full_dataframe, calibratR_model_path)
#}
print("Printing results to files...")
create_output(selected_fragments,
paste(path_to_output_result, paste("BOriS_", strsplit(strsplit(path_to_data, "/")[[1]]
[length(strsplit(path_to_data, "/")[[1]])], "\\.fa")[[1]][1], "_fa.fasta", sep=""), sep = "/"),
paste(path_to_output_result, paste("BOriS_", strsplit(strsplit(path_to_data, "/")[[1]]
[length(strsplit(path_to_data, "/")[[1]])], "\\.fa")[[1]][1], "_df.csv", sep=""), sep = "/"),
paste(path_to_output_result, paste("BOriS_", strsplit(strsplit(path_to_data, "/")[[1]]
[length(strsplit(path_to_data, "/")[[1]])], "\\.fa")[[1]][1], "_ambiguous_df.csv", sep=""), sep = "/")
)
print("Done!")
}
data <- read.csv("/home/sperlea/projects/publications/BOriS/data/BOriS_NCBI/assembly_summary.txt", sep = "\t")
data <- data[which(data[["excluded_from_refseq"]] == ""), ]
data <- data[which(data[["release_type"]] == "Major"), ]
data_Chromosome <- data[which(data[["assembly_level"]] == "Chromosome"), ]
data_Complete <- data[which(data[["assembly_level"]] == "Complete Genome"), ]
data_Contig <- data[which(data[["assembly_level"]] == "Contig"), ]
data_Scaffold <- data[which(data[["assembly_level"]] == "Scaffold"), ]
data_Chromosome[1, ]
?gunzip
call_BOriS_with_ftp_file <- function(ftpurl, source, target){
library(RCurl)
library(R.utils)
filenames = strsplit(getURL(ftpurl, dirlistonly = TRUE), "\n")[[1]]
correctname <- filenames[grep("_genomic.fna.gz", filenames)][-grep("_from_", filenames[grep("_genomic.fna.gz", filenames)])]
print(correctname)
print(paste(target, correctname, sep = ""))
download.file(paste(ftpurl, correctname, sep = ""), paste(target, correctname, sep = ""), quiet=T)
gunzip(paste(target, correctname, sep = ""))
newfilename = sub(".gz", "", paste(target, correctname, sep = ""))
call_BOriS(path_to_data = newfilename,
path_to_seeds = paste(source, "seeds.fa", sep = ""),
path_to_lsgkmpredict = "gkmpredict",
path_to_model = paste(source, "ls-gkm_train_model_1250_L10_K6_D4_T4.model.txt", sep = ""),
path_to_output_result = target,
calibratR_model_path = paste(source, "ls-calibratR_model_1250_train_selected.RDS", sep = ""),
1250, FALSE, single_chromosome = TRUE)
}
ff <- "ffefe"
substr(ff, nchar(ff))
substr(ff, nchar(ff)-1, nchar(ff))
substr(ff, nchar(ff), nchar(ff))
substr(ff, nchar(ff), nchar(ff)) == "/"
data <- read.csv("/home/sperlea/projects/publications/BOriS/data/BOriS_NCBI/assembly_summary.txt", sep = "\t")
data <- data[which(data[["excluded_from_refseq"]] == ""), ]
data <- data[which(data[["release_type"]] == "Major"), ]
data_Chromosome <- data[which(data[["assembly_level"]] == "Chromosome"), ]
shiny::runApp('projects/publications/BOriS/programs/BOriS_v0.6')
runApp('projects/publications/BOriS/programs/BOriS_v0.6')
runApp('projects/publications/BOriS/programs/BOriS_v0.6')
runApp('projects/publications/BOriS/programs/BOriS_v0.6')
runApp('projects/publications/BOriS/programs/BOriS_v0.6')
call_BOriS_from_shiny <- function(path, linear, single_chromosome){
##these paths will need to be changed!!!
#path = "/home/sperlea/projects/publications/BOriS/data/genomes_tmp/NC_000913.3.fa"
path_to_seeds = "/home/sperlea/projects/publications/BOriS/programs/BOriS_v0.6/seeds.fa"
path_to_lsgkmpredict = "/home/sperlea/projects/publications/BOriS/programs/BOriS_v0.6/gkmpredict"
path_to_model = "/home/sperlea/projects/publications/BOriS/programs/BOriS_v0.6/ls-gkm_train_model_1250_L10_K6_D4_T4.model.txt"
path_to_output_result = "/home/sperlea/projects/publications/BOriS/programs/BOriS_v0.6/ff"
fragmentsize = 1250
linear = FALSE
cores = 4
#call_BOriS(path, path_to_seeds, path_to_lsgkmpredict, path_to_model, path_to_output_result,
#                       fragmentsize, linear, cores, single_chromosome)
}
runApp('projects/publications/BOriS/programs/BOriS_v0.6')
runApp('projects/publications/BOriS/programs/BOriS_v0.6')
runApp('projects/publications/BOriS/programs/BOriS_v0.6')
runApp('projects/publications/BOriS/programs/BOriS_v0.6')
runApp('projects/publications/BOriS/programs/BOriS_v0.6')
runApp('projects/publications/BOriS/programs/BOriS_v0.6')
runApp('projects/publications/BOriS/programs/BOriS_v0.6')
runApp('projects/publications/BOriS/programs/BOriS_v0.6')
runApp('projects/publications/BOriS/programs/BOriS_v0.6')
runApp('projects/publications/BOriS/programs/BOriS_v0.6')
