☰

別のファイルの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からそのデータ、つまり「Start Gene」と「End Gene」の間にあるデータを抽出しようとします。また、希望の出力から「Coding Sequence」を除外したいと思います。

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

私はパールを使ってみました。

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;

おすすめ記事