обработка файлов multiFASTA

Question

Мне было любопытно узнать, существует ли какой-либо инструмент биоинформатики, способный обрабатывать файл multiFASTA, предоставляющий мне такую ​​информацию, как количество последовательностей, длина, содержание нуклеотидов/аминокислот и т. д.и, возможно, автоматически рисовать описательные сюжеты.Также подошло бы решение R BIOconductor или модуль BioPerl, но мне ничего не удалось найти.

Вы можете помочь мне?Большое спасибо :-)

Answer 1

Некоторые инструменты для тиснения представляют собой набор небольших инструментов, которые могут вам помочь.

• seqstats возвращает длину последовательности
• pepstats должно давать вам содержание аминокислот и т.д.Некоторые инструменты также предлагают функции построения графиков.Очень удобно.http://emboss.sourceforge.net/apps/release/5.0/emboss/apps/groups.html

Чтобы подсчитать количество записей fasta, я использую: grep -c '^>' mySequences.fasta.

Чтобы убедиться, что ни одна из записей не дублируется, я проверяю, что при выполнении этого я получаю один и тот же номер: grep '^>' mySequences.fasta | sort | uniq | wc -l

Answer 2

Вас также может заинтересовать Фасовать, который является инструментом из Исходное Дерево Кента, хотя это требует немного больше усилий (вы должны выполнить загрузку и компиляцию), чем просто использование grep...вот некоторый пример вывода:

me@my-lab ~/data $ time faSize myfile.fna
215400419 bases (104761 N's 215295658 real 215295658 upper 0 lower) in 731620 sequences in 1 files
Total size: mean 294.4 sd 138.5 min 30 (F5854LK02GG895) max 1623 (F5854LK01AHBEH) median 307
N count: mean 0.1 sd 0.4
U count: mean 294.3 sd 138.5
L count: mean 0.0 sd 0.0
%0.00 masked total, %0.00 masked real

real    0m3.710s
user    0m3.541s
sys     0m0.164s

Answer 3

Следует отметить (для тех, кто наткнулся на это, как я только что сделал), что существует надежная библиотека Python, специально разработанная для решения этих задач, называемая Биопитон.С помощью нескольких строк кода вы можете быстро получить ответы на все вышеперечисленные вопросы.Вот несколько очень простых примеров, в основном адаптированных по ссылке.В руководстве также есть стандартные графики GC% и графики длины последовательности.

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(), '*'))