From 120cb2a7c8b1c36f830b7928d30c6d74c2129f54 Mon Sep 17 00:00:00 2001
From: Laurent Modolo <laurent@modolo.fr>
Date: Mon, 8 Oct 2018 16:41:28 +0200
Subject: [PATCH] intersect_SNP.R: add restriction sequence
---
src/intersect_SNP.R | 27 +++++++++++++++++++++++++--
1 file changed, 25 insertions(+), 2 deletions(-)
diff --git a/src/intersect_SNP.R b/src/intersect_SNP.R
index 6d2649f..6cd39fb 100755
--- a/src/intersect_SNP.R
+++ b/src/intersect_SNP.R
@@ -1,5 +1,6 @@
#!/usr/bin/Rscript
library("tidyverse")
+library("seqinr")
args <- commandArgs(trailingOnly = TRUE)
snp_a <- read_delim(args[1], delim = "\t") %>%
@@ -12,8 +13,30 @@ only_b <- snp_b %>%
select(cords) %>%
setdiff(snp_a %>% select(cords)) %>%
pull(cords)
+snp <- snp_b %>%
+ filter(cords %in% only_b)
-snp_b %>%
- filter(cords %in% only_b) %>%
+fastafile <- read.fasta(file = args[3],
+ as.string = TRUE)
+
+snp$seq_list <- snp %>%
+ apply(1, FUN = function(x, fastafile, POS, CHROM){
+ begin <- as.integer(x[ POS ]) - 10
+ end <- as.integer(x[ POS ]) + 10
+ chrom <- x[ CHROM ]
+ print(paste0(chrom, ":", begin, "-", end))
+ seq_restric <- fastafile[[ chrom ]] %>%
+ c2s() %>%
+ substr(begin, end) %>%
+ s2c()
+ seq_restric[11] <- toupper(seq_restric[11])
+ return(seq_restric %>% c2s())
+ },
+ fastafile = fastafile,
+ POS = which(colnames(snp) %in% "POS"),
+ CHROM = which(colnames(snp) %in% "CHROM")
+ )
+
+snp %>%
write_csv(paste0("only_", args[2]))
--
GitLab