#####
# Gaëlle Lelandais <gaelle.lelandais@universite-paris-saclay.fr>
#
#
# Ce contenu est sous Licence CC BY-SA 4.0
# https://creativecommons.org/licenses/by-sa/4.0/
######


# Chargement des librairies
library(seqinr)

# -- Sequences des chromosomes
genomeSeq = read.fasta(file = "sequences_all_genome.fasta")

# Fichier BED avec les coordonnées des pics
# pour lesquels les séquences doivent être récupérées.
bedFile = read.table("res_peaks.bed")

listSeq   = NULL
listNames = NULL

# Pour chaque pics
for(i in 1:nrow(bedFile)){
  
  # nom du chromose
  chrom = bedFile[i,1]
  # récupération de la séquence
  seq  = paste0(getSequence(genomeSeq)[[which(getAnnot(genomeSeq) == paste0(">", chrom))]][bedFile[i,2]:bedFile[i,3]], collapse = "")
  # création d'un nom pour le fichier FASTA
  name = paste(as.character(bedFile[i,4]), chrom, "from:", bedFile[i,2], "to:", bedFile[i,3])  
  
  # Ajout aux listes
  listSeq = c(listSeq, list(seq))
  listNames = c(listNames, list(name))
}

# Ecriture du fichier FASTA
write.fasta(listSeq, listNames, "seq_peaks.fasta", nbchar = 20)

####
# Et voilà :) 
####