2021-02-21 20:59:39 +01:00
|
|
|
library(immuneSIM)
|
2021-02-26 02:19:40 +01:00
|
|
|
library(Biostrings)
|
2021-02-21 20:59:39 +01:00
|
|
|
|
2021-03-02 20:08:14 +01:00
|
|
|
generate_repertoire <- function(number_of_sequences) {
|
2021-03-10 12:34:20 +01:00
|
|
|
return(immuneSIM(
|
2021-02-25 20:00:35 +01:00
|
|
|
number_of_seqs = number_of_sequences,
|
|
|
|
species = "hs",
|
|
|
|
receptor = "tr",
|
2021-03-10 12:34:20 +01:00
|
|
|
chain = "b"
|
|
|
|
))
|
2021-02-25 20:00:35 +01:00
|
|
|
}
|
2021-02-21 20:59:39 +01:00
|
|
|
|
2021-02-28 02:23:58 +01:00
|
|
|
# TODO save also v_call and j_call
|
|
|
|
preprocess_data <- function(repertoire, sequencing_runs) {
|
2021-02-26 02:19:40 +01:00
|
|
|
sequences <- as.character(repertoire$sequence)
|
2021-02-28 02:23:58 +01:00
|
|
|
reads <- Biostrings::DNAStringSet(rep(sequences, sequencing_runs))
|
2021-02-26 02:19:40 +01:00
|
|
|
names(reads) <- seq_len(length(reads))
|
|
|
|
reverse_complement <- Biostrings::reverseComplement(reads)
|
|
|
|
return(reverse_complement)
|
|
|
|
}
|
|
|
|
|
2021-02-28 02:23:58 +01:00
|
|
|
save_data <- function(repertoire) {
|
2021-03-02 20:08:14 +01:00
|
|
|
file_name <- "data/sequence.fasta"
|
|
|
|
Biostrings::writeXStringSet(repertoire, file_name, format = "fasta")
|
2021-02-26 02:20:11 +01:00
|
|
|
}
|
2021-02-21 20:59:39 +01:00
|
|
|
|
2021-02-25 20:02:11 +01:00
|
|
|
parse_cli_arguments <- function(args) {
|
2021-02-28 02:23:58 +01:00
|
|
|
if (length(args) != 2) {
|
|
|
|
stop("usage: repertoire.r <number of sequences> <sequencing_runs>")
|
2021-02-25 20:02:11 +01:00
|
|
|
}
|
2021-02-28 02:23:58 +01:00
|
|
|
return(c(args[1], args[2]))
|
2021-02-25 20:02:11 +01:00
|
|
|
}
|
|
|
|
|
|
|
|
args <- commandArgs(trailingOnly = TRUE)
|
2021-03-02 20:08:14 +01:00
|
|
|
arguments <- parse_cli_arguments(args)
|
|
|
|
number_of_sequences <- as.integer(arguments[1])
|
|
|
|
sequencing_runs <- as.integer(arguments[2])
|
|
|
|
repertoire <- generate_repertoire(number_of_sequences)
|
2021-02-28 02:23:58 +01:00
|
|
|
processed_data <- preprocess_data(repertoire, sequencing_runs)
|
2021-02-26 02:20:11 +01:00
|
|
|
save_data(processed_data)
|