完成脚本在两个文件中搜索并提取一段数据的想法

Ideas to complete a script to search in two files and extract a section of data

我一直在做一个脚本,它需要两个文件来提取数据的特定部分以创建一个新文件。 如果你想看完整的文件,这里有一个 GitHub link: enter link description here

文件一(报告文件)是一种文件,当一个值>=0.5时报告我(第6列是我感兴趣的值)。 这个文件是这样的(这只是一部分):

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        .       .      
...

文件二(fasta文件)是生物信息学中使用的一种文件,下面是它的样子(这只是一部分):

>AGY29650.2|NA spike protein
MTYSVFPLMCLLTFIGANAKIVTLPGNDA...EEYDLEPHKIHVH*

我的脚本的目的是当第 N°6 列中的值 >=0.5 时取第 1 列和第 4 列,例如,第 N°4 行是#POSITIVE 值,因此脚本采用N°1 列中的值(AGY29650_2_NA,这是一个 ID)和 N°4, 23 列中的值(位置)。 然后脚本搜索将文件二(fasta文件)中的ID(AGY29650_2_NA)与本文件AGY29650.2中的ID匹配,然后在数据中查找位置23,例如字母T在位置23 :

MTYSVFPLMCLLTFIGANAKIV T LP

然后,脚本打印位置23,左2个字母,右2个字母,输出:

IVTLP

脚本不完整,但是,这是我还没有解决的第一个问题。文件之间的 ID 有一些不同,例如:

AGY29650_2_NA (file one) and AGY29650.2 (file two)

为了解决这个问题,同事建议我使用正则表达式来select每个文件中的ID,例如:

s/^\s*([^_]+)_([0-9]+)_([a-zA-Z0-9]+)/.|/

我的第二个问题是我无法解决如何将此正则表达式合并到脚本中的问题,我可能在考虑 foreach 循环。 我的第三个问题是证书,如果脚本真的在搜索位置(第 4 列)并取相邻的残基(左边两个字母和右边两个字母)作为最终输出。 这是不完整的脚本:

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;
​

我正在寻求帮助以将正则表达式合并到脚本中,然后打印我正在寻找的输出。

欢迎任何想法或评论。

未经测试(因为您没有 post MRE),但这应该有效:

    my $position = $cols[3];
    my $id_list = $cols[0];
    $id_list =~ s/^\s*([^_]+)_([0-9]+)_([a-zA-Z0-9]+)/.|/;   # Add this line
    if (exists($fastadata{$id_list})) {

这会修改 $id_list 变量,使其与您的哈希键兼容。

perl -lane 'if (!$#ARGV) { $h{$F[0]}=$F[3]; next }
            $F[0]=~s/^(\w+)_(\d+)_(\w+)/.|/;
            if ( $h{$F[0]} && $F[5] > 0.5 ) 
            { print substr $h{$F[0]}, $F[4]-3, 5 }' fasta.txt report.txt

-a 将记录的字段拆分为空格并填充特殊数组 @F$F[0] 是第一个字段,$F[1] 是第二个字段,依此类推.

如果这是第一个文件,fasta.txtif (!$#ARGV)

创建散列 %hfasta.txtAGY29650.2|NA 的第一个字段作为键:$h{$F[0]}

给它赋值fasta.txt的第4个字段,MTYSVFPLMCLLT...作为值:$h{$F[0]}=$F[3]

***这里结束fasta.txt的处理,开始下一个文件的处理,report.txt***

使 report.txtAGY29650_2_NA 的第一个字段看起来与上一个文件 AGY29650.2|NA 的相同:$F[0]=~s/^(\w+)_(\d+)_(\w+)/.|/

如果report.txt的第一个字段$F[0]作为键存在于%h中并且report.txt的第6个字段大于0.5:if ( $h{$F[0]} && $F[5] > 0.5 )

取子串MTYSVFPLMCLLT...substr $h{$F[0]}report.txt 的第 5 个字段指定的偏移量减去 3:$F[4]-3,对于接下来的 5 个字符并打印它们:IVTLP.