☰

如何使用另一个文件中的 Id 从 gff3 文件中提取数据？

tag-icon • tag-icon perl
如何使用另一个文件中的 Id 从 gff3 文件中提取数据？

我有一个包含多个 ID 的文件

File 1:
g24007
g51692

和一个gff3文件如下

File2:
    # start gene g24007
    scaffold591 method  gene    3322458 3376057 0.41    -   .   ID=g24007
    scaffold591 method  transcript  3322458 3376057 0.41    -   .   ID=g24007.t1;Parent=g24007
    scaffold591 method  transcription_end_site  3322458 3322458 .   -   .   Parent=g24007.t1
    scaffold591 method  CDS 3323084 3323326 1   -   0   ID=g24007.t1.cds;Parent=g24007.t1
    # coding sequence = [atggaaaaagctaaagatggcgaagagagcccaagtgaggcatctcctccagcccaggtggggcttgaaaatatccctg
    # cgacggtgtctggggaggagggccagctgctgtatcacgaggagactatcgatcttggtggagacgagtttgggtctgaagagaatgaggaaccctca
    --
    # end gene g24007
    # start gene g20000
    scaffold591 method  gene    3322458 3376057 0.41    -   .   ID=g20000
    scaffold591 method  transcript  3322458 3376057 0.41    -   .   ID=g20000.t1;Parent=g20000
    ffold591    method  intron  3356166 3369049 1   -   .   Parent=g20000.t1
    scaffold591 method  CDS 3323084 3323326 1   -   0   ID=g20000.t1.cds;Parent=g20000.t1
    # coding sequence = [atggaaaaagctaaagatggcgaagagagcccaagtgaggcatctcctccagcccaggtggggcttgaaaatatccctg
    --
    # end gene g20000

在这里，我试图映射 file1 中的 Id 并从 file2 中提取相应的数据，即位于“起始基因”和“结束基因”之间的数据，同时我想从所需的输出中排除“编码序列”。

Expected output:
# start gene g24007
scaffold591 method  gene    3322458 3376057 0.41    -   .   ID=g24007
scaffold591 method  transcript  3322458 3376057 0.41    -   .   ID=g24007.t1;Parent=g24007
scaffold591 method  transcription_end_site  3322458 3322458 .   -   .   Parent=g24007.t1
scaffold591 method  CDS 3323084 3323326 1   -   0   ID=g24007.t1.cds;Parent=g24007.t1
# end gene g24007

我一直在尝试使用perl。

My code:
use strict;
use warnings;
use Data::Dumper;

my $file1 = 'IDs.txt';
open FILE1, "<", $file1 or die $!;
my $file2 = 'gff3.txt';

open FILE2, "<", $file2 or die $!;
my %id;
my @array;

while(<FILE1>)
{
    $id{$_} = 1;
}
#print Dumper \%id;

my $gene_id = 0;
while (<FILE2>)
{
    if($_ !~ /^#/)
    {
        @array = split(/\t/,$_);
        $array[8] =~ s/ID=//g;
        if($id{$_})
        {
            print $_, @array;
        }
    }
}
close FILE1;
close FILE2;

答案1

@Hari 没有查看您的预期输出，我尝试使用标准 gff3 文件。但是，我的脚本不打印“#startgene”和“#endgene”行。希望这对您有帮助

Code:
#!/usr/local/perl

use strict;
use warnings;

my $file1 = $ARGV[0];
my $file2 = $ARGV[1];
my $output_file = $ARGV[2];

my %id;
my $ctr = 0;
open(IN, $file1);
while(<IN>)
{

    $_ =~ s/\n|\r//g;
    $ctr++;
    $id{$_} = $ctr;
}
close IN;

open(IN, $file2);
open(OUT, ">".$output_file);
while(<IN>)
{
    $_ =~ s/\n|\r//g;

    if($_ !~ /^#/)
    {
        my @tmp = split(/\t/, $_);

        if($tmp[8] =~ /ID=g(\d+)/)
        {
            my $gene_id = "g".$1;
            if(exists $id{$gene_id})
            {
                print OUT $_."\n";
            }
        }
        elsif($tmp[8] =~ /Parent=g(\d+)\.t(\d+)/)
        {
            my $gene_id = "g".$1;
            if(exists $id{$gene_id})
            {
                print OUT $_."\n";
            }
        }
    }
}
close IN;
close OUT;

相关内容