I have been doing a script that takes two files to extract a specific part of the data to make a new file. If you want to see the complete file, here's a GitHub link: enter link description here
File one (report file) is a type of file that reports me when a value is >=0.5 (column N°6 is the value that interests me). This file is something like this (this is only a part):
AGY29650_2_NA netOGlyc-4.0.0.13 CARBOHYD 2 2 0.0804934 . .
AGY29650_2_NA netOGlyc-4.0.0.13 CARBOHYD 4 4 0.0925522 . .
AGY29650_2_NA netOGlyc-4.0.0.13 CARBOHYD 13 13 0.0250116 . .
AGY29650_2_NA netOGlyc-4.0.0.13 CARBOHYD 23 23 0.565981 . .
...
File two (fasta file) is a type of file used in bioinformatics, here is how looks it (this is only a part):
>AGY29650.2|NA spike protein
MTYSVFPLMCLLTFIGANAKIVTLPGNDA...EEYDLEPHKIHVH*
The purpose of my script is to take the column 1 and 4 when a value is >=0.5 in the column N°6, for example, the row N°4 is a #POSITIVE value, therefore the script takes the value in the column N°1 (AGY29650_2_NA, this is an ID) and the value in column N°4, 23 (position). Then the script searches to match the ID (AGY29650_2_NA) in file two (fasta file) with the ID in this file AGY29650.2 and then looking for position 23 in the data, for example, the letter T in position 23:
MTYSVFPLMCLLTFIGANAKIV T LP
Then, the script print position 23 with 2 letters to the left and 2 letters to the right, the output:
IVTLP
The script is incomplete, however, here is the first problem that I couldn’t solve yet. The IDs between files have a few differences, for example:
AGY29650_2_NA (file one) and AGY29650.2 (file two)
To resolve this, a colleague recommends that I use regular expressions to select the ID in each file, for example:
s/^\s*([^_]+)_([0-9]+)_([a-zA-Z0-9]+)/$1.$2|$3/
My second problem is that I couldn't solve how to incorporate this regex into the script, I´m thinking in a foreach loop, maybe. My third problem is a certificate if the script is really searching for the position (column N°4) and take the residues adjacent (two-letter to left and two to the right) for the final output. Here is the incomplete script:
use strict;
use warnings;
use Bio::SeqIO;
my $file = $ARGV[0];
my $in = $ARGV[1];
my %fastadata = ();
my @array_residues = ();
my $seqio_obj = Bio::SeqIO->new(-file => $in,
-format => "fasta" );
while (my $seq_obj = $seqio_obj->next_seq ) {
my $dd = $seq_obj->id;
my $ss = $seq_obj->seq;
###my $ee = $seq_obj->desc;
$fastadata{$dd} = "$ss";
}
my $thres = 0.5; ### Selection of values in column N°5 with the following condition: >=0.5
# Open file
open (F, $file) or die; ### open the file or end the analyze
while(my $one = <F>) {### readline => F
$one =~ s/\n//g;
$one =~ s/\r//g;
my @cols = split(/\s+/, $one); ### split columns
next unless (scalar (@cols) == 7); ### the line must have 7 columns to add to the array
my $val = $cols[5];
if ($val >= 0.5) {
my $position = $cols[3];
my $id_list = $cols[0];
if (exists($fastadata{$id_list})) {
my $new_seq = $fastadata{$id_list};
my $subresidues = substr($new_seq, $position -3, 6);
}
}
}
close F;
I´m looking for help to incorporate the regex in the script and then to print the output that I´m looking for.
Any idea or comment is welcome.