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

user_15444515 0
Viewed 166

I'm still quite a beginner and my friend (who's not answering so far) provided me a code for downloading a genomic sequence from Ensembl.org and writing it to a *.csv file using dictionaries. Unluckily, the file contains only one column and 89870 rows, I'm not sure how to fix it. It would ease my work with counting because it acts weird when doing plots. I don't know where could be a mistake. Here's the code:

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)

The output shows a scrolling table of contents but it's a bit shifted:

    ID  A   G   C   T   N



ENSECAT00000046986.1    67  64  83  71  0



ENSECAT00000031957.1    81  83  75  85  0

etc. etc., imagine over 80 thousand lines. Then I would like to count the sum of all "N's" (it's not always zero) and I have no idea how to do it with this format... Thanks in advance!

EDIT: I've downloaded the sequence from here: http://ftp.ensembl.org/pub/release-103/fasta/equus_caballus/cds/, unzipped it:

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

And then the code I've posted follows. The *.csv file always contains a name of a specific gene (ID - ENSECAT000... etc.) and then nitrogen bases (A, T, G, C) and also unknown bases (N). This whole file then has 8k lines but only one column, I would like to have it properly separated (each base to one column, if possible) because then it would be easier to count how many of each base is in the whole file (how many Ns to be specific). The reason I want to know this is when I'm making a plot, I'm comparing two sequences, cds (coding sequences) and cDNA (complementary DNA) and after subtracting N the plot acts weird, cds gets bigger than cDNA and that's nonsense. Here's the code for the plot:

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

I'm analysing mammals for my thesis, I'm not encountering this problem with every species but it's still enough. I hope my question is more understandable now.

EDIT 2:

I'm either really bad at maths or it's too late at night here or the file acts weird... How come that the sums of N bases are different?

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
biopython pandas bioinformatics csv python
2 Answers

The script you have is creating a TAB delimited output file, not a comma separated one. If you remove the delimiter='\t' parameter, it will default to a comma.

Secondly, you appear to be getting extra blank rows. These are removed by adding the newline='' parameter when opening the output file. This is specified in the 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)

This should then produce something like:

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

You can decompress your .gz file with Python as follows:

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

SeqIO has a to_dict method. If you use that in combination with collections.Counter you can write your code more succinctly. We'll also put everything in a pandas.DataFrame directly and not go through the intermediate step of writing out a CSV file.

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

Our dataframe looks like:

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
...

We can now easily filter out only the records which have unknown bases with 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

Or count the total number of N bases with df['N'].sum():

3504

We can now calculate the GC percentage

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'] looks like:

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

And the plot now looks as follows:

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

cds GC percentage plot

Edit

Regarding your latest update. Define df['cds_wo_N'] as follows:

df['cds_wo_N'] = df["A"]+df["G"]+df["C"]+df["T"]
user_500650
Related