Extraindo subconjunto do arquivo fasta

Question

Eu tenho um arquivo fasta que se parece com isso:

>chr1
ACGGTGTAGTCG
>chr2
ACGTGTATAGCT
>chrUn
ACGTGGATATTT
>chr21
ACGTTGATGAAA
>chrX
GTACGGGGGTGG
>chrUn5
TGATAGCTGTTG

Eu só quero extrair chr1 , chr2 , chr21 , chrX junto com suas seqüências. Então a saída que eu quero é:

>chr1
ACGGTGTAGTCG
>chr2
ACGTGTATAGCT
>chr21
ACGTTGATGAAA
>chrX
GTACGGGGGTGG

Como posso fazer isso na linha de comando do Unix?

    

Answer 1

com sed :

sed -n '/^>chr1$\|^>chr2$\|^>chr21$\|^>chrX$/{p;n;p}' file

• -n suprime a saída automática.
• /.../ da expressão regular para corresponder a >chr1 , >chr2 , >chr21 ou >chrX .
• {p;n;p} se a expressão corresponder, imprima a linha, leia a próxima linha de entrada para o espaço padrão e imprima essa linha também.

Se for awk , é quase o mesmo mecanismo:

awk '/^>chr1$|^>chr2$|^>chr21$|^>chrX$/{print;getline;print;}' file

    

Answer 2

Para o exemplo simples que você mostra, onde todas as sequências cabem em uma única linha, você pode usar apenas grep (se o seu grep não suportar a opção --no-group-separator , passe a saída por grep -v -- '--' ) :

$ grep -wEA1 --no-group-separator 'chr1|chr2|chr21|chrX' file.fa 
>chr1
ACGGTGTAGTCG
>chr2
ACGTGTATAGCT
>chr21
ACGTTGATGAAA
>chrX
GTACGGGGGTGG

Supondo que você tenha sequências de várias linhas (o que parece provável se você estiver lidando com cromossomos), será necessário concatená-las em uma linha primeiro. Isso complica consideravelmente as coisas. Você poderia usar um awk one-liner:

$ awk -vRS=">" 'BEGIN{t["chr1"]=1;t["chr2"]=1;t["chr21"]=1;t["chrX"]=1}
                {if($1 in t){printf ">%s",$0}}' file.fa
>chr1
ACGGTGTAGTCG
>chr2
ACGTGTATAGCT
>chr21
ACGTTGATGAAA
>chrX
GTACGGGGGTGG

O script acima define o separador de registro como > ( vRS=">" ). Isso significa que "linhas" são definidas por >~ e não \n . Em seguida, o bloco BEGIN configura uma matriz em que cada um dos IDs de sequência de destino é uma chave. O resto simplesmente verifica cada "linha" (sequência) e, se o primeiro campo estiver na matriz ( $i in t ), ele imprime a "linha" atual ( $0 ) precedida por > .

Se você estiver fazendo esse tipo de coisa com frequência, escrever o array rapidamente se tornará irritante. Pessoalmente, eu uso os dois scripts abaixo, que eu herdei de um antigo colega de laboratório, para converter o FASTA para o formato tbl ( sequence_name<TAB>sequence , uma sequência por linha) e de volta:

• FastaToTbl

  #!/usr/bin/awk -f
  {
          if (substr($1,1,1)==">")
          if (NR>1)
                      printf "\n%s\t", substr($0,2,length($0)-1)
          else 
              printf "%s\t", substr($0,2,length($0)-1)
          else 
                  printf "%s", $0
  }END{printf "\n"}
  

• TblToFasta

  #! /usr/bin/awk -f
  {
    sequence=$NF
  
    ls = length(sequence)
    is = 1
    fld  = 1
    while (fld < NF)
    {
       if (fld == 1){printf ">"}
       printf "%s " , $fld
       if (fld == NF-1){
          printf "\n"
        }
        fld = fld+1
    }
    while (is <= ls){
      printf "%s\n", substr(sequence,is,60)
      is=is+60
    }
  }
  

Se você salvá-los no seu $PATH e torná-los executáveis, você pode simplesmente grep para suas sequências de destino (e isso funcionará para seqüências de várias linhas, ao contrário do que foi dito acima):

$ FastaToTbl file.fa | grep -wE 'chr1|chr2|chr21|chrX' | TblToFasta
>chr1 
ACGGTGTAGTCG
>chr2 
ACGTGTATAGCT
>chr21 
ACGTTGATGAAA
>chrX 
GTACGGGGGTGG

Isso é muito mais fácil de estender, pois você pode transmitir grep de um arquivo de destinos de pesquisa:

$ cat ids.txt 
chr1
chr2
chr21
chrX

$ FastaToTbl file.fa | grep -wFf ids.txt | TblToFasta
>chr1 
ACGGTGTAGTCG
>chr2 
ACGTGTATAGCT
>chr21 
ACGTTGATGAAA
>chrX 
GTACGGGGGTGG

Finalmente, se você estiver trabalhando com grandes seqüências, você pode considerar uma das várias ferramentas que podem fazer esse tipo de coisa para você. Eles serão muito mais rápidos e eficientes a longo prazo. Por exemplo, fastafetch do conjunto exonerate . Está disponível nos repositórios da maioria das distribuições Linux. Nos sistemas baseados em Debian, você pode instalá-lo com sudo apt-get install exonerate , por exemplo. Depois de instalar, você pode fazer:

## Create the index    
fastaindex -f file.fa -i file.in
## Loop and retrieve your sequences
for seq in chr1 chr2 chr21 chrX; do
     fastafetch -f file.fa -i file.in -q "$seq" 
done 
>chr1
ACGGTGTAGTCG
>chr2
ACGTGTATAGCT
>chr21
ACGTTGATGAAA
>chrX
GTACGGGGGTGG

Como alternativa, você pode usar meu próprio retrieveseqs.pl , que tem alguns outros detalhes funções:

$ retrieveseqs.pl -h

    retrieveseqs.pl will take one or more lists of ids and extract their sequences from 
    multi FASTA file

    USAGE : retrieveseqs.pl [-viofsn] <FASTA sequence file> <desired IDs, one per line>

    -v : verbose output, print a progress indicator (a "." for every 1000 sequences processed)
    -V : as above but a "!" for every desired sequence found.
    -f : fast, takes first characters of name "(/^([^\s]*)/)" given until the first space as the search string
         make SURE that those chars are UNIQUE.
    -i : use when the ids in the id file are EXACTLY identical
         to those in the FASTA file
    -h : Show this help and exit.
    -o : will create one fasta file for each of the id files
    -s : will create one fasta file per id
    -n : means that the last arguments (after the sequence file)
         passed are a QUOTED list of the names desired.
    -u : assumes uniprot format ids (separated by |)

No seu caso, você faria:

$ retrieveseqs.pl -fn file.fa "chr1 chr2 chr21 chrX"
[7 (4/4 found]
>chr1
ACGGTGTAGTCG
>chr2
ACGTGTATAGCT
>chr21
ACGTTGATGAAA
>chrX
GTACGGGGGTGG

Note que isto foi algo que escrevi para o meu próprio trabalho e não está muito bem documentado nem suportado. Ainda assim, eu tenho usado isso felizmente por anos.