Estou tentando analisar um arquivo grande em que cada duas linhas consecutivas têm o mesmo tamanho (o texto é completamente diferente). Eu procurei e meu primeiro post aqui. Eu encontrei um script e tentei modificá-lo, mas sem alegria. Arquivo é um arquivo de saída de seqüenciamento. Eu já analisei a sequência e as pontuações de qualidade, então o arquivo é assim:
CCTCGNAACCCAAAAACTTTGATTTCTNATAAGGTGCCAGCGGAGTCCTAAAAGCAACATCCGCTGATCCCTGGT
AAAAA#EEEEEEEEEEEEEEEEEEEEE#EEEEEEEEAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE
CCCCANCCAAACTCCCCACCTGACAATNTCCTCCGCCCGGATCGACCCGCCGAAGCGAGTCTTGGGTCTAAA
AAAAA#EEEEEEEEEEEAEEEEEEEEE#EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE
ATCGTNTATGGTTGAGACTAGGACGGTNTCTGATCGTCTTCGAGCCCCCAACTTTCGTTCTTGATTAATGAAAAC
AAAAA#EEEEEEEEEEEEEEEEEEEEE#EEEEEEEEEEEEAEEEEEEAEEEAEEEEEEEEEEEEEEEEEEEEEEE
CCCACNTGGAGCTCTCGATTCCGTGGGNTGGCTCAACAAAGCAGCCACCCCGTCCTACCTATTTAAAGTTTG
AAAAA#EEEEEEEEEEEEEEEEEEEEE#EEEEEEEEEEEEEEEEEEAEEEEEEEEEEEEEEEEEEEEEEEEE
GCATCNTTTATGGTTGAGACTAGGACGNTATCTGATCGTCTTCGAGCCCCCAACTTTCGTTCTTGATTAATGAA
6AAAA#EEEEEAAAEEEEEEAEEAEEE#EEEEEEEAEAEEEEAEEAAA/EAEEEEAEEAEEAEEAEAAEEEEEE
O problema: Em algum lugar existe um par corrompido de linhas de tal forma que cada seqüência base não tem uma pontuação correspondente, ou seja, os comprimentos de cada par de duas linhas devem ser iguais, como posso analisar o par que está incorreto ? O arquivo tem 100 milhões de linhas.
Eu tentei este código chamado parser.sh:
{ curr = $0 }
(NR%2)==0 {
currLgth = length(curr)
prevLgth = length(prev)
maxLgth = (currLgth > prevLgth ? currLgth : prevLgth)
if (prevLgth==currLgth) {
print ""
print prevLgth
print currLgth
for (i=1; i<=maxLgth; i++) {
}
}
}
{ prev = curr }
e executaria awk -f parser.sh filename
mas isso imprimiu todos os comprimentos de linhas, embora eu estivesse usando "não igual" ('==').
75
75
72
72
75
75
72
72
Não sou um programador, então peça desculpas antecipadamente, mas preciso de ajuda com isso. Geralmente, pode encontrar código e modificá-lo para funcionar, mas não neste caso. -p
Os arquivos Fastq possuem quatro linhas para cada leitura. Leia # 1 e, g, terá as seguintes 4 linhas:
@sample1
CGGCATCGTTTATGGTTGAGACTAGGACG
+
AAAAAEEEEEEEEEEEEEEEEEEEEEEEE
A primeira linha é o nome da amostra, a segunda linha é a sequência real, a terceira linha é um símbolo '+' e a quarta linha é um conjunto de "pontuações" ASCII para cada base na sequência. Cada base tem exatamente uma pontuação, portanto, o comprimento da linha 2 deve ser igual ao comprimento da linha quatro. Eu tinha analisado as linhas 2 e 4, procurando por pares de linha com comprimento desigual. Em vez disso, recebi o que parece que o par foi perdido.
Veja um exemplo de como um arquivo FASTQ pode parecer, com os pontos de interrogação representando as pontuações de qualidade perdidas ou não analisadas:
@sample1
CGGCATCGTTTATGGTTGAGACTAGGACG
+
AAAAAEEEEEEEEEEEEEEEEEEEEEEEE
@sample2
CCGGCTTCCGGTTCATCCCGCATCGCCAGTTC
+
@sample3
AAAA6E6/EEEEEEEE6/EE/EEAEEAA//E/
+
@sample4
ATTTCGGGGGGGGGGGGGG
+
??????????????????????????????????
@Sample5
GGTTAGCGCGCAGTTGGGCACCGTAACCCGGCTT
+
AAAAAEEEEEAEEEEEEEEEEEEEEEEEE//<EE
@sample6
CTAACCTGTCTCACGACGGTCTAAACCCAGCTCA
+
AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEE
Aqui está o que meus arquivos parsed (line2 + line4) pareciam:
CGGCATCGTTTATGGTTGAGACTAGGACG
AAAAAEEEEEEEEEEEEEEEEEEEEEEEE
CCGGCTTCCGGTTCATCCCGCATCGCCAGTTC
AAAA6E6/EEEEEEEE6/EE/EEAEEAA//E/
ATTTCGGGGGGGGGGGGGG
GGTTAGCGCGCAGTTGGGCACCGTAACCCGGCTT
AAAAAEEEEEAEEEEEEEEEEEEEEEEEE//<EE
CTAACCTGTCTCACGACGGTCTAAACCCAGCTCA
AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEE
Existem duas linhas de sequência consecutivas sem linha de pontuação de qualidade entre elas:
ATTTCGGGGGGGGGGGGGG
GGTTAGCGCGCAGTTGGGCACCGTAACCCGGCTT
Usando o código que você me deu:
awk 'NR%2==0 && length($0)!=last{print "Bad pair at lines",NR-1,"and",NR}{last=length($0)}' Fastq-seq-qual-parsed.txt
Bad pair at lines 5 and 6
OU: ./new-try.awk
Eu sugeriria
awk '
{ first = $0; getline; second = $0 }
length(first) != length(second) {
print "Error at line", NR-1
print first
print second
}
' file
Poderia usar o bash simples também, mas será muito mais lento:
nr=1
while IFS= read -r first; IFS= read -r second; do
if (( ${#first} != ${#second} )); then
printf "%s\n" "problem at line $nr" "$first" "$second"
fi
((nr+=2))
done < file
Tente:
awk 'NR%2==0 && length($0)!=last{print "Bad pair at lines",NR-1,"and",NR} {last=length($0)}' file
Vamos considerar isso como um arquivo de teste:
$ cat file
good123
good345
bad12
bad123
good_again
good_also1
Usando nosso comando, o par incomparável é identificado corretamente:
$ awk 'NR%2==0 && length($0)!=last{print "Bad pair at lines",NR-1,"and",NR} {last=length($0)}' file
Bad pair at lines 3 and 4
NR%2==0 && length($0)!=last{print "Bad pair at lines",NR-1,"and",NR}
Quando estamos em uma linha com numeração par, NR%2==0 , verificamos se o comprimento da linha é o mesmo da linha anterior. Se não for o mesmo, length($0)!=last , imprimimos uma mensagem.
last=length($0)
Isso salva o comprimento da linha atual na variável last .
Para aqueles que preferem o código espalhado por várias linhas:
awk '
NR%2==0 && length($0)!=last {
print "Bad pair at lines",NR-1,"and",NR
}
{
last=length($0)
}' file
Para imprimir, por exemplo, a linha 3 de um arquivo, podemos usar:
$ awk 'NR==3' file
bad12
Para imprimir um intervalo, digamos todas as linhas de 3 a 6, podemos usar:
$ awk 'NR>=3 && NR<=6' file
bad12
bad123
good_again
good_also1
Como alternativa, podemos obter resultados semelhantes do sed usando:
$ sed -n '3p' file
bad12
$ sed -n '3,6p' file
bad12
bad123
good_again
good_also1
Considere este arquivo de entrada:
$ cat File
@sample1
CGGCATCGTTTATGGTTGAGACTAGGACG
+
AAAAAEEEEEEEEEEEEEEEEEEEEEEEE
@sample2
CCGGCTTCCGGTTCATCCCGCATCGCCAGTTC
+
@sample3
AAAA6E6/EEEEEEEE6/EE/EEAEEAA//E/
+
@sample4
ATTTCGGGGGGGGGGGGGG
+
??????????????????????????????????
@Sample5
GGTTAGCGCGCAGTTGGGCACCGTAACCCGGCTT
+
AAAAAEEEEEAEEEEEEEEEEEEEEEEEE//<EE
@sample6
CTAACCTGTCTCACGACGGTCTAAACCCAGCTCA
+
AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEE
@sample7
CTAACCTGTCTCACGACGGTCTAAACCCAGCTCA
+
AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEE
Podemos detectar amostras inválidas, ou seja, amostras com comprimentos de linha desiguais ou com linhas secundárias que começam com ? , da seguinte forma:
$ awk '/^\+/{next} /^@/{s=$0;n=NR;next} prev{if(/^\?/ || length(prev)!=length($0)) printf "Sample %s (line %s) is bad:\n%s\n%s\n",s,n,prev,$0;prev="";next} {prev=$0}' File
Sample @sample4 (line 11) is bad:
ATTTCGGGGGGGGGGGGGG
??????????????????????????????????
Sample @sample7 (line 23) is bad:
CTAACCTGTCTCACGACGGTCTAAACCCAGCTCA
AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEE
Como alternativa, se quisermos ignorar amostras cuja segunda linha ('qualidade') começa com ? , então:
$ awk '/^\+/{next} /^@/{s=$0;n=NR;next} prev{if(!/^\?/ && length(prev)!=length($0)) printf "Sample %s (line %s) is bad:\n%s\n%s\n",s,n,prev,$0;prev="";next} {prev=$0}' File
Sample @sample7 (line 23) is bad:
CTAACCTGTCTCACGACGGTCTAAACCCAGCTCA
AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEE
Primeiro, faça um arquivo de teste em que as linhas 5 e 6 tenham comprimentos diferentes, então há algo para encontrar (o "
printf '%s\n' aaa aaa bbbb bbbb cccc ccc ddd ddd > foo
Resumo foo em dois arquivos virtuais usando bash substituição de processo e sed , onde cada caractere é substituído por . :
... então diff desses arquivos:
diff <(sed 's/././g' foo) <(sed -n '1~2{s/././g;p;p}' foo)
A saída mostra que a linha 6 não corresponde:
6c6
< ...
---
> ....
Se a saída acima for muito detalhada, os programas diff e kindred terão muitas opções ou poderão ser filtrados conforme necessário. Para mostrar apenas os números de linha:
diff <(sed 's/././g' foo) <(sed -n '1~2{s/././g;p;p}' foo) |
sed -n 's/c.*//p'
Saída:
6
Ou com mais detalhes, ou seja, linhas de arquivos originais não correspondentes numeradas:
f=foo
diff <(sed 's/././g' $f) <(sed -n '1~2{s/././g;p;p}' $f) |
sed -n 's/^\(.*\)c.*//p' | grep -B 1 -wf - <(cat -n $f)
Saída:
5 cccc
6 ccc