Склеивание строк из длинной ссылочной строки с учетом их координат в этой ссылочной строке - PullRequest
Новая
       11

Склеивание строк из длинной ссылочной строки с учетом их координат в этой ссылочной строке

0 голосов
/ 26 января 2019

Хотя это вопрос genomics, поскольку он имеет дело с объединением (получением подмножеств) строк, я думаю, что это актуально для этой аудитории, а не только Bioconductor.

Довольно просто, у меня есть списокдлинных струн (хромосомы генома).Например, я создаю и сохраняю 10 хромосом, используя пакет Bioconductor Biostrings: 

set.seed(1)
set <- NULL
for (i in 1:10) set <- c(set,paste(sample(Biostrings::DNA_ALPHABET[1:4],10000,replace=T),collapse=""))

genome.set <- Biostrings::DNAStringSet(set)
names(genome.set) <- paste0("chr",1:10)

И затем у меня есть data.frame координат транскрипта (из файла GTF), где каждыйТранскрипт может иметь несколько строк: 

library(dplyr)
gtf.df <- data.frame(seqnames = sample(names(genome.set),100,replace=T),
                     strand = sample(c("+","-"),100,replace=T),
                     start = sample(1:9000,100,replace=F)) %>%
  dplyr::mutate(end = start+sample(1:1000,100,replace = F))

gtf.df <- gtf.df %>% dplyr::group_by(seqnames) %>%
  dplyr::arrange(start,end) %>%
  dplyr::mutate(transcript_id = paste0(seqnames,"_",sample(1:8,length(seqnames),replace=T))) %>%
  dplyr::ungroup()

И я хочу, чтобы каждый транскрипт соединял свои последовательности, выделяя их из genome.set.

Повторное использование Biostrings Iдобиться этого следующим образом: 

transcript_ids <- unique(gtf.df$transcript_id)
transcript.seqs <- sapply(1:length(transcript_ids),function(t){
  transcript.gtf.df <- gtf.df %>% dplyr::filter(gtf.df$transcript_id == transcript_ids[t])
  transcript.seq <- paste(sapply(1:nrow(transcript.gtf.df),function(e)
    unname(as.character(Biostrings::subseq(genome.set[which(names(genome.set) == transcript.gtf.df$seqnames[1])],start=transcript.gtf.df$start[e],end=transcript.gtf.df$end[e])))
  ),collapse="")
  if(transcript.gtf.df$strand[1] == "-") transcript.seq <- unname(as.character(Biostrings::reverseComplement(Biostrings::DNAString(transcript.seq))))
  return(transcript.seq)
})

Моя проблема в том, что у меня есть 4520919 стенограммы в моих реальных данных, и последняя часть занимает много времени.Поэтому мой вопрос заключается в том, можно ли и как это сделать быстрее, используя Biostrings или любым другим способом.

1 Ответ

0 голосов
/ 26 января 2019

Я переписал ваши sapply методы с двумя основными изменениями:

  • Во-первых, я использовал vapply, что в целом гораздо быстрее
  • Во-вторых, я использовал многоот .subset2 до поднабора фреймов данных

Редактировать

  • Мне удалось избавиться от внутреннего цикла (vapply)
  • Заменена функция Biostrings::reverseComplement

Вот код 

names_genome.set <- names(genome.set)
transcript_ids <- unique(gtf.df$transcript_id)
transcript_seqs <- vapply(seq_along(transcript_ids), function (t) {
  ind_id <- which(.subset2(gtf.df, 5L) == transcript_ids[t])
  x <- unname(as.character(genome.set[names_genome.set == .subset2(gtf.df, 1L)[ind_id[1L]]]))
  out <- paste0(substring(text = x, first = .subset2(gtf.df, 3L)[ind_id], last = .subset2(gtf.df, 4L)[ind_id]), collapse = '')
  if (.subset2(gtf.df, 2L)[ind_id[1L]] == '-') {
    out <- unlist(strsplit(out, ''))
    ind_A <- out == 'A'
    ind_T <- out == 'T'
    ind_C <- out == 'C'
    ind_G <- out == 'G'
    out[ind_A] <- 'C'
    out[ind_T] <- 'G'
    out[ind_G] <- 'T'
    out[ind_C] <- 'A'
    out <- paste(out, collapse = '')
  }
  out
}, character(1))

Вот некоторые контрольные цифры с предоставленными образцами данных 

# Unit: milliseconds
#       expr       min        lq      mean    median        uq      max neval cld
#     sapply 160.94296 169.97698 180.13836 175.20474 182.58224 400.3273   100   c
# vapply_old  66.20113  69.59185  72.96804  71.45861  74.56051 120.0023   100  b 
# vapply_new  47.45224  49.51573  52.87001  50.97023  54.52104 109.3320   100 a  

microbenchmark::microbenchmark(
  'sapply' = {
    transcript.seqs <- sapply(1:length(transcript_ids),function(t){
      transcript.gtf.df <- gtf.df %>% dplyr::filter(gtf.df$transcript_id == transcript_ids[t])
      transcript.seq <- paste(sapply(1:nrow(transcript.gtf.df),function(e)
        unname(as.character(Biostrings::subseq(genome.set[which(names(genome.set) == transcript.gtf.df$seqnames[1])],start=transcript.gtf.df$start[e],end=transcript.gtf.df$end[e])))
      ),collapse="")
      if(transcript.gtf.df$strand[1] == "-") transcript.seq <- unname(as.character(Biostrings::reverseComplement(Biostrings::DNAString(transcript.seq))))
      return(transcript.seq)
    })
  },
  'vapply_old' = {
    transcript_seqs <- vapply(seq_along(transcript_ids), function (t) {
      ind_id <- which(.subset2(gtf.df, 5L) == transcript_ids[t])
      x <- unname(as.character(genome.set[names_genome.set == .subset2(gtf.df, 1L)[ind_id[1L]]]))
      out <- vapply(ind_id, 
                    function (e) substr(x = x, start = .subset2(gtf.df, 3L)[e], stop = .subset2(gtf.df, 4L)[e]),
                    character(1))
      out <- paste(out, collapse = '')
      if (.subset2(gtf.df, 2L)[ind_id[1L]] == '-') {
        out <- unname(as.character(Biostrings::reverseComplement(Biostrings::DNAString(out))))
      }
      out
    }, character(1))
  },
  'vapply_new' = {
    transcript_seqs <- vapply(seq_along(transcript_ids), function (t) {
      ind_id <- which(.subset2(gtf.df, 5L) == transcript_ids[t])
      x <- unname(as.character(genome.set[names_genome.set == .subset2(gtf.df, 1L)[ind_id[1L]]]))
      out <- paste0(substring(text = x, first = .subset2(gtf.df, 3L)[ind_id], last = .subset2(gtf.df, 4L)[ind_id]), collapse = '')
      if (.subset2(gtf.df, 2L)[ind_id[1L]] == '-') {
        out <- unlist(strsplit(out, ''))
        ind_A <- out == 'A'
        ind_T <- out == 'T'
        ind_C <- out == 'C'
        ind_G <- out == 'G'
        out[ind_A] <- 'C'
        out[ind_T] <- 'G'
        out[ind_G] <- 'T'
        out[ind_C] <- 'A'
        out <- paste(out, collapse = '')
      }
      out
    }, character(1))
  }
)

Я постараюсь найти способы еще улучшить его (возможна векторизация). Я еще не понял, к примеру, что делает функция reverseComplement - возможно, это может быть выполнено более эффективно.

Вы можете попробовать с большими наборами данных и посмотреть, есть ли улучшения,Также, Rcpp может быть идеей, если эффективность действительно поставлена ​​на карту.

...