python - texto - dividir un archivo fasta y renombrar sobre la base de la primera línea

strip en python

Question

Tengo un gran archivo con el siguiente contenido:

nombre de archivo: input.txt

>chr1
jdlfnhl
dh,ndh
dnh.

dhjl

>chr2
dhfl
dhl
dh;l

>chr3

shgl
sgl

>chr2_random
dgld

Necesito dividir este archivo de tal manera que obtenga cuatro archivos separados como a continuación:

archivo 1: chr1.fa

>chr1
jdlfnhl
dh,ndh
dnh.

dhjl

archivo 2: chr2.fa

>chr2
dhfl
dhl
dh;l

archivo 3: chr3.fa

>chr3

shgl
sgl

archivo 4: chr2_random.fa

>chr2_random
dgld

Intenté csplit en Linux, pero no pude cambiar el nombre por el texto inmediatamente después de ">".

csplit -z input.txt ''/>/'' ''{*}''

Answer 1

Su mejor opción sería usar el programa fastaexplode de la suite exonerada :

$ fastaexplode -h
fastaexplode from exonerate version 2.2.0
Using glib version 2.30.2
Built on Jan 12 2012
Branch: unnamed branch

fastaexplode: Split a fasta file up into individual sequences
Guy St.C. Slater. [email protected]. 2000-2003.

Synopsis:
--------
fastaexplode <path>

General Options:
---------------
-h --shorthelp [FALSE] <TRUE>
   --help [FALSE] 
-v --version [FALSE] 

Sequence Input Options:
----------------------
-f --fasta [mandatory]  <*** not set ***>
-d --directory [.] 

--

Answer 2

Script ligeramente desordenado, pero debería funcionar en un archivo grande ya que solo lee una línea a la vez

Para ejecutar, haz python thescript.py input.txt (o leerá de stdin, como cat input.txt | python thescript.py )

import sys
import fileinput

in_file = False

for line in fileinput.input():
    if line.startswith(">"):
        # Close current file
        if in_file:
            f.close()

        # Make new filename
        fname = line.rstrip().partition(">")[2]
        fname = "%s.fa" % fname

        # Open new file
        f = open(fname, "w")
        in_file = True

        # Write current line
        f.write(line)

    elif in_file:
        # Write line to currently open file
        f.write(line)

    else:
        # Something went wrong, no ">chr1" found yet
        print >>sys.stderr, "Line %r encountered, but no preceeding > line found"

Answer 3

Si desea hacer algo más complicado con los archivos FASTA / FASTQ, debería considerar Biopython.

Aquí hay una publicación sobre cómo modificar y volver a escribir archivos FASTQ: http://news.open-bio.org/news/2009/09/biopython-fast-fastq/

Y otra sobre la división de archivos FASTA: http://lists.open-bio.org/pipermail/biopython/2012-July/008102.html

Answer 4

with open(''data.txt'') as f:
    lines=f.read()
    lines=lines.split(''>'')
    lines=[''>''+x for x in lines[1:]]
    for x in lines:
        file_name=x.split(''/n'')[0][1:]  #use this variable to create the new file
        fil=open(file_name+''.fa'',''w'')
        fil.write(x)
        fil.close()

Answer 5

Si específicamente quieres probar esto con Python, puedes usar este código

f2 = open("/dev/null", "r")
f = open("input.txt", "r")
for line in f:
    if ">" in line:
        f2.close()
        f2 = open(line.split(">")[1]),"w")
    else:
        f2.write(line)

f.close()

Answer 6

Como indica que está en una caja de Linux, ''awk'' parece ser la herramienta adecuada para el trabajo.

USO:
./foo.awk your_input_file

foo.awk:

#!/usr/bin/awk -f

/^>chr/ {
    OUT=substr($0,2) ".fa"
}

OUT {
    print >OUT
}

Puedes hacer eso también en una línea:

awk ''/^>chr/ {OUT=substr($0,2) ".fa"}; OUT {print >OUT}'' your_input

Answer 7

Alternativamente, BioPython podría haber sido utilizado. Instalarlo en un virtualenv es fácil:

virtualenv biopython_env
source biopython_env/bin/activate
pip install numpy
pip install biopython

Y una vez hecho esto, dividir el archivo fasta es fácil. Supongamos que tiene la ruta al archivo fasta en la variable fasta_file :

from Bio import SeqIO

parser = SeqIO.parse(fasta_file, "fasta")
for entry in parser:
   SeqIO.write(entry, "chr{}.fa".format(entry.id), "fasta")

Tenga en cuenta que esta versión de formato funciona en Python2.7, pero podría no funcionar en versiones anteriores.

En cuanto al rendimiento, utilicé esto para dividir la referencia del genoma humano del proyecto 1000 Genomes en un tiempo insignificante, pero no sé cómo funcionaría para archivos más grandes.

Answer 8

#!/usr/bin/perl-w
use strict;
use warnings;


my %hash =();
my $key = '''';
open F, "input.txt", or die $!;
while(<F>){
    chomp;
    if($_ =~ /^(>.+)/){
        $key = $1;
    }else{
       push @{$hash{$key}}, $_ ;
    }   
}

foreach(keys %hash){
    my $key1 = $_;
    my $key2 ='''';
    if($key1 =~ /^>(.+)/){
        $key2 = $1;
    }
    open MYOUTPUT, ">","$key2.fa", or die $!;
    print MYOUTPUT join("/n",$_,@{$hash{$_}}),"/n";
    close MYOUTPUT;
}