-
Notifications
You must be signed in to change notification settings - Fork 0
/
evolutionAlignment.R
30 lines (30 loc) · 1.4 KB
/
evolutionAlignment.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
#!/usr/bin/R
# Universitat Potsdam
# Author Gaurav Sablok
# date: 2024-3-17
evolutionAlignment <- function(fasta_file, path){
library(ape)
library(Biostrings)
library(seqinr)
library(msa)
# a evolutionaryAlignment function that calculates the evolutionary
# rates across the fasta file, given a fasta file, it will automatically
# align the sequences and then will estimate the evolutionary rates on
# the clean alignments.To use this
# evolutionAlignment("/Users/gauravsablok/Desktop/CodeRelease/fasta_sample_datasets/read_check.fasta",
# path = "/Users/gauravsablok")
# it will also plot the evolutionary rates and it uses two methods ka/ks and dn/ds
setwd(path)
fasta_file <- readDNAStringSet(fasta_file)
fasta_alignment <- msa(fasta_file)
fasta_convert <- msaConvert(fasta_alignment, type = "ape::DNAbin")
write.FASTA(fasta_convert, file = "fasta_convert.fasta")
fasta_re_check <- read.alignment(paste(getwd(),
"fasta_convert.fasta", sep = "/"), format = "fasta")
kaks_fasta <- kaks(fasta_re_check, rmgap = TRUE)
dnds_fasta <- dnds(fasta_re_check)
ka_histogram <- hist(kaks_fasta$ka, main = paste("Histogram of ka values"))
ka_histogram <- hist(kaks_fasta$ka, main = paste("Histogram of ka values"))
save.image(file = "ka_histogram", safe = TRUE)
save.image(file = "ks_histogram", safe = TRUE)
}