Hola,
una pregunta habitual cuando analizas un fichero de nucleótidos, por ejemplo un ensamblaje de un genoma, es qué porcentaje GC tiene. Asumiendo que el fichero está en formato FASTA, podemos obtener fácilmente ese valor con un mini-programa (one-liner) escrito en lenguaje perl. Por ejemplo, para el genoma comprimido de Brachypodium distachyon obtenido de Ensembl Plants, podríamos obtenerlo así:
zcat Brachypodium_distachyon.Brachypodium_distachyon_v3.0.dna.toplevel.fa.gz | \
perl -lne 'if(!/^>/){ $SQ=uc($_); while($SQ =~ /([ACTG])/g){ $stat{$1}++; $tot++ } }
END{ printf("%%GC=%1.1f\n",100*($stat{"G"}+$stat{"C"})/$tot);
foreach $nt (keys(%stat)){ print "$nt\t$stat{$nt}" } }'
%GC=46.4
A 72549289
T 72561114
C 62839311
G 62789747
Si quieres calcular el %GC solamente para ciertas regiones del genoma entonces puedes codificarlas en un fichero BED y usar bedtools nuc, como se explica en https://www.biostars.org/p/47047
Hasta pronto,
Bruno
El código es correcto, pero cuanto más código se pueda sacar de perl y dárselo a comandos del SO más rápido será. Aquí va una sugerencia quitando las cabeceras con grep, pasando todo a mayúsculas con tr, dividiendo las líneas en letras con fold, y eliminando las Ns con grep, dejando que perl sólo cuente las letras y rellene un hash para el resumen final:
ResponderEliminartime zgrep -v "^>" Brachypodium_distachyon.Brachypodium_distachyon_v3.0.dna.toplevel.fa.gz \
| tr '[:lower:]' '[:upper:]' | fold -w1 | grep -v N \
| perl -lne '$s{$_}++; $c++; END { print "@{[%s]}";
printf "%%GC = %.2f", 100*($s{G}+$s{C})/$c }'
Gracias
ResponderEliminar