[reut]

you can run cufflinks with the .bam file

Thanks chadn737,
In order to run Cufflink in default, is it I must include or exclude "-h" option?
eg.

samtools view input.bam > output.sam

Thanks again.

You don't have to convert the accepted_hits.bam to .sam for cufflinks, it works with the bam file as well.
(which is better, since the bam file is compressed and therefore a lot smaller than the sam file)

[edge]

Thanks chadn737,
In order to run Cufflink in default, is it I must include or exclude "-h" option?
eg.

samtools view input.bam > output.sam

Thanks again.

[chadn737]

Its fairly simple to convert bam to sam using samtools.

$ samtools view -h -o accepted_hits.sam accepted_hits.bam

[edge]

Hi telos,

Do you know that how to specify Tophat produce accepted_hits.sam?
After I run Tophat, why it only generate accepted_hits.bam
Thanks for advice.

[edge]

Hi telos,

Do you know that how to specify Tophat produce accepted_hits.sam?
After I run Tophat, why it only generate accepted_hits.bam
Thanks for advice.

[telos]

OK, fair enough.. I encountered the problem when comparing the SAM output with the GFF file from your script. Nothing a regexp can't solve, but it would be nice nevertheless if the file produced by your script were entirely consistent with the TopHat SAM output.

[genec]

Yeah, the MT/M thing is always an issue. Both MT and M will work, so there's not one that's right, you just have to be consistent from the beginning.

Gene

[telos]

MT -> chrM

You've omitted changing MT in the Ensembl GTF not to chrMT but to chrM for compatibility with TopHat.

[genec]

Bug fix

That was a good catch, Michelle. I'm attaching a fixed gtf_to_gff.pl. The previous version dropped the very last gene in the gtf file.

Gene

[mdimon]

thank you! (and a little bug?)

Thanks for the script! The C. elegans version is great for other GTF files downloaded from UCSC also.

I did notice what appears to be a little bug:
push @trs, [@exons];
should be added before the final
process(@trs);

(I am not a perl expert, I'm more of a python type, so I may be wrong, but until I added this line the last record from the GTF file didn't get printed to the GFF3 file.)

-- Michelle

[seqfast]

thank you!

Thanks very much, this works well. I had something similar but was getting hung up in the details. much appreciate people making these most useful scripts available - Thanks Gene,

-sf

[genec]

See the attached updated script. I modified it to work with your C elegans file. I believe it works, but give the output a good look to make sure that everything is processed correctly.

Gene

[seqfast]

script looks great, need help for c elegans

Thanks for the script, looks great and works well for the human gtf. I'm working on c.elegans gtf files (from ensembl), and the ENSG* strings aren't there ... i'm not a regex expert and figured I'd ask if it was an easy fix to use the c.elegans gtf files. I like this script for it's simplicity, I could use the other one mentioned in this thread if need be. Here is a snippet, i've also attached it in case of formatting issues. Thanks!

-sf

I snoRNA exon 3747 3909 . - . gene_id "Y74C9A.6"; transcript_id "Y74C9A.6"; exon_number "1"; gene_name "Y74C9A.6"; transcript_name "NR_001477.2";
I protein_coding exon 10095 10232 . - . gene_id "Y74C9A.3"; transcript_id "Y74C9A.3.1"; exon_number "1"; gene_name "Y74C9A.3"; transcript_name "Y74C9A.3.1";
I protein_coding CDS 10095 10148 . - 0 gene_id "Y74C9A.3"; transcript_id "Y74C9A.3.1"; exon_number "1"; gene_name "Y74C9A.3"; transcript_name "Y74C9A.3.1"; protein_id "Y74C9A.3.1";
I protein_coding start_codon 10146 10148 . - 0 gene_id "Y74C9A.3"; transcript_id "Y74C9A.3.1"; exon_number "1"; gene_name "Y74C9A.3"; transcript_name "Y74C9A.3.1";
I protein_coding exon 9727 9846 . - . gene_id "Y74C9A.3"; transcript_id "Y74C9A.3.1"; exon_number "2"; gene_name "Y74C9A.3"; transcript_name "Y74C9A.3.1";
I protein_coding CDS 9727 9846 . - 0 gene_id "Y74C9A.3"; transcript_id "Y74C9A.3.1"; exon_number "2"; gene_name "Y74C9A.3"; transcript_name "Y74C9A.3.1"; protein_id "Y74C9A.3.1";
I protein_coding exon 6037 6327 . - . gene_id "Y74C9A.3"; transcript_id "Y74C9A.3.1"; exon_number "3"; gene_name "Y74C9A.3"; transcript_name "Y74C9A.3.1";
I protein_coding CDS 6037 6327 . - 0 gene_id "Y74C9A.3"; transcript_id "Y74C9A.3.1"; exon_number "3"; gene_name "Y74C9A.3"; transcript_name "Y74C9A.3.1"; protein_id "Y74C9A.3.1";
I protein_coding exon 5195 5296 . - . gene_id "Y74C9A.3"; transcript_id "Y74C9A.3.1"; exon_number "4"; gene_name "Y74C9A.3"; transcript_name "Y74C9A.3.1";
I protein_coding CDS 5195 5296 . - 0 gene_id "Y74C9A.3"; transcript_id "Y74C9A.3.1"; exon_number "4"; gene_name "Y74C9A.3"; transcript_name "Y74C9A.3.1"; protein_id "Y74C9A.3.1";
I protein_coding exon 4124 4358 . - . gene_id "Y74C9A.3"; transcript_id "Y74C9A.3.1"; exon_number "5"; gene_name "Y74C9A.3"; transcript_name "Y74C9A.3.1";
I protein_coding CDS 4224 4358 . - 0 gene_id "Y74C9A.3"; transcript_id "Y74C9A.3.1"; exon_number "5"; gene_name "Y74C9A.3"; transcript_name "Y74C9A.3.1"; protein_id "Y74C9A.3.1";

[HTS]

I see. Thanks for the explanation! The reason gtf2gff3 doesn't work for you is probably because you forgot to convert chromosome names in the Ensembl convention to the UCSC convention? I forgot that I also wrote a small script to do that (among other things to filter the downloaded GTF file to suit my needs) before running gtf2gff3 (with the default configuration). I guess the real difference is that gtf2gff3 doesn't assume any particular ordering of the lines so it loads everything into memory and tries to figure out appropriate gene models from there. Since Ensmbl GTF files do group things according to genes/transcripts, it is good to explore that property.

[genec]

Yes, I had tried that gtf2gff3 script, but it wasn't working right for me. Maybe I didn't configure it correctly.

The script I posted has trivial memory requirements since it only holds one gene's worth of data in memory at once. All the exons for a gene are assumed to be located together in the gtf file, which seems to hold true for the Ensembl file. This script won't work for non-Ensembl gtf files without modification.

Gene