Entering edit mode
Hi everyone, I'm trying to obtain a GC content track for the M. tuberculosis genome to plot on IGV. Below there is the code I have used to do this. However, when calculating the GC content, R blocks without giving an output. How can I solve the problem?
library(Biostrings)
library(Gviz)
#load FASTA file
genome_seq<-readDNAStringSet("M.tuberculosis_NC000962.3.fasta")
# Calculate genome length
genome_length <- sum(width(genome_seq))
# GC content for a certain window size
window_size <- 10000
gc_content <- numeric(genome_length)
for (i in 1:(genome_length - window_size + 1)) {
window_seq <- subseq(genome_seq, start = i, end = i + window_size - 1)
gc_content[i] <- sum(unlist(strsplit(as.character(window_seq), "")) %in% c("G", "C")) / window_size * 100
}
Thank everyone for help