1. Sobre

Extrair sequências em um arquivo FASTA acaba sendo necessário mesmo para aqueles que não são especialistas em Bioinformática. Imagina que você tem o genoma viral e deseja extrair sequencias montadas a partir de uma lista de proteínas/genes de interesse. E, para isso, vamos utilizar o R!

Além disso, vamos verificar algumas opções de manipulações de arquivos FASTA. Um arquivo FASTA é composto por um sinmbolo '>' seguido de uma descrição, e, logo abaixo, a sua sequência, podendo ela ser em DNA ou AA. Por exemplo:

>Seq1
ATATGCAGCGATTCTACTCAGAGCGCAACTCATCACAGCAGCAC
>Seq2
ATCTATATCGGACATCATCTTACTACATCATTCAGCGAGGACTTATCAGCCAGAGC
>Seq3
ATCAGCGACTACTATCACGGCAGCAGCATCATCTCTCAGACGC
>Seq4
ATCTACGAGCAGCGACATTA

2. Instalação

Para tais análises, vamos utilizar o pacote seqinr no R. Você pode instalá-lo diretamente no RStudio.

#install.packages("seqinr")
library("seqinr")

3. Manipulações

Vamos a algumas manipulações e, no final, como extrair sequências de um FASTA. Estou supondo que você quer carregar um arquivo chamado sequencias_genoma.fasta, e ele está no seu diretório atual (verificar com o getwd()).

#carrega o fasta
yeast <- read.fasta(file = "sequencias_genoma.fasta")

Você pode buscar pelo nome indexado da sequência. A sequencia retornada vai estar separa por nucleotídeo, e você pode filtrá-la. Segue abaixo:

#primeira sequencia fasta pelo indexador ou pelo indice 
yeastseq <- getSequence(yeast$Seq1)
yeastseq <- getSequence(yeast['Seq1'])
yeastseq <- yeast[1](/lmigueel/Bioinformatica/wiki/1)

Caso você queira retornar a sequência em uma string diretamente, use:

#Como string
yeastseq2 <- getSequence(yeast$Seq1,as.string = TRUE)
yeastseq2

Agora vamos a algumas manipulações, e vou usar a sequência quebrada por nucletídeo:

#tamanho
length(yeastseq)

#frequencia
table(yeastseq)

#conteudo GC (vindo de yeastseq <- yeast[1](/lmigueel/Bioinformatica/wiki/1))
GC(yeastseq)

#numero de dinucleotideos
count(yeastseq,2)

# subset
yeastseq[c(1:100)]

#procurar o nome das sequencias
getName(yeast)

Agora vamos extrair sequências de um FASTA. Para isso vou criar um arquivo chamado list.txt que contém os IDs de interesse que quero extrair. O header também existe e deixarei como genes. Tomando como base o FASTA na Seção 1, o arquivo list.txt será:

genes
Seq2
Seq4

Here we go:

## Subset de sequencias por uma lista
myfasta <- read.fasta(file = "sequencias_genoma.fasta", seqtype = "DNA",as.string = TRUE, set.attributes = FALSE)
subsetlist <- read.table("list.txt", header=TRUE)

#faz a interseccao entre os nomes no FASTA e sua lista
my_fasta_sub <- myfasta[names(myfasta) %in% subsetlist$genes]

#salva arquivo
write.fasta(sequences = my_fasta_sub, names = names(my_fasta_sub), nbchar = 80, file.out = "subset_list.fasta")