le traitement des fichiers multiFASTA

StackOverflow https://stackoverflow.com/questions/1789307

Question

Je suis curieux de savoir s'il y a un outil bio-informatique là capable de traiter un fichier multiFASTA me donnant comme infos nombre de séquences, la longueur, le contenu nucléotide / acide aminé, etc., et peut-être tirer automatiquement des parcelles descriptives. Aussi une solution R BioConductor ou un module BioPerl ferait, mais je n'a pas réussi à trouver quoi que ce soit.

Pouvez-vous me aider? Merci beaucoup: -)

Était-ce utile?

La solution

Certains des outils sont gaufrer une collection de petits outils qui peuvent vous aider.

Pour compter le nombre d'entrées FASTA, j'utilise:  grep -c '^>' mySequences.fasta.

Pour vous assurer qu'aucun des entrées sont en double, je vérifie que je reçois le même nombre lors de cette opération: grep '^>' mySequences.fasta | sort | uniq | wc -l

Autres conseils

Vous pouvez également être intéressé par faSize , qui est un outil de la

Il convient de noter (pour tous ceux qui tombent sur cela, comme je l'ai fait) qu'il ya une bibliothèque Python robuste spécifiquement conçu pour gérer ces tâches appelées Biopython . En quelques lignes de code, vous pouvez accéder rapidement des réponses pour toutes les questions ci-dessus. Voici quelques exemples très basiques, principalement adaptés à partir du lien. Il y a des graphiques GC%-plaque de la chaudière et des graphiques de longueur de séquence dans le tutoriel également. 

In [1]: from Bio import SeqIO

In [2]: allSeqs = [seq_record for seq_record in SeqIO.parse('/home/kevin/stack/ls_orchid.fasta', """fasta""")]

In [3]: allSeqs[0]
Out[3]: SeqRecord(seq=Seq('CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGATGAGACCGTGG...CGC', SingleLetterAlphabet()), id='gi|2765658|emb|Z78533.1|CIZ78533', name='gi|2765658|emb|Z78533.1|CIZ78533', description='gi|2765658|emb|Z78533.1|CIZ78533 C.irapeanum 5.8S rRNA gene and ITS1 and ITS2 DNA', dbxrefs=[])

In [4]: len(allSeqs) #number of unique sequences in the file
Out[4]: 94

In [5]: len(allSeqs[0].seq) # call len() on each SeqRecord.seq object
Out[5]: 740

In [6]: A_count = allSeqs[0].seq.count('A')
        C_count = allSeqs[0].seq.count('C')
        G_count = allSeqs[0].seq.count('G')
        T_count = allSeqs[0].seq.count('T')

        ​print A_count # number of A's

        144

In [7]: allSeqs[0].seq.count("AUG") # or count how many start codons
Out[7]: 0

In [8]: allSeqs[0].seq.translate() # translate DNA -> Amino Acid
Out[8]: Seq('RNKVSVGEPAEGSLMRPWNKRSSESGGPVYSAHRGHCSRGDPDLLLGRLGSVHG...*VY', HasStopCodon(ExtendedIUPACProtein(), '*'))
Licencié sous: CC-BY-SA avec attribution
Non affilié à StackOverflow
scroll top