Блог

Я все еще довольно новичок, и мой друг (который пока не отвечает) предоставил мне код для загрузки геномной последовательности с Ensembl.org и записать его в файл *.csv с помощью словарей. К несчастью, файл содержит только один столбец и 89870 строк, я не уверен, как это исправить. Это облегчило бы мою работу с подсчетом, потому что это странно действует при составлении сюжетов. Я не знаю, где может быть ошибка. Вот код:

 from Bio.SeqIO.FastaIO import FastaIterator

record_ids = []
records = []

with open("equus_cds.fa") as handle:
     for record in FastaIterator(handle):
            record_ids.append(record.id)
            records.append(record)

data_cds = {}

for record in records:
    data_cds[record.id] = {'A': 0, 'G': 0, 'C': 0, 'T': 0, 'N': 0}
    for letter in str(record.seq):
        data_cds[record.id][letter]  = 1

import csv

with open('data_cds.csv', 'w') as csvfile:
    writer = csv.writer(csvfile, delimiter = "t")
    writer.writerow(['ID', 'A', 'G', 'C', 'T', 'N'])
    for key, values in data_cds.items():
        writer.writerow([key, values['A'], values['G'], values['C'], values['T'], values['N']])

with open ("data_cds.csv") as file:
    print (file.readline())
    for lines in file.readlines():
        print(lines)
 

На выходе отображается прокручивающееся оглавление, но оно немного сдвинуто:

     ID  A   G   C   T   N



ENSECAT00000046986.1    67  64  83  71  0



ENSECAT00000031957.1    81  83  75  85  0
 

и т. д. и т. Д., представьте себе более 80 тысяч строк.
Затем я хотел бы посчитать сумму всех «N» (она не всегда равна нулю), и я понятия не имею, как это сделать с этим форматом…
Заранее спасибо!

ИЗМЕНИТЬ: Я скачал последовательность отсюда: http://ftp.ensembl.org/pub/release-103/fasta/equus_caballus/cds/, расстегнул молнию:

 handle = gzip.open('file1.fa.gz')
with open('equus_cds.fa', 'wb') as out:
    for line in handle: 
        out.write(line)
 

А затем следует код, который я опубликовал. Файл *.csv всегда содержит имя определенного гена (ID — ENSECAT000… и т.д.), а затем азотистые основания (A, T, G, C), а также неизвестные основания (N). Затем весь этот файл содержит 8 тысяч строк, но только один столбец, я хотел бы, чтобы он был правильно разделен (каждая база на один столбец, если это возможно), потому что тогда было бы проще подсчитать, сколько каждой базы во всем файле (сколько Ns, если быть точным).
Причина, по которой я хочу это знать, заключается в том, что, когда я составляю сюжет, я сравниваю две последовательности, cd (кодирующие последовательности) и кДНК (комплементарную ДНК), и после вычитания N сюжет ведет себя странно, cd становится больше, чем кДНК, и это ерунда. Вот код для сюжета:

     data1 = pd.read_csv ("data_cds.csv", delimiter="t")

data1['x'] = (data1['G']   data1['C'] - data1['N']) / (data1['A']       data1['G']   data1['C']   data1['T'] - data1['N'])
data1['x'].plot.hist(bins=2000)
plt.xlim([0, 1])
plt.xlabel("cds GC percentage")
plt.title("Equus caballus", style="italic")
 

Я анализирую млекопитающих для своей диссертации, я не сталкиваюсь с этой проблемой у каждого вида, но этого все равно достаточно. Надеюсь, теперь мой вопрос более понятен.

ПРАВКА 2:

Я либо очень плохо разбираюсь в математике, либо здесь слишком поздно ночью, либо файл ведет себя странно… Как получилось, что суммы N оснований различны?

 df['N'].sum()
3504.0

df['cds_wo_N'] = df["A"] df["G"] df["C"] df["T"]-df["N"]
df['cds_wo_N'].sum()
88748562.0

df['cds_w_N'] = df["A"] df["G"] df["C"] df["T"] df["N"]
df['cds_w_N'].sum()
88755570.0

df['N_subt'] = df['cds_w_N']-df['cds_wo_N']
df['N_subt'].sum()
7008.0
 

SeqIO есть to_dict способ. Если вы используете это в сочетании с collections.Counter тем, что вы можете написать свой код более лаконично. Мы также поместим все в pandas.DataFrame файл напрямую и не будем проходить промежуточный этап написания CSV-файла.

 from collections import Counter
from Bio import SeqIO
import pandas as pd
import matplotlib.pyplot as plt

record_dict = SeqIO.to_dict(SeqIO.parse("Equus_caballus.EquCab3.0.cds.all.fa", "fasta"))
record_dict = {record_id: Counter(record_seq) for record_id, record_seq in record_dict.items()}
df = pd.DataFrame.from_dict(record_dict, orient='index')
 

Наш фрейм данных выглядит так:

Теперь мы можем легко отфильтровать только те записи, которые имеют неизвестные базы с df[df['N'].notnull()]

Или подсчитайте общее количество N оснований с помощью df['N'].sum() :

 3504
 

Теперь мы можем рассчитать процент GC

 df = df.fillna(0) # replace the NaNs with zero
df['cds GC percentage'] = (df['G']   df['C'] - df['N']) / (df['A']   df['G']   df['C']   df['T'] - df['N'])
 

df['cds GC percentage'] выглядит как:

И сюжет теперь выглядит следующим образом:

 df['cds GC percentage'].plot.hist(bins=2000)
plt.xlim([0, 1])
plt.xlabel("cds GC percentage")
plt.title("Equus caballus", style="italic");
 

Редактировать

Относительно вашего последнего обновления. Определите df['cds_wo_N'] следующим образом:

 df['cds_wo_N'] = df["A"] df["G"] df["C"] df["T"]
 

Сценарий, который у вас есть, создает выходной файл с разделителями табуляции, а не разделенный запятыми. Если вы удалите delimiter='t' параметр, по умолчанию он будет равен запятой.

Во — вторых, вы, похоже, получаете дополнительные пустые строки. Они удаляются путем добавления newline='' параметра при открытии выходного файла. Это указано в документации.

 from Bio.SeqIO.FastaIO import FastaIterator
import csv

record_ids = []
records = []

with open("equus_cds.fa") as handle:
     for record in FastaIterator(handle):
            record_ids.append(record.id)
            records.append(record)

data_cds = {}

for record in records:
    data_cds[record.id] = {'A': 0, 'G': 0, 'C': 0, 'T': 0, 'N': 0}
    for letter in str(record.seq):
        data_cds[record.id][letter]  = 1

with open('data_cds.csv', 'w', newline='') as csvfile:
    writer = csv.writer(csvfile, delimiter = "t")
    writer.writerow(['ID', 'A', 'G', 'C', 'T', 'N'])
    
    for key, values in data_cds.items():
        writer.writerow([key, values['A'], values['G'], values['C'], values['T'], values['N']])

with open("data_cds.csv") as file:
    for line in file:
        print(line)
 

Затем это должно привести к чему-то вроде:

 ID,A,G,C,T,N
ENSECAT00000046986.1,67,64,83,71,0
ENSECAT00000031957.1,81,83,75,85,0
 

Вы можете распаковать свой .gz файл с помощью Python следующим образом:

 import shutil
import gzip

with gzip.open('Equus_caballus.EquCab3.0.cds.all.fa.gz', 'rb') as f_in, 
    open('equus_cds.fa', 'wb') as f_out:
    shutil.copyfileobj(f_in, f_out)