Basic genes codon count in R with seqinr

Hello everyone, im new to the r and i have some problems. I need to use just seqinr, so my question is ;

I need to read my genes in FASTA into the list genes so it should be like ;

   library(seqinr)
    genes <- read.fasta('geness.txt')
    length(genes)
## [1] 4120

and i can see my number of genes with the length.

but I need to apply codon_count() to each gene and store the results in the list gene_codon_count.

    codon_count <- function(gene) {
  answer <- rep(0, 64)
  names(answer) <- c("aaa", "aac", "aag", "aat", "aca", "acc", "acg", "act",
      "aga",  "agc", "agg", "agt", "ata", "atc", "atg", "att", "caa", "cac", 
      "cag", "cat", "cca", "ccc", "ccg", "cct", "cga", "cgc", "cgg",  "cgt",
      "cta", "ctc", "ctg", "ctt", "gaa", "gac", "gag", "gat",  "gca", "gcc",
      "gcg", "gct", "gga", "ggc", "ggg", "ggt", "gta",  "gtc", "gtg", "gtt",
      "taa", "tac", "tag", "tat", "tca", "tcc",  "tcg", "tct", "tga", "tgc",
      "tgg", "tgt", "tta", "ttc", "ttg",  "ttt")
  for(i in seq(from=1, to=length(gene), by=3)) {
    codon <- tolower(paste0(gene[i], gene[i+1], gene[i+2]))
    answer[codon] <- answer[codon] + 1
  }
  return(answer)
}

and the answer should be

length(genes_codon_count)
## [1] 4120
length(genes_codon_count[[1]])
## [1] 64

But i cant get any results so im missing or doing something wrong.

gene sequence genome

• 174 views • link

written 6.2 years ago by monlawski •

From your code, it seems to me that you are trying to count the codons and extract codon usage in sequences from a file. There is easy way to do this:

Example fasta file with dummy sequences:

$ cat test.fa
>a
atgcgc
>b
gcggcgg
>c
gtcgtcg

Load this in to R:

## Load Biostrings library
> library(Biostrings)

## Load fasta file
> test=readDNAStringSet("test.fa")

## Print test file
> test
DNAStringSet object of length 3:
    width seq                                  names               
[1]     6 ATGCGC                               a
[2]     7 GCGGCGG                              b
[3]     7 GTCGTCG                              c

## Print codon (tri nucleotide) summary (with one base overlapping) and remove the non-existing codons in the sequence
> trinucleotideFrequency(test)[,colSums(trinucleotideFrequency(test))!=0]
     ATG CGC CGG CGT GCG GGC GTC TCG TGC
[1,]   1   1   0   0   1   0   0   0   1
[2,]   0   0   2   0   2   1   0   0   0
[3,]   0   0   0   1   0   0   2   2   0

## Count number of codons over all
> sum(trinucleotideFrequency(test))
[1] 14

## Print codon (tri nucleotide) summary (with no base overlapping) and remove the non-existing codons in the sequence

> as.data.frame(trinucleotideFrequency(test,step = 3), row.names = names(test))[,colSums(trinucleotideFrequency(test,step = 3))!=0]
  ATG CGC GCG GTC
a   1   1   0   0
b   0   0   2   0
c   0   0   0   2

## Count number of codons
> sum(trinucleotideFrequency(test, step = 3))
[1] 6

• 0 views • link

written 6.2 years ago by cpad0112

0 answers

No answers yet.