Читайте нуклеотиды в FASTA без использования BioPython - PullRequest
Новая
       39

Читайте нуклеотиды в FASTA без использования BioPython

0 голосов
/ 26 декабря 2018

Мне нужно получить тот же вывод, полученный с помощью следующего кода, но без использования BioPython.Я застрял ... Кто-нибудь может мне помочь?Спасибо !!! 

from Bio import SeqIO
records = SeqIO.parse("data/assembledSeqs.fa", "fasta")
for i, seq_record in enumerate(records):
    print("Sequence %d:" % i)
    print("Number of A's: %d" % seq_record.seq.count("A"))
    print("Number of C's: %d" % seq_record.seq.count("C"))
    print("Number of G's: %d" % seq_record.seq.count("G"))
    print("Number of T's: %d" % seq_record.seq.count("T"))
    print()

Файл FASTA выглядит так: 

>chr12_9180206_+:chr12_118582391_+:a1;2 total_counts: 115 Seed: 4 K:    20 length: 79
TTGGTTTCGTGGTTTTGCAAAGTATTGGCCTCCACCGCTATGTCTGGCTGGTTTACGAGC
AGGACAGGCCGCTAAAGTG
>chr12_9180206_+:chr12_118582391_+:a2;2 total_counts: 135 Seed: 4 K: 20 length: 80
CTAACCCCCTACTTCCCAGACAGCTGCTCGTACAGTTTGGGCACATAGTCATCCCACTCG
GCCTGGTAACACGTGCCAGC
>chr1_8969882_-:chr1_568670_-:a1;113 total_counts: 7600 Seed: 225 K: 20 length: 86
CACTCATGAGCTGTCCCCACATTAGGCTTAAAAACAGATGCAATTCCCGGACGTCTAAAC
CAAACCACTTTCACCGCCACACGACC
>chr1_8969882_-:chr1_568670_-:a2;69 total_counts: 6987 Seed: 197 K: 20   length: 120
TGAACCTACGACTACACCGACTACGGCGGACTAATCTTCAACTCCTACATACTTCCCCCA
TTATTCCTAGAACCAGGCGACCTGCGACTCCTTGACGTTGACAATCGAGTAGTACTCCCG

Я пробовал это, но не работает вообще 

def count_bases (fasta_file_name):
    with open(fasta_file_name) as file_content:
        for seqs in file_content:
            if seqs.startswith('>'):
                for i, seq in enumerate('>'):
                    print("Sequence %d:" % i)
            else:
                print("Number of A's: %d" % seqs.count("A"))
                print("Number of C's: %d" % seqs.count("C"))
                print("Number of G's: %d" % seqs.count("G"))
                print("Number of T's: %d" % seqs.count("T"))
                print()
    return bases

result = count_bases('data/assembledSeqs.fa')

Ответы [ 3 ]

0 голосов
/ 26 декабря 2018

Я бы использовал для этого диктовку Counter следующим образом: 

from collections import Counter

with open('temp.fasta', 'r') as f:
    seq_name, d = '', ''
    for line in f:
        if line.startswith('>'):
            print(seq_name, d)
            seq_name = line.rstrip('\n')
            d = Counter()
        else:
            for i in line.rstrip('\n'):
                d[i] += 1

Я взял приведенный вами пример и сохранил его как temp.fasta и после запуска кода, показанного выше,вывод выглядит следующим образом: 

>chr1_8969882_-:chr1_568670_-:a2;69 total_counts: 6987 Seed: 197 K: 20   length: 120 Counter({'C': 41, 'A': 31, 'T': 28, 'G': 20})
>chr12_9180206_+:chr12_118582391_+:a1;2 total_counts: 115 Seed: 4 K:    20 length: 79 Counter({'T': 24, 'G': 24, 'C': 17, 'A': 14})
>chr12_9180206_+:chr12_118582391_+:a2;2 total_counts: 135 Seed: 4 K: 20 length: 80 Counter({'C': 30, 'T': 17, 'A': 17, 'G': 16})
>chr1_8969882_-:chr1_568670_-:a1;113 total_counts: 7600 Seed: 225 K: 20 length: 86 Counter({'C': 31, 'A': 27, 'T': 16, 'G': 12})

Я оставлю форматирование на ваше усмотрение.

0 голосов
/ 30 декабря 2018

Я чувствую, что предложение @ Kayvee о Counter из библиотеки коллекций является правильным, но пример кода, который поставляется с этим предложением, терпит неудачу по ряду причин, то есть не работает и не работаетиспользуйте Counter хорошо.Вот как я мог бы подойти к этой проблеме: 

from collections import Counter

def display(index, counter):
    print("Sequence {}:".format(index))
    for nucleotide in 'ACGT':
        print("Number of {}'s: {}".format(nucleotide, counter[nucleotide]))
    print()

with open('assembledSeqs.fa') as file:
    index = 0
    counter = Counter()

    for line in file:
        string = line.rstrip('\n')

        if string.startswith('>'):

            if index > 0:
                display(index, counter)

            index += 1
            counter.clear()
        else:
            counter.update(string)

    if counter:  # pop out the last one in the now closed file
        display(index, counter)

Этот вывод должен соответствовать тому, что вы получаете от BioPython: 

> python3 test.py
Sequence 1:
Number of A's: 14
Number of C's: 17
Number of G's: 24
Number of T's: 24

Sequence 2:
Number of A's: 17
Number of C's: 30
Number of G's: 16
Number of T's: 17

Sequence 3:
Number of A's: 27
Number of C's: 31
Number of G's: 12
Number of T's: 16

Sequence 4:
Number of A's: 31
Number of C's: 41
Number of G's: 20
Number of T's: 28

>
0 голосов
/ 26 декабря 2018

Это будет работать.Сначала мы хотим разобрать ваш файл fasta в заголовки и последовательности, затем сгладить список списков последовательностей в список строк, затем подсчитать количество каждого нуклеотида в строке, а затем вывести: 

import sys

fasta = sys.argv[1]

def fastaParser(infile):
    seqs = []
    headers = []
    with open(infile) as f:
        sequence = ""
        header = None
        for line in f:
            if line.startswith('>'):
                headers.append(line[1:-1])
                if header:
                    seqs.append([sequence])
                sequence = ""
                header = line[1:]
            else:
                sequence += line.rstrip()
        seqs.append([sequence])
    return headers, seqs

headers, seqs = fastaParser(fasta)

flat_seqs = [item for sublist in seqs for item in sublist]

def countNucs(instring):
    # will count upper and lower case sequences, if do not want lower case remove .upper()
    g = instring.upper().count('G') 
    c = instring.upper().count('C')
    a = instring.upper().count('A')
    t = instring.upper().count('T')
    return 'G = {}, C = {}, A = {}, T = {}'.format(g, c, a, t)

for header, seq in zip(headers, flat_seqs):
    print(header, countNucs(seq))

используя ваш пример fasta для запуска в командной строке: 

[me]$ python3.5 fasta.py infasta.fa

output: 

chr12_9180206_+:chr12_118582391_+:a1;2 total_counts: 115 Seed: 4 K:    20 length: 79 G = 24, C = 17, A = 14, T = 24
chr12_9180206_+:chr12_118582391_+:a2;2 total_counts: 135 Seed: 4 K: 20 length: 80 G = 16, C = 30, A = 17, T = 17
chr1_8969882_-:chr1_568670_-:a1;113 total_counts: 7600 Seed: 225 K: 20 length: 86 G = 12, C = 31, A = 27, T = 16
chr1_8969882_-:chr1_568670_-:a2;69 total_counts: 6987 Seed: 197 K: 20   length: 120 G = 20, C = 41, A = 31, T = 28
...