Объединить несколько текстовых файлов последовательностей ДНК в Python или R?

Мне было интересно, как объединить файлы exon/DNA fasta с помощью Python или R.

Примеры файлов:

До сих пор мне очень нравилось использовать пакет 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

До сих пор у меня возникли проблемы с попыткой применить эту технику, когда у меня есть несколько файлов exon (пытаясь найти цикл для открытия и выполнения метода cbind для всех файлов, заканчивающихся на .fa в каталоге), и иногда не все файлы имеют экзоны все одинаковые по длине, поэтому 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.


python r concatenation dna-sequence
person CuriousDude    schedule 09.06.2017    source источник
comment
Можете ли вы объяснить больше о том, что он должен делать?   -  person JakeD    schedule 09.06.2017
comment
В двух словах, он объединяет 2 строки ДНК файлов Fasta.   -  person Iván C.    schedule 10.06.2017
comment
Он объединяет два файла fasta, поскольку сейчас у меня есть гены, разбитые на файлы exon, но мне нужно их объединить для запуска определенной программы.   -  person CuriousDude    schedule 11.06.2017

Ответы (2)


arrow_upward
0
arrow_downward

Я только что дал этот ответ в 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

Я использовал код Лукаша Рогальского из ответа на Чтение формата файла fasta в Python dict .

person Iván C.    schedule 09.06.2017
comment
Благодарю вас! Я попробую этот код и немного адаптирую его, чтобы объединить все файлы в определенном каталоге, поскольку было бы утомительно использовать с файлами open for 20+ exon. Я также знаю только Python 2.7, но еще раз спасибо за помощь! Я ценю это! - person CuriousDude; 11.06.2017
comment
Я думал, что это будет работать в Python 2, но я не проверял: P - person Iván C.; 11.06.2017
comment
Он работал в Python 2.7. Еще раз спасибо. Ваши комментарии и комментарии Лакшми помогли мне придумать код для оптимизации моей работы. Ваше здоровье! - person CuriousDude; 16.06.2017

arrow_upward
1
arrow_downward

Если вы предпочитаете использовать один лайнер Linux, у вас есть

      cat exon1.txt exon2.txt > outfile

если вы хотите использовать только уникальные записи из внешнего файла

      awk '/^>/{f=!d[$1];d[$1]=1}f' outfile > sorted_outfile
person Lakshmi KrishnaKumaar    schedule 12.06.2017
comment
Я пробовал это, но по какой-то причине это испортило форматирование Fasta, и я действительно не знаю, почему. Мой файл заканчивается ›именами, которые также соединяются случайными строками.... - person CuriousDude; 14.06.2017
comment
Предыдущий комментарий сортирует данные как таковые без сохранения формата fasta, извините за неправильную команду. Я отредактировал вторую команду, которая может решить вашу задачу @DNAngel - person Lakshmi KrishnaKumaar; 15.06.2017
comment
Благодарю вас! Это очень помогло :) - person CuriousDude; 16.06.2017