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

Get V and J sequences from sequence ID

This commit is contained in:
coolneng 2021-04-21 21:29:03 +02:00
parent fb5d781c66
commit 659f0097d8
Signed by: coolneng
GPG Key ID: 9893DA236405AF57
1 changed files with 33 additions and 31 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

64
src/alignment.r
View File

@ -1,36 +1,39 @@
library(Biostrings) library(Biostrings)
library(parallel) library(parallel)
parse_data <- function(files) { parse_data <- function(file) {
reversed_sequences <- Biostrings::readQualityScaledDNAStringSet(files[1]) reversed_sequences <- Biostrings::readQualityScaledDNAStringSet(file)
sequences <- Biostrings::reverseComplement(reversed_sequences) sequences <- Biostrings::reverseComplement(reversed_sequences)
vdj_metadata <- read.csv(files[2])
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_metadata)) return(list(sequences, vj_segments))
} }
get_vj_sequence <- function(identifier, names, sequences) { match_id_sequence <- function(names, sequences, id) {
matches <- grep(names, pattern = identifier) matches <- grep(names, pattern = id)
row <- matches[1] row <- matches[1]
return(as.character(sequences[row])) return(as.character(sequences[row]))
} }
construct_full_sequences <- function(vdj_segments, metadata) { get_vj_sequence <- function(sequences, names, vdj_segments) {
v_sequences <- lapply(metadata$v_call, metadata <- unlist(strsplit(sequences, split = " "))
v_identifier <- metadata[2]
j_identifier <- metadata[3]
v_sequence <- match_id_sequence(names, vdj_segments, id = v_identifier)
j_sequence <- match_id_sequence(names, vdj_segments, id = j_identifier)
return(c(v_sequence, j_sequence))
}
fetch_vj_sequences <- function(vdj_segments, sequences) {
vj_sequences <- mclapply(names(sequences),
names(vdj_segments), names(vdj_segments),
vdj_segments, vdj_segments,
FUN = get_vj_sequence FUN = get_vj_sequence,
mc.cores = detectCores()
) )
j_sequences <- lapply(metadata$j_call, return(c(vj_sequences[1], vj_sequences[2]))
names(vdj_segments),
vdj_segments,
FUN = get_vj_sequence
)
full_sequence <- paste(v_sequences, metadata$junction, j_sequences, sep = "")
return(Biostrings::DNAStringSet(full_sequence))
} }
align_sequence <- function(sequence, vdj_segment) { align_sequence <- function(sequence, vdj_segment) {
@ -38,25 +41,24 @@ align_sequence <- function(sequence, vdj_segment) {
subject = sequence, subject = sequence,
pattern = vdj_segment, pattern = vdj_segment,
type = "global-local", type = "global-local",
gapOpening = 1, gapOpening = 1
)) ))
} }
perform_alignment <- function(sequences, vdj_segments, metadata) { get_hvr_sequences <- function(sequences, vdj_segments) {
vj_sequences <- construct_full_sequences(vdj_segments, metadata) vj_sequences <- fetch_vj_sequences(vdj_segments, sequences)
sequence_alignment <- mcmapply(vj_sequences, v_alignment <- parallel::mcmapply(sequences,
vdj_segments, vj_sequences[1],
FUN = align_sequence, FUN = align_sequence
mc.cores = detectCores() )
j_alignment <- parallel::mcmapply(sequences,
vj_sequences[2],
FUN = align_sequence
) )
return(sequence_alignment)
} }
input_files <- c("data/curesim_sequence.fastq", "data/vdj_metadata.csv") data <- parse_data(file = "data/curesim_sequence.fastq")
data <- parse_data(files = input_files) hvr_sequences <- get_hvr_sequences(
alignment <- perform_alignment(
sequences = data[[1]], sequences = data[[1]],
vdj_segments = data[[2]], vdj_segments = data[[2]]
metadata = data[[3]] )
)
print(alignment)