Beefy Boxes and Bandwidth Generously Provided by pair Networks
Your skill will accomplish
what the force of many cannot
 
PerlMonks  

bioperl module to extract specific nucleotides given chromosome and exact location of that nucleotide

by xxArwaxx (Novice)
on Mar 25, 2014 at 22:47 UTC ( #1079751=perlquestion: print w/replies, xml ) Need Help??

xxArwaxx has asked for the wisdom of the Perl Monks concerning the following question:

Hello, folks;

so, I'm working on a project that I'm trying to understand more about the genotyping data that was given to us(.gen files), I wrote a perl script to extract the lines I'm interested in from the given genotyping data and saved them in lines in a comma-delimited txt file(CVS format), the lines in .txt file are in the format of:

physical_location, Allele 1, Allele2

34787638,A,C

34787800,A,G

the big question is:

are the given locations in .gen files do correspond to the same alleles(1,2) in my txt file or not? to do that, I need to search for that exact physical_location (in each line) against any genome database(NCBI, Ensembl..etc) to retrieve the alleles(1,2) that fall in that exact location (for specific chromosome).

Using genome browser to do that manually is time consuming, and I believe it is a common task so there should be a BioPerl module to retrieve the alleles in a specific location for a given chromosome.

any ideas if there is a BioPerl module that can do that, or how to approach this problem ?

EDIT(more info):

the genotyping data that was given to us was in .gen files in which each file represents the genotype info for ONE chromosome, so I have a total of 22 .gen files (sizes between 700 MB up to 4 GB). Each line in a .gen file is in the format:

(Chromsome, MarkerID, Ph_location, Allele1, Allele2, .......some other irrelevant info)

Then there is another small .txt file that have the 'list of genes' that I'm interested in(50 genes) where each line is a gene, in the following format (GeneName, Chromosome, StartPosition, EndPosition).

The Perl script I wrote was to extract the lines I'm interested in from the .gen file (which are the lines that have a Ph_location (now experimenting on the .gen for chr1) that falls in the range from StartPosition to EndPosition from 'list of gene' .txt file (looking to those lines for chr1, it happens to be 10 genes) and saved them as lines in a comma-delimited txt file(CVS format) that have the format I mentioned in the very beginning of this post before the editing part(with location, allele1, allele2 format).

Now, I have a .txt file with the lines I'm interested in, again with location, allele1, allele2 format. NOW, I'm searching for a BioPerl module to retrieve (from any genome database(NCBI, Ensembl..etc)) to retrieve the alleles(1,2) that correspond to that location (from the .txt file)and report these corresponding alleles in another text file along with their location, so I finally end up with two .txt files, one (from genotype extracted lines), that I'll pull the locations from to use it to search for alleles1,2, and another .txt file that will have the location and the corresponding alleles after the search. Our goal is to validate the alleles given to us by retrieving those from ncbi or wherever using their locations, to eventually try to draw a different approach to classify the data moving forward.

I apologize for the lengthy post, I'm trying to make it more clear to get the best out of this discussion. I hope it is more clear now. Thanks in advance! :)

FINAL OUTPUT

two .txt files, one with the lines I'm interested in in the format (physical_location, Allele 1, Allele2), file 1 looks like (this is in chromosome 1, in case anyone wanted to run this against ncbi or any genome database):

34787638,A,C

34788686,A,G

34789549,C,T

34789695,C,G

34789808,C,T

347890859,C,G

then another .txt file, file2.txt that has the same first column(the locations), I need the other two colomns(Allele1,Allele2) for all the locations in this file these need to be retrieved from any genome database(NCBI for example) by a bioperl module that goes through each line in file1 extract the first column, go to ncbi, fetch the corresponding alleles. this module should accept two arguments; the location, and (chromosome number); in which chromosome it should fetch the corresponding alleles from. That's our goal, so we can doublecheck if the given alleles in file1 actually falls in the given locations.

  • Comment on bioperl module to extract specific nucleotides given chromosome and exact location of that nucleotide

Replies are listed 'Best First'.
Re: bioperl module to extract specific nucleotides given chromosome and exact location of that nucleotide
by biohisham (Priest) on Mar 26, 2014 at 02:06 UTC
    So you basically have some coordinates and want to check if anyone of the alleles matches the location given by the coordinates right ? But your coordinates don't tell about which chromosome these alleles are expected but what you can probably do is download the chromosomes in say FastA format and then use Bio::SeqIO using the subseq function.

    here is the documentation to the module Bio::SeqIO

    The subseq function requires a beginning and an end to what it is going to extract so probably you want to reformat your coordinates as follows

    $seq->subseq(coordinate-1, coordinate+1)

    to extract three bases with the middle one being the actual nucleotide you want to compare your alleles against.

    Note that BioPerl can be a bit slow and therefore consider using something like sfetch or esl-sfetch which will again require providing beginning and end positions to what it is to extract


    David R. Gergen said "We know that second terms have historically been marred by hubris and by scandal." and I am a four y.o. monk...
      I'll try this out and see! Still, more suggestions are welcomed! I was searching earlier and was able to find something that has to do with exact location being a method for a 'SeqFeature' object, and didn't have anything to do with chromosome number. I'm interested in finding allele1 and allele2 for that position.
Re: bioperl module to extract specific nucleotides given chromosome and exact location of that nucleotide
by kcott (Bishop) on Mar 26, 2014 at 07:52 UTC

    G'day xxArwaxx,

    "any ideas if there is a BioPerl module that can do that, or how to approach this problem ?"

    I'm not familiar with the BioPerl modules; however, based on your problem description, a script to achieve this is very straightforward:

    #!/usr/bin/env perl -l use strict; use warnings; use autodie; use Text::CSV; my $ref_file = 'pm_1079751_reference.csv'; my @gen_files = qw{pm_1079751_source1.gen pm_1079751_source2.gen}; my $allele_sep = '-'; my %ref_data; my $ref_csv = Text::CSV::->new; open my $ref_fh, '<', $ref_file; while (my $row = $ref_csv->getline($ref_fh)) { $ref_data{$row->[0]} = join $allele_sep, @$row[1,2]; } close $ref_fh; for my $gen_file (@gen_files) { my $gen_csv = Text::CSV::->new({allow_whitespace => 1}); open my $gen_fh, '<', $gen_file; while (my $row = $gen_csv->getline($gen_fh)) { next unless exists $ref_data{$row->[2]}; if ($ref_data{$row->[2]} eq join $allele_sep, @$row[3,4]) { print join(',' => @$row[2,3,4]), "\n\tfound in '$gen_file':\n\t\t", join ', ' => @$ +row; } } close $gen_fh; }

    You stated the file you've already created is in CSV format, so I've used the Text::CSV module to parse it.

    [Side note: CSV files have, by convention, a .csv extension, not a .txt extension (as you indicated you used).]

    Your description of the .gen files seems to indicate a CSV format so I've used Text::CSV for these as well. If that's a false assumption on my part, you'll need to make some code changes; however, the logic I've used shouldn't need reworking.

    Using these test files:

    $ cat pm_1079751_reference.csv 34787638,A,C 34787800,A,G
    $ cat pm_1079751_source1.gen C1, M1, 99999999, A, C, ... source1 info ... C1, M2, 34787638, A, G, ... source1 info ... C1, M3, 34787800, A, C, ... source1 info ... C1, M4, 34787800, A, G, ... source1 info ...
    $ cat pm_1079751_source2.gen C1, M1, 99999999, A, C, ... source2 info ... C1, M2, 34787638, A, G, ... source2 info ... C1, M4, 34787800, A, G, ... source2 info ... C1, M3, 34787638, A, C, ... source2 info ...

    I get this output:

    34787800,A,G found in 'pm_1079751_source1.gen': C1, M4, 34787800, A, G, ... source1 info ... 34787800,A,G found in 'pm_1079751_source2.gen': C1, M4, 34787800, A, G, ... source2 info ... 34787638,A,C found in 'pm_1079751_source2.gen': C1, M3, 34787638, A, C, ... source2 info ...

    -- Ken

      I appreciate your response, but I believe that it should be much more simpler than this, check the FINAL OUT portion in the post that I've just added, it explains it much more. Thanks!
        "I appreciate your response, but I believe that it should be much more simpler than this, check the FINAL OUT portion in the post that I've just added, it explains it much more. Thanks!"

        So, two days after you posted your question you decide to tell us the exact output format you want. You then have a problem with my sample output which, in the absence of any other information, I had to guess at.

        If you're incapable of modifying

        print join(',' => @$row[2,3,4]), "\n\tfound in '$gen_file':\n\t\t", join ', ' => @$row;

        which outputs

        34787800,A,G found in 'pm_1079751_source1.gen': C1, M4, 34787800, A, G, ... source1 info ...

        to code which just outputs the first line of that, i.e.

        34787800,A,G

        I suggest you read "Perl introduction for beginners".

        -- Ken

Re: bioperl module to extract specific nucleotides given chromosome and exact location of that nucleotide
by AnomalousMonk (Bishop) on Mar 26, 2014 at 04:48 UTC

    You bio-guys sure love underlining. You've managed to underline almost the entire page, something I'm a bit surprised is possible. In any event, I suspect the perpetrator is the non-closing underline tag at the end of  <u><b>EDIT(more info):</b><u> which you might consider fixing.

Re: bioperl module to extract specific nucleotides given chromosome and exact location of that nucleotide
by frozenwithjoy (Priest) on Mar 25, 2014 at 22:53 UTC
    Can you provide more information? Are you looking up the genotype data in a bam file, fasta file, or?
      Done! Check more info section in the post.

Log In?
Username:
Password:

What's my password?
Create A New User
Domain Nodelet?
Node Status?
node history
Node Type: perlquestion [id://1079751]
Approved by Old_Gray_Bear
Front-paged by biohisham
help
Chatterbox?
and the web crawler heard nothing...

How do I use this? | Other CB clients
Other Users?
Others wandering the Monastery: (2)
As of 2021-12-01 08:39 GMT
Sections?
Information?
Find Nodes?
Leftovers?
    Voting Booth?

    No recent polls found

    Notices?