SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Help needed to generate Consensus from alignment sarbashis Illumina/Solexa 0 08-10-2011 09:58 PM
how to generate ini file cjose Bioinformatics 4 07-22-2011 10:18 AM
Looking for tool to generate GFF splice graph files from RNA seq data desperado Bioinformatics 0 07-07-2011 03:17 AM
Generate subset from BAM format ForeignMan Genomic Resequencing 4 02-28-2011 06:20 AM
generate CIGAR string from 2 sequences? bbimber Bioinformatics 0 03-20-2010 09:44 AM

Reply
 
Thread Tools
Old 10-26-2011, 02:22 PM   #1
ykdang
Junior Member
 
Location: Dallas

Join Date: Oct 2011
Posts: 7
Default How to generate GFF?

I am working on a fungus call Neurosprora. Previous lab member built a gbrowser (1.70) to visualize the deep-seq data with GFF2 file format. Now I obtain some new deep-seq data and want to load it onto that browser. I heard that he used bowtie to do alignment and convert into gff2, but I cannot find the script to convert. Does anyone have similar script which could share? Or just indicate any tools can do this trick? I am not a perl person and don't think I could learn it in a short while.
ykdang is offline   Reply With Quote
Old 10-27-2011, 07:39 AM   #2
gringer
David Eccles (gringer)
 
Location: Wellington, New Zealand

Join Date: May 2011
Posts: 836
Default

If you have a BAM file (which bowtie can produce), this tutorial suggests that GBrowse can handle BAM directly:

http://gmod.org/wiki/GBrowse_NGS_Tut..._the_SAM_Files

You could also have a look at the admin tutorial:

http://gmod.org/gbrowse2/tutorial/tutorial.html#graph
gringer is offline   Reply With Quote
Old 10-27-2011, 09:33 AM   #3
gringer
David Eccles (gringer)
 
Location: Wellington, New Zealand

Join Date: May 2011
Posts: 836
Default

FWIW, I've just today ended up writing a python script to convert from SAM/BAM files to GFF3 files using pysam. The code may be useful for you if you can't find anything else suitable for your conversion.

For a bit of context, I broke up genomic contigs into 100bp fragments, naming the sequences <contig>#<start>-<end>, then used bowtie2 to map them to the genome -- that's why I've got the 'read.qname.find' bits in the code. My contigs started with 'v', which bowtie replaced with 'N' for some odd reason, so I had to do a bit of extra fiddling to add the 'v' back in.

Here's the relevant part of my code which does the SAM->GFF conversion:
Code:
samFile = pysam.Samfile(samFileName, "r")
totalCount = 0
sys.stderr.write("Getting reads from SAM file...")
sys.stdout.write("##gff-version 3\n")
gffWriter = csv.writer(sys.stdout, delimiter = '\t')
for read in samFile:
    totalCount += 1
    if(read.tid > 0):
        qContig = 'v' + read.qname[1:read.qname.find("#")]
        qContigStart = read.qname[read.qname.find("#")+1:read.qname.find("-")]
        qContigEnd = read.qname[read.qname.find("-")+1:]
        tContig = samFile.getrname(read.tid)
        strand = "-" if read.is_reverse else "+"
        score = "."
        for tag in read.tags:
            if(tag[0] == "XS"):
                score = str(tag[1]+1000)
        gffWriter.writerow((qContig,
                            "sam2gff3-"+os.path.basename(samFileName),
                            "nucleotide_match", qContigStart, qContigEnd,
                            score, strand, ".",
                            "Name=SAM_%s,ID=%s-%s;Target=%s %d %d" %
                            (tContig, qContig, tContig, tContig,
                             read.pos, read.aend)))
    if(totalCount % 100000 == 0):
        sys.stderr.write(".")

Last edited by gringer; 11-02-2011 at 12:17 PM. Reason: fixed up GFF3 parsing errors, added Name attribute to avoid contig clashes
gringer is offline   Reply With Quote
Old 11-02-2011, 12:14 PM   #4
ykdang
Junior Member
 
Location: Dallas

Join Date: Oct 2011
Posts: 7
Default

Thanks a lot, that's could be very useful
ykdang 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 06:55 PM.


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