-1

I have a file that contains several IDs

File 1:
g24007
g51692

And a gff3 file as follows

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

Here I am trying to map the Ids from file1 and extract the corresponding data from file2 i.e lying between "start gene" and "end gene" and at the same time I want to exclude the "coding sequence" from my desired output.

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

I have been trying with 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;
  • can you accept other solutions except perl? – RomanPerekhrest Nov 5 '17 at 11:12
  • @RomanPerekhrest yes. awk or Perl commands will be really helpful. – Hari Nov 5 '17 at 11:15
2

@Hari without looking at your expected output I tried with a standard gff3 file. However, my script doesn't print the "# start gene" and "# end gene" lines. Hope this will be helpful for you

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;
  • I wanted to use it on standard gff3 file. It's working. – Hari Nov 5 '17 at 18:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.