El problema más habitual y sencillo es la redundancia de secuencias, cuando en un mismo fichero aparecen secuencias idénticas con identificadores diferentes. En otras ocasiones el problema es el contrario, identificadores idénticos apuntan a secuencias parciales de una misma secuencia (especialmente en datos de proteómica o transcriptómica).
El siguiente script filtra un fichero de sequencias en formato FASTA, dos tipos de filtro son posibles: 1) filtrar las secuencias redundantes (repetidas), 2) filtrar las secuencias de transcritos parciales (recuperando únicamente los transcritos de mayor secuencia en ficheros de proteínas)
#!/usr/bin/perl
my ($infile) = @ARGV;
my $outfile = filter_fasta_file($infile,['longtrans','nonredundant']);
# Create a FASTA file with custom filtered sequences
sub filter_fasta_file {
my ($infile, $filter, $outfile) = @_;
if (!defined($outfile)){$outfile=$infile.'.filtered'}
open(INFILE,"$infile");
open(OUTFILE,">$outfile");
my ($sequences,$names,$headers) = read_FASTA_file($infile);
my (@previous_sequences,@previous_names,@previous_headers);
# my (@filtered_sequences,@filtered_names,@fitered_headers);
for (my $i=0; $i<=$#{$headers}; $i++){
my $write=1;
if (in_array($filter,'nonredundant')){
my ($found,$posic)=in_array($sequences, $sequences->[$i], 1);
foreach my $pos (@{$posic}){
if ($i != $pos && $i > $pos){
$write = 0; last;
}
}
}
if (in_array($filter,'longtrans')){
$names->[$i] =~ /([^\.]+)\.?(\d?)/;
my $name = $1;
my $variant = $2;
my ($found,$posic)=in_array($names, "/$1/", 1);
foreach my $pos (@{$posic}){
if ($i != $pos){
if (length($sequences->[$i]) < length($sequences->[$pos])) {
$write = 0; last;
}
if (length($sequences->[$i]) == length($sequences->[$pos]) && $sequences->[$i] eq $sequences->[$pos]){
$names->[$pos] =~ /([^\.]+)\.?(\d?)/;
my $name2 = $1;
my $variant2 = $2;
if ($variant2 < $variant){
$write = 0; last;
}
}
}
}
}
if ($write) {
print OUTFILE $headers->[$i]."\n".$sequences->[$i]."\n";
}
}
close INFILE;
close OUTFILE;
return $outfile;
}
#################################################################################
# Create 3 array references with sequences, names and headers from a FASTA file
sub read_FASTA_file
{
my($FASTAfile,$full_header,$allow_empty_sequences,$allow_multi_rows) = @_;
if(!$full_header){ $full_header = $allow_empty_sequences = 0 }
if(!$allow_empty_sequences){ $allow_empty_sequences = 0 }
if(!$allow_multi_rows){ $allow_multi_rows = 0 }
my (@sequences,@names,@headers,$seq,$seqname,$header);
my $headerOK = 0;
open(FASTA,$FASTAfile) || die "# read_FASTA_file : cannot find $FASTAfile\n";
$seq = '';
while(<FASTA>){
next if(/^#/); # allow comments
s/\012\015?|\015\012?|^\s+|\s+$//;
if(/>/){
$headerOK++;
if($seq || ($allow_empty_sequences && $seqname)){
if(!$allow_multi_rows){
$seq =~ s/\d|\s//g;
}
push(@sequences,$seq); # store this FASTA sequence
push(@names,$seqname);
push(@headers,$header);
$seq='';
}
if($full_header){
$seqname = substr($_,1);
}else { # get next sequence name
if (/^>\s*([^\t|^\||^;]+)/){
$header = $_;
$seqname = $1;
if (length($seqname)>32 && $seqname =~ /^([^\s]+)+\s/){
$seqname = $1;
}
$seqname =~ s/^\s+|\s+$// ;
}
}
}else{
if(!$allow_multi_rows){
s/\s+//g;
}
$seq .= $_;
}
}
close(FASTA);
# take care of last FASTA sequence
push(@sequences,$seq);
push(@names,$seqname);
push(@headers,$header);
if(!$headerOK){ @sequences = @names = @headers = (); } #Aug2007
return (\@sequences,\@names,\@headers);
}
#################################################################################
# Searchs a string inside an array
sub in_array {
my ($array,$search_for,$search_posic) = @_;
# my %items = map {$_ => 1} @$array; # create a hash out of the array values
# my $exist = (exists($items{$search_for}))?1:0;
# return (exists($items{$search_for}))?1:0;
my @posic;
# If search a pattern
if ($search_for =~ /^\/.+\/$/){
$search_for =~ s/^\/|\/$//g;
@posic = grep { $array->[$_] =~ /$search_for/ } 0 .. $#{$array};
} else {
@posic = grep { $array->[$_] eq $search_for } 0 .. $#{$array};
}
if (defined($search_posic) && $search_posic){
(scalar @posic)?return (1,\@posic):return (0);
} else {
(scalar @posic)?return 1:return 0;
}
}
An example FASTA file to test the script:
>AT1G51370.2 | Symbols: | F-box family protein | chr1:19045615-19046748 FORWARD
MVGGKKKTKICDKVSHEEDRISQLPEPLISEILFHLSTKDSVRTSALSTKWRYLWQSVPG
LDLDPYASSNTNTIVSFVESFFDSHRDSWIRKLRLDLGYHHDKYDLMSWIDAATTRRIQH
LDVHCFHDNKIPLSIYTCTTLVHLRLRWAVLTNPEFVSLPCLKIMHFENVSYPNETTLQK
LISGSPVLEELILFSTMYPKGNVLQLRSDTLKRLDINEFIDVVIYAPLLQCLRAKMYSTK
NFQIISSGFPAKLDIDFVNTGGRYQKKKVIEDILIDISRVRDLVISSNTWKEFFLYSKSR
PLLQFRYISHLNARFYISDLEMLPTLLESCPKLESLILVMSSFNPS*
>AT1G70790.1 | Symbols: | C2 domain-containing protein | chr1:26700724-26702127 FORWARD
MEDKPLGILRVHVKRGINLAIRDATTSDPYVVITLANQKLKTRVINNNCNPVWNEQLTLS
IKDVNDPIRLTVFDKDRFSGDDKMGDAEIDFRPFLEAHQMELDFQKLPNGCAIKRIRPGR
TNCLAEESSITWSNGKIMQEMILRLKNVECGEVELMLEWTDGPGCKGLGREGSKKTPWMP
TKRLD*
>AT1G50920.1 | Symbols: | GTP-binding protein-related | chr1:18870555-18872570 FORWARD
MVQYNFKRITVVPNGKEFVDIILSRTQRQTPTVVHKGYKINRLRQFYMRKVKYTQTNFHA
KLSAIIDEFPRLEQIHPFYGDLLHVLYNKDHYKLALGQVNTARNLISKISKDYVKLLKYG
DSLYRCKCLKVAALGRMCTVLKRITPSLAYLEQIRQHMARLPSIDPNTRTVLICGYPNVG
KSSFMNKVTRADVDVQPYAFTTKSLFVGHTDYKYLRYQVIDTPGILDRPFEDRNIIEMCS
ITALAHLRAAVLFFLDISGSCGYTIAQQAALFHSIKSLFMNKPLVIVCNKTDLMPMENIS
EEDRKLIEEMKSEAMKTAMGASEEQVLLKMSTLTDEGVMSVKNAACERLLDQRVEAKMKS
KKINDHLNRFHVAIPKPRDSIERLPCIPQVVLEAKAKEAAAMEKRKTEKDLEEENGGAGV
YSASLKKNYILQHDEWKEDIMPEILDGHNVADFIDPDILQRLAELEREEGIREAGVEEAD
MEMDIEKLSDEQLKQLSEIRKKKAILIKNHRLKKTVAQNRSTVPRKFDKDKKYTTKRMGR
ELSAMGLDPSSAMDRARSKSRGRKRDRSEDAGNDAMDVDDEQQSNKKQRVRSKSRAMSIS
RSQSRPPAHEVVPGEGFKDSTQKLSAIKISNKSHKKRDKNARRGEADRVIPTLRPKHLFS
GKRGKGKTDRR*
>AT1G70795.1 | Symbols: | C2 domain-containing protein | chr1:26700724-26702127 FORWARD
MEDKPLGILRVHVKRGINLAIRDATTSDPYVVITLANQKLKTRVINNNCNPVWNEQLTLS
IKDVNDPIRLTVFDKDRFSGDDKMGDAEIDFRPFLEAHQMELDFQKLPNGCAIKRIRPGR
TNCLAEESSITWSNGKIMQEMILRLKNVECGEVELMLEWTDGPGCKGLGREGSKKTPWMP
TKRLD*
>AT1G70790.2 | Symbols: | C2 domain-containing protein | chr1:26700724-26702127 FORWARD
MEDKPLGILRVHVKRGINLAIRDATTSDPYVVITLANQKLKTRVINNNCNPVWNEQLTLS
IKDVNDPIRLTVFDKDRFSGDDKMGDAEIDFRPFLEAHQMELDFQKLPNGCAIKRIRPGR
TNCLAEESSITWSNGKIMQEMILRLKNVECGEVELMLEWTDGPGCKGLGREGSKKTPWMP
TKRLD*
>AT1G36960.1 | Symbols: | unknown protein | chr1:14014796-14015508 FORWARD
MTRLLPYKGGDFLGPDFLTFIDLCVQVRGIPLPYLSELTVSFIAGTLGPILEMEFNQDTS
TYVAFIRVKIRLVFIDRLRFFRREEAAASNTITDQTHMTSSNSSDISPASPISQPPLPAS
LPSHDSYFDAGIQASRLVNPRAISQHHFSSSYSDFKGKEKAKIKIGECSKRKKDKQVDSG
T*
>AT1G51370.1 | Symbols: | F-box family protein | chr1:19045615-19047141 FORWARD
MVGGKKKTKICDKVSHEEDRISQLPEPLISEILFHLSTKDSVRTSALSTKWRYLWQSVPG
LDLDPYASSNTNTIVSFVESFFDSHRDSWIRKLRLDLGYHHDKYDLMSWIDAATTRRIQH
LDVHCFHDNKIPLSIYTCTTLVHLRLRWAVLTNPEFVSLPCLKIMHFENVSYPNETTLQK
LISGSPVLEELILFSTMYPKGNVLQLRSDTLKRLDINEFIDVVIYAPLLQCLRAKMYSTK
NFQIISSGFPAKLDIDFVNTGGRYQKKKVIEDILIDISRVRDLVISSNTWKEFFLYSKSR
PLLQFRYISHLNARFYISDLEMLPTLLESCPKLESLILEMVKNQSTRRHGEKREPNVMVS
TVPWCLVSSLKFVELKRSIPRYEGEMELVRYVLTNSTVLKKLRLNVYYTKKAKCAFLTEL
VAIPRCSSTCVVLVL*
Hay un error en el código, creo.
ResponderEliminarLa expresión regular s/\012\015?|\015\012?//; se aplica a la variable $_, pero luego, ese nuevo valor transformado, no es utilizado.
Es posible que esa línea deba aplicarse antes de $seq .= $_;
Si antes no se ha notado es porque, realmente, esa línea no tiene efecto alguno en el programa: la línea que filtra los caracteres en blanco, incluidos los caracteres \015 y \012, es la de $seq =~ s/\s+//g;. Incluso sobra la de chomp $seq;.
Lo he corregido, la expresión 's/\012\015?|\015\012?//' la suelo usar en lugar de 'chomp' para eliminar saltos de línea (así evito problemas con archivos creados en algunos editores de windows), pero como bien dices había algo de código redundante.
ResponderEliminarEs mejor escribirlo como
ResponderEliminar$seq =~ s/\s+$//;
según se indica en el perlfaq4, en la cuestión «How do I strip blank space from the beginning/end of a string?».
O hacerlo solo una vez. Es decir, vas acumulando en $seq .= $_; y luego, cuando termine el bucle, haces un
$seq =~ s/\s+$//gm;
y eliminará todos los finales de línea de toda la secuencia.
Caso extremo: si sabemos que la $seq no ha de tener ningún espacio, pues ya lo sabemos hacer:
$seq =~ s/\s+//g;
y solo hay que hacerlo una vez, al final.