daily pastebin goal
44%
SHARE
TWEET

Untitled

a guest May 16th, 2018 111 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. library("seqinr")
  2. library("Biostrings")
  3. library("MASS")
  4. library("GenomicRanges")
  5. library("dndscv")
  6.  
  7. ## Helper methods copied from dndscv
  8. nt = c("A","C","G","T")
  9. compnt = setNames(rev(nt), nt)
  10.  
  11. chr2cds = function(pos,cds_int,strand) {
  12.   if (strand==1) {
  13.     return(which(pos==unlist(apply(cds_int, 1, function(x) x[1]:x[2]))))
  14.   } else if (strand==-1) {
  15.     return(which(pos==rev(unlist(apply(cds_int, 1, function(x) x[1]:x[2])))))
  16.   }
  17. }
  18.  
  19. rest_nt <- function(n) {
  20.   rest_nts <- nt[-which(nt == n)]
  21.  
  22.   sample(rest_nts, 1)
  23. }
  24.  
  25. data("refcds_hg19", package="dndscv")
  26.  
  27. TP53 <- RefCDS[which(sapply(RefCDS, function(x) x$gene_name) == "TP53")]
  28.  
  29. cds_int <- TP53[[1]]$intervals_cds
  30. cds_pos <- unlist(apply(cds_int, 1, function(x) x[1]:x[2]))
  31. splice_pos <- TP53[[1]]$intervals_splice
  32.  
  33. # Use mutations from example dataset but remove all mapping to TP53
  34. data("dataset_simbreast", package="dndscv")
  35. my_mutations <- subset(mutations, !(chr == "17" & pos %in% cds_pos) & nchar(ref) == 1 & !(pos %in% splice_pos))
  36.  
  37. # Sample 50 random positions in the TP53 CDS
  38. my_mut_pos <- sample(cds_pos, 50)
  39.  
  40. # Make a data frame with random mutations at sampled positions
  41. my_mut_cds_pos <- sapply(my_mut_pos, function(x) (chr2cds(x, TP53[[1]]$intervals_cds, TP53[[1]]$strand)))
  42. my_mut_ref <- sapply(my_mut_cds_pos,
  43.                      function(x) as.character(TP53[[1]]$seq_cds[x]))
  44. my_muts <- data.frame(sampleID=paste0("S_", 1:length(my_mut_pos)),
  45.                       chr=TP53[[1]]$chr,pos=my_mut_pos,
  46.                       ref=compnt[my_mut_ref],
  47.                       mut=sapply(my_mut_ref, rest_nt),
  48.                       stringsAsFactors=F)
  49.  
  50. ## Apply dndscv on the combined mutations
  51. dndsout = dndscv(rbind(my_mutations[, 1:5], my_muts))
  52.  
  53. sel_cv = dndsout$sel_cv
  54. signif_genes = sel_cv[sel_cv$qglobal_cv<0.1, c("gene_name","qglobal_cv")]
  55. rownames(signif_genes) = NULL
  56. print(signif_genes)
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top