coolneng
/
locigenesis
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Remove redundant HVR sequence construction

This commit is contained in:
coolneng 2021-04-07 19:49:44 +02:00
parent e694ee3292
commit dd9f7ffde4
Signed by: coolneng
GPG Key ID: 9893DA236405AF57
2 changed files with 8 additions and 24 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

22
src/alignment.r
View File

@ -1,27 +1,15 @@
library(Biostrings) library(Biostrings)
library(parallel) library(parallel)
construct_hvr_sequence <- function(data) {
vdj_dataframe <- as.data.frame(data)
hvr_sequence <- paste(vdj_dataframe$v_sequence,
vdj_dataframe$d_sequence,
vdj_dataframe$j_sequence,
sep = ""
)
vdj_dataframe$hvr <- hvr_sequence
return(vdj_dataframe)
}
parse_data <- function(files) { parse_data <- function(files) {
reversed_sequences <- Biostrings::readQualityScaledDNAStringSet(files[1]) reversed_sequences <- Biostrings::readQualityScaledDNAStringSet(files[1])
sequences <- Biostrings::reverseComplement(reversed_sequences) sequences <- Biostrings::reverseComplement(reversed_sequences)
vdj_alignment <- read.csv(files[2]) vdj_metadata <- read.csv(files[2])
vdj_dataframe <- construct_hvr_sequence(vdj_alignment)
vj_segments <- union( vj_segments <- union(
readRDS("data/v_segments.rds"), readRDS("data/v_segments.rds"),
readRDS("data/j_segments_phe.rds") readRDS("data/j_segments_phe.rds")
) )
return(list(sequences, vj_segments, vdj_dataframe)) return(list(sequences, vj_segments, vdj_metadata))
} }
get_vj_sequence <- function(identifier, names, sequences) { get_vj_sequence <- function(identifier, names, sequences) {
@ -41,7 +29,7 @@ construct_full_sequences <- function(vdj_segments, metadata) {
vdj_segments, vdj_segments,
FUN = get_vj_sequence FUN = get_vj_sequence
) )
full_sequence <- paste(v_sequences, metadata$hvr, j_sequences, sep = "") full_sequence <- paste(v_sequences, metadata$junction, j_sequences, sep = "")
return(Biostrings::DNAStringSet(full_sequence)) return(Biostrings::DNAStringSet(full_sequence))
} }
@ -64,8 +52,8 @@ perform_alignment <- function(sequences, vdj_segments, metadata) {
return(sequence_alignment) return(sequence_alignment)
} }
input_files <- c("data/curesim_sequence.fastq", "data/vdj_alignment.csv") input_files <- c("data/curesim_sequence.fastq", "data/vdj_metadata.csv")
data <- parse_data(input_files) data <- parse_data(files = input_files)
alignment <- perform_alignment( alignment <- perform_alignment(
sequences = data[[1]], sequences = data[[1]],
vdj_segments = data[[2]], vdj_segments = data[[2]],

8
src/repertoire.r
View File

@ -14,15 +14,11 @@ save_data <- function(data, reads) {
Biostrings::writeXStringSet(data$sequence, "data/sequence.fasta") Biostrings::writeXStringSet(data$sequence, "data/sequence.fasta")
vdj_sequences <- data[-1] vdj_sequences <- data[-1]
amplified_vdj <- vdj_sequences[rep(seq_len(nrow(vdj_sequences)), reads), ] amplified_vdj <- vdj_sequences[rep(seq_len(nrow(vdj_sequences)), reads), ]
write.csv(amplified_vdj, "data/vdj_alignment.csv", row.names = FALSE) write.csv(amplified_vdj, "data/vdj_metadata.csv", row.names = FALSE)
} }
process_data <- function(repertoire, reads) { process_data <- function(repertoire, reads) {
columns <- c( columns <- c("sequence", "junction", "v_call", "j_call")
"sequence", "v_sequence_alignment",
"d_sequence_alignment", "j_sequence_alignment",
"v_call", "j_call"
)
data <- repertoire[, columns] data <- repertoire[, columns]
dna_sequence <- Biostrings::DNAStringSet(data$sequence) dna_sequence <- Biostrings::DNAStringSet(data$sequence)
data$sequence <- Biostrings::reverseComplement(dna_sequence) data$sequence <- Biostrings::reverseComplement(dna_sequence)