从文件名、contig标识符和序列长度创建数据帧

Question

我试图从fasta文件中创建一个数据帧，其中包含一个头文件（contig的名称）和一个DNA序列。在我的数据帧的第一列中，我想知道文件名，在第二行中，我想知道contig名称，在第三列中，我想知道contig序列的长度（碱基对的数量-我不必计算这个-它也在contig ID中，所以我可以稍后拆分它）。你知道吗

在jupyter笔记本（嵌入bash shell）中，我尝试了以下方法：

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

我希望这会用filename1、filename2、filename3填充files=[] 标识符=[]使用>；contig\u id\u 1\u length=309，>；contig\u id\u 2\u length=400，>；contig\u id\u 3\u length=40009等。然后我可以使用split（）拆分contig id以检索contig的长度，并将所有3个系列添加到pd数据帧。你知道吗

Answer 1

所以我生成了一些虚拟数据：

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

下面的代码是一种“纯”python方法，它不需要任何其他模块（除了panda，对于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

或者，您可以使用来自^{}的^{}：

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

它们的工作原理是相同的，即构建一个list（read_list）的tuples。由于每个元组充当一个记录pandas可以很容易地将它们转换为数据帧。你知道吗