CODEWITHSUNDEEP
pythonbioinformaticsfasta
MeggyClay
MeggyClay

Reputation: 21

Dataframe creation from filename, contig identifier, and length of sequence

I am attempting to create a dataframe from fasta files which contain a header (the name of the contig) and a DNA sequence. In the first column of my dataframe I would like to have the name of the file, on the second line I would like to have the contig name, and in the third column I would like to have the length of the contig sequence (number of basepairs - I don't have to count this- it's also in the contig ID so I could split that later).

Within the jupyter notebook (embedded in a bash shell) I have tried the following: 

files = []
identifiers = []
# r=root, d=directories, f = files
for r, d, f in os.walk(path):
    for file in f:
        if '.fasta' in file:
            files.append(os.path.join(file)) #this grabs my file names and appends them to files - works
            open(file, "r")
            for line in file:
                identifiers.append(line) # this would grab the identifier - found on the first line of the file

I would expect this to fill files=[] with filename1, filename2, filename3 and identifiers=[] with >contig_id_1_length=309, >contig_id_2_length=400, >contig_id_3_length=40009 etc. Then I could split the contig ids with split() to retrieve the length of the contigs and add all 3 series to a pd dataframe.

Upvotes: 2

Views: 142

Answers (1)

Alex
Alex

Reputation: 7045

So I have generated some dummy data:

f1.fasta

>ctg_1_length=147
TCGTGGTCACCGATCGAAGATCCAATATCCGGAGATCGTCTACCTGTATGTAGTAAGCGCAAGGCCCGTTTACTGCGTCACCCTAGCAGAACGCCGACCAGGTCTCCTATAGTCACCGGCCTCGCACCTTTAAGTATGTATAGACGG
>ctg_2_length=141
GCTTGGGTGGGAACGGCTCGTGGCGGAGTACCCGAGAGTGGTTTCGGTATCTGGTGTCGTGCCAGGTTTAATTGAAAATTCAAGATTTTAAGTATCGCTTCAGATAGATTACTTACTGCGAGTGCCTTGTCACAGGGCGGG
>ctg_3_length=124
CCTTCGACCATGGATATCCTAACTCAGCCCCAGCCAGCTAACTCTGGACCAACCGAGAGCGTCTTTCTTTGATGTAACTAAGCTGGCGTTGGGCCCCCCGGTGTTCTAACGTATCTGAAGCCAA
>ctg_4_length=124
CGCGAACTTATCTTGTTATCGAAGATAGCTGTAGGAACTCGGCCAGCCCGACTATTTCGTTCGCCGCTTTCCCCTGGCTCTAGATGCAGTCCACAGATTCTTCTCAGGTGATGCGAGGAACAGG
>ctg_5_length=137
CCAACCCCTGCTCTAGGCTTACCGCCAAGCTACTCAATGGTTCGGTCGATGCAGAACGTATTACTATGTTCTCGACTCTCTGAAACCGCTGTCTACGAGGCAAGCCCCAAAATAGATGGAGGGGCCTCGCCTGTGGG

f2.fasta

>ctg_1_length=106
TCGATATTGGTTAAGGCGCGCAGCAATTTGGGAGTTGACGCACAACGTTCGGATGCGAGAGTGAGCATACGGTAGAGCCGAACCCACAATGGGTAACCGAACGACA
>ctg_2_length=60
CTACGATCTGAAATCCACTTCACGTGATCCGCGAGATGGGTTATTCGGTTTTTAGAACAT
>ctg_3_length=145
ACACTTATATCCACGATTGAGTGGCTCATCGGTGTGACACTCTGACGTCGTTTGAATACCTGCCCGGACAGGGTTTTCGTCAAACTCCCCGCGACGGTTCGTAACTGTCTGTACCCGTCGGCTGGACGAAGTTTAGATATAAAAC
>ctg_4_length=88
GAGCCGCTACATTACTTAATAACTTACAAAGGGCGAAGTCACATATTTCGTAAGAAGCATTCCTCGTCAGAATCCATTCCAAACCCCA
>ctg_5_length=87
CTACGCTAAGCTGCGGTACGACGGGGATATTACACGTACTAATCCATACCAACTAAATGGCATGTTGTTGAAGATAGCACTTTGAGG

The following code is a "pure" python approach, it doesn't require any other modules (except pandas, for the DataFrame):

import pandas as pd
from pathlib import Path

files = [x for x in Path().iterdir() if x.suffix == ".fasta"]
# [PosixPath('f1.fasta'), PosixPath('f2.fasta')]
read_list = []
for file in files:
    with file.open("r") as handle:
        for line in handle:
            if line.startswith(">"):
                line = line.strip()
                read_list.append((file.name,  # Change to file.resolve() for the absolute path
                                  *line[1:].split("=")
                                ))

df = pd.DataFrame(read_list, columns=["file", "ctg", "len"])

#        file           ctg  len
# 0  f1.fasta  ctg_1_length  147
# 1  f1.fasta  ctg_2_length  141
# 2  f1.fasta  ctg_3_length  124
# 3  f1.fasta  ctg_4_length  124
# 4  f1.fasta  ctg_5_length  137
# 5  f2.fasta  ctg_1_length  106
# 6  f2.fasta  ctg_2_length   60
# 7  f2.fasta  ctg_3_length  145
# 8  f2.fasta  ctg_4_length   88
# 9  f2.fasta  ctg_5_length   87

Alternatively, you could use SeqIO from biopython:

import pandas as pd
from pathlib import Path
from Bio import SeqIO 

files = [x for x in Path().iterdir() if x.suffix == ".fasta"]
read_list = []
for file in files:
    with file.open("r") as handle:
        for record in SeqIO.parse(handle, "fasta"):
            read_list.append((file.name, record.id, len(record.seq)))

df = pd.DataFrame(read_list, columns=["file", "ctg", "len"])

#        file               ctg  len
# 0  f1.fasta  ctg_1_length=147  147
# 1  f1.fasta  ctg_2_length=141  141
# 2  f1.fasta  ctg_3_length=124  124
# 3  f1.fasta  ctg_4_length=124  124
# 4  f1.fasta  ctg_5_length=137  137
# 5  f2.fasta  ctg_1_length=106  106
# 6  f2.fasta  ctg_2_length=60    60
# 7  f2.fasta  ctg_3_length=145  145
# 8  f2.fasta  ctg_4_length=88    88
# 9  f2.fasta  ctg_5_length=87    87

These both work on the same principle of building a list (read_list) of tuples. As each tuple acts as a record pandas can turn them into a DataFrame very easily.

Upvotes: 1

Related Questions