#!/perl/bioinfo: GC

9 de marzo de 2021

contenido GC de un fichero FASTA

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