複数の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;