如何将 *.csv 文件的列从一列分隔为多列 (Python、BioPython、Pandas)?

How to separate columns of *.csv file from one to more (Python, BioPython, Pandas)?

我还是个初学者,我的朋友(到目前为止还没有回答)为我提供了一个代码,用于从 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)

输出显示滚动 table 的内容,但有一点偏移:

    ID  A   G   C   T   N



ENSECAT00000046986.1    67  64  83  71  0



ENSECAT00000031957.1    81  83  75  85  0

等等等,想象一下超过 8 万行。 然后我想计算所有“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)。然后整个文件有 8k 行但只有一列,我想将它正确分隔(如果可能的话,每个碱基到一列)因为这样会更容易计算整个文件中每个碱基的数量(如何许多 Ns 是具体的)。 我想知道这一点的原因是当我制作一个情节时,我正在比较两个序列,cds(编码序列)和 cDNA(互补 DNA)并且在减去 N 之后情节表现得很奇怪,cds 变得比 cDNA 更大,这就是废话。这是情节的代码:

    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

您的脚本正在创建一个制表符分隔的输出文件,而不是逗号分隔的文件。如果删除 delimiter='\t' 参数,它将默认为逗号。

其次,您似乎得到了额外的空白行。这些通过在打开输出文件时添加 newline='' 参数来删除。这在 documentation.

中指定
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

您可以使用 Python 解压缩 .gz 文件,如下所示:

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)

SeqIO有一个to_dict method. If you use that in combination with 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')

我们的数据框看起来像:

A G C T N
ENSECAT00000046986.1 67 64 83 71 NaN
ENSECAT00000031957.1 81 83 75 85 NaN
ENSECAT00000038711.1 85 59 82 59 NaN
ENSECAT00000058645.1 74 66 82 78 NaN
ENSECAT00000058952.1 69 63 82 71 NaN
...

我们现在可以使用 df[df['N'].notnull()]

轻松过滤掉碱基未知的记录
A G C T N
ENSECAT00000016113.2 155 264 245 135 20
ENSECAT00000048238.2 274 247 166 196 20
ENSECAT00000052603.2 370 280 283 374 1000
ENSECAT00000074965.1 654 1081 545 586 20
ENSECAT00000049830.1 177 486 458 194 20
...
ENSECAT00000029115.3 94 191 167 92 20
ENSECAT00000050439.2 734 1358 1296 717 20
ENSECAT00000058713.2 728 1353 1294 715 20
ENSECAT00000046294.1 694 1362 1341 729 20
ENSECAT00000064068.1 248 501 539 330 20

或者用df['N'].sum()统计N个碱基总数:

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'] 看起来像:

cds GC percentage
ENSECAT00000046986.1 0.515789
ENSECAT00000031957.1 0.487654
ENSECAT00000038711.1 0.494737
ENSECAT00000058645.1 0.493333
ENSECAT00000058952.1 0.508772
...

现在的情节如下所示:

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"]