SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
SAM format MD tag with gaps in reference? nZyMe Bioinformatics 4 10-22-2014 12:16 AM
What tools can convert sequence file from tabular format to fasta format? yangjianhunt Bioinformatics 5 03-26-2014 01:48 PM
Appropriate Reference Format HSG_student1 Genomic Resequencing 1 04-25-2011 11:23 PM
Bowtie MD tag format reference mendl7 Bioinformatics 2 02-23-2011 01:47 PM
SAMtool's pileup format - reference base question mfischer Bioinformatics 11 10-19-2010 04:12 PM

Reply
 
Thread Tools
Old 09-22-2015, 07:17 AM   #1
mcrawford
Member
 
Location: Brighton

Join Date: Jun 2015
Posts: 10
Default Need my reference sequence to be in a specific format

I need my reference genome (S. cerevisiae strain S288c) to be in a specific format for some analysis but can't find anything matching what I need online.

i.e. the data I have is in this format:
>I [org=Saccharomyces cerevisiae] [strain=S288C] [chromosome=I]
CCACACCACACCCACACACCCAC...etc

But I need it to be in a table like:

chr pos seq
1 1 C
1 2 C
1 3 A
1 4 C
etc.

I'm hoping that one of the following will be true, please let me know if you can help... Thanks!
1. The data already exists somewhere online
2. There's a program that can make this file for me out of my reference sequence (possibly PySAMstats?)
3. Some hints to write my own R script to make the file?

Last edited by mcrawford; 09-22-2015 at 07:53 AM.
mcrawford is offline   Reply With Quote
Old 09-22-2015, 08:19 AM   #2
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,982
Default

There will always be cooler ways of doing this. What follows happens to be one way. I am assuming chromosome 1 is followed by 2 etc in the genome file.
Use at your own risk

Save the code below in a file (e.g. table.pl).

Code:
use strict;
use warnings;

my $datafile = $ARGV[0];
my $outfile = $ARGV[1];

open (IN, $datafile) or die "can't open the datafile: $datafile\n";
open (OUT, ">$outfile") or die "can't open the outputfile: $outfile\n";

my $count = 1;
my $chrom = 0;

while(my $line = <IN>){
        if ($line =~ /^>/){
                $chrom++;
                $count = 1;
                }
        else {
                chomp $line;
                my @nuc = split(//, $line);
                while (@nuc){
                        my $a = pop(@nuc);
                        print OUT $chrom." ".$count." ".$a."\n";
                        $count++;
                        }
                }

}
close IN;
close OUT;
Run like this

Code:
$ perl table.pl genome.fa table.txt
GenoMax is offline   Reply With Quote
Old 09-23-2015, 01:06 AM   #3
mcrawford
Member
 
Location: Brighton

Join Date: Jun 2015
Posts: 10
Default

Quote:
Originally Posted by GenoMax View Post
There will always be cooler ways of doing this. What follows happens to be one way. I am assuming chromosome 1 is followed by 2 etc in the genome file.
Use at your own risk
Thanks so much GenoMax, this worked perfectly and produced exactly what I wanted! You're the best
mcrawford is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 08:24 PM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2019, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO