filtrar o comprimento

Question

Eu tenho um arquivo que contém 3940 seqüências e se parece com

>TCONS_00000066 gene=XLOC_000030
CCGCCGGCTGCTGCGCGCACCGACTTGTCACCACCCCAGCACGTCCTCCACGTATACAAG
CGCTACGGTCCACCGCGGCAGCGTCGACGTCCTTGTCCGCAAACATGGTGGTGGCAGCTT
CCTCATCGAGCAGCAGCAACTCATCCTCGAGGGGAAGGGCCCAGAGCTTCTAATCCTACA
CGGCAACAACACTTTATACTTGTGTATAATTTCTCTTCGTTTCTGAGTTCATGGCTATCT
TTGTCTCTCTTATCTTCTCCCTTTTGCTATCTCTATATTTGTGATTGCCATGGAAATACA
>TCONS_00006042 gene=XLOC_003448
GCCACTAGCCAGCCCAGCCAGGGGAAGGGGAGGAGCTGCAAGCCCAACCCCCTGCTCAAC
CCTAAATTGCTTCCGCCGATCGGTGAGAGCTCCGATGCCTTCTTCTTCTTCTTCTTCCTC
CCCCTCTACCTGTTCCTTCTCCGAGATAACTGCAACATTTTCAGCACTTTTTCTGGCCAT
TCTCAAGTCCCCAGCCCAGGGACTAGAGTGTTACTATGGCTAGAGCAAATGAGATGGTCA
GGGCAGACTCAAGGATGATGGTTGTCTTTAGTGCCCTGGCATCTAAATCAGGGCCACTGA
CATTTGAAGACTCGCTCAGATTTGTCAAGAAAGTGAAGGCTTGTAACTACATGTTGTATT

Eu quero as sequências com mais de 200 caracteres em outro arquivo

    

Answer 1

Com awk , você precisa definir > como um separador de registro primeiro:

awk 'BEGIN{RS=">";ORS=""}length($0)>200{print ">"$0}' input > output

Outra opção com pcregrep :

pcregrep -M '^>[^>]{201,}' input > output

Ou para contar apenas sequências de DNA, não caracteres no cabeçalho:

pcregrep -M '^>[^>]*\n[^>]{201,}' input > output

    

Answer 2

Você pode usar sed, awk ou grep.

awk 'length($0)>200' file > newfile

OR

grep '^.\{201\}' file > newfile

    

Answer 3

Python ( split.py ):

import sys

# call with the file as parameter

base = 0
line = ''
with open(sys.argv[-1]) as fp:
    with open('shorter', 'w') as fps:
        with open('longer', 'w') as fpl:
            for x in fp:
                if line and x.startswith('>'):
                    print len(line), base
                    if (len(line) - base) >= 200:
                        fpl.write(line)
                    else:
                        fps.write(line)
                    line = x
                    base = len(x)  # lenght of the ">..." line
                    continue
                if x.startswith('>'):  # very first one
                    base = len(x)
                line += x
            if line:
                if len(line) >= 200:
                    fpl.write(line)
                else:
                    fps.write(line)
                line = ""

chame com python split.py inputfile e, em seguida, mv shorter inputfile (depois de verificar se os arquivos estão OK)

    

Answer 4

Você pode fazer:

f=0
for arg in -e -v
do  f=$((f+1))
    tr '>\n' '\n>' <infile |
    grep "$arg" '>.\{200\}'|
    tr '>\n' '\n>' >"./outfile$f"
done

Que irá gravar todas as sequências de 200 caracteres ou mais para outfile1 e todas as mais curtas para outfile2.

Ele traduz todos os > para \n ewlines e vice-versa - assim, seu primeiro exemplo acima se torna:

TCONS_00000066 gene=XLOC_000030>CCGCCGGCTGCTGCGCGCACCGACTTGTCACCACCCCAGCACGTCCTCCACGTATACAAG>CGCTACGGTCCACCGCGGCAGCGTCGACGTCCTTGTCCGCAAACATGGTGGTGGCAGCTT>CCTCATCGAGCAGCAGCAACTCATCCTCGAGGGGAAGGGCCCAGAGCTTCTAATCCTACA>CGGCAACAACACTTTATACTTGTGTATAATTTCTCTTCGTTTCTGAGTTCATGGCTATCT>TTGTCTCTCTTATCTTCTCCCTTTTGCTATCTCTATATTTGTGATTGCCATGGAAATACA>

Em seguida, corresponde a (ou não) 200 caracteres do primeiro > para cada linha de entrada. Para a primeira iteração do loop for - o que corresponde a linhas com 200 caracteres ou mais - ele inverte a conversão para as linhas que grep correspondem e grava os resultados para outfile1. A segunda iteração obtém sequências mais curtas e as grava no outfile2.

Você deve observar que essa contagem incluirá as novas linhas que estavam nos dados.

Aqui está outra maneira que não tem esse problema:

sed -n '/^>/!{H;$!d
};   x;s/\n[ACGT]\{20\}/&/4p
'    <infile >outfile

    

Answer 5

cat file | while read -r line; do
  if [ ${#line} -gt 200 ]; then
    echo "${line}"
  fi
done

EDIT A pergunta foi atualizada: não o comprimento de uma linha, mas um conjunto de linhas é desejado.

No script a seguir, eu echo > TCONS, caso contrário, o script pulará o último hit.

multiline=""
(cat input; echo ">TCONS string for last token") | while read line; do
        if [[ "$(echo "${line}"| cut -c1-6)" = ">TCONS" ]]; then
                if [ ${#multiline} -gt 200 ]; then
                        echo "${multiline}"
                fi
                multiline=""
        else
                multiline="${multiline}${line}"
        fi
done

    

Answer 6

Já dei a você meus scripts FastaToTbl e TblToFasta em outra resposta . Usando eles, você poderia fazer

FastaToTbl sequences.fa | awk 'length($2)>200' | TblToFasta > newfile.fa

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
  }
}