在 Python 或 R 中连接多个 DNA 序列文本文件？

Question

我想知道如何使用 Python 或 R 连接 exon/DNA fasta 文件。

示例文件：

到目前为止，我非常喜欢将 R ape 包用于 cbind 方法，完全是因为 fill.with.gaps=TRUE 属性。当一个物种缺少外显子时，我真的需要插入缺口。

我的代码：

ex1 <- read.dna("exon1.txt", format="fasta")
ex2 <- read.dna("exon2.txt", format="fasta")
output <- cbind(ex1, ex2, fill.with.gaps=TRUE)
write.dna(output, "Output.txt", format="fasta")

示例：

exon1.txt

>sp1
AAAA 
>sp2
CCCC

exon2.txt

>sp1
AGG-G
>sp2
CTGAT
>sp3
CTTTT

输出文件：

>sp1
AAAAAGG-G
>sp2
CCCCCTGAT
>sp3
----CTTTT

到目前为止，当我有多个外显子文件时，我在尝试应用此技术时遇到了麻烦（试图找出一个循环来打开目录中所有以 .fa 结尾的文件并执行 cbind 方法），有时并非所有文件都具有长度相同的外显子 - 因此 DNAbin 停止工作。

到目前为止我有：

file_list <- list.files(pattern=".fa") 

myFunc <- function(x) {
   for (file in file_list) {
     x <- read.dna(file, format="fasta")
     out <- cbind(x, fill.with.gaps=TRUE)
     write.dna(out, "Output.txt", format="fasta")
   }
}

然而，当我运行这个并检查我的输出文本文件时，它遗漏了许多外显子，我认为这是因为并非所有文件都具有相同的外显子长度......或者我的脚本在某处失败并且我想不通:(

有什么想法吗？我也可以试试 Python.

Answer 1

我刚刚在 Python 3:

中得出这个答案

def read_fasta(fasta): #Function that reads the files
  output = {}
  for line in fasta.split("\n"):
    line = line.strip()
    if not line:
      continue
    if line.startswith(">"):
      active_sequence_name = line[1:]
      if active_sequence_name not in output:
        output[active_sequence_name] = []
      continue
    sequence = line
    output[active_sequence_name].append(sequence)
  return output

with open("exon1.txt", 'r') as file: # read exon1.txt
  file1 = read_fasta(file.read())
with open("exon2.txt", 'r') as file: # read exon2.txt
  file2 = read_fasta(file.read())

finaldict = {}                                     #Concatenate the
for i in list(file1.keys()) + list(file2.keys()):  #both files content
  if i not in file1.keys():
    file1[i] = ["-" * len(file2[i][0])]
  if i not in file2.keys():
    file2[i] = ["-" * len(file1[i][0])]
  finaldict[i] = file1[i] + file2[i]

with open("output.txt", 'w') as file:  # output that in file 
  for k, i in finaldict.items():       # named output.txt
    file.write(">{}\n{}\n".format(k, "".join(i))) #proper formatting

很难完整地评论和解释它，它可能对你没有帮助，但这总比没有好:P

我使用了 Łukasz Rogalski 对的回答中的代码。

Answer 2

如果您更喜欢使用 Linux 一个衬垫，您有

      cat exon1.txt exon2.txt > outfile

如果您只需要输出文件中的唯一记录，请使用

      awk '/^>/{f=!d[];d[]=1}f' outfile > sorted_outfile

在 Python 或 R 中连接多个 DNA 序列文本文件？

Concatenate multiple text files of DNA sequences in Python or R?

python

r

concatenation

dna-sequence

示例文件：