如何将 *.csv 文件的列从一列到多列（Python、BioPython、Pandas）？

Question

我还是个初学者，我的朋友（到目前为止还没有回答）为我提供了从 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

等等。等等，想象一下超过 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

Answer 1

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')

我们的数据框如下所示：

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

Answer 2

您拥有的脚本正在创建一个制表符分隔的输出文件，而不是逗号分隔的输出文件。如果您删除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)