Cufflinks gffread utility
1
0
Entering edit mode
5.3 years ago
Thomas ▴ 110

Hello,

This is a question regarding the following utility:

http://cole-trapnell-lab.github.io/cufflinks/file_formats/#the-gffread-utility

If I use gffread for genomic features labelled as existing on the negative strand, will gffread find the reverse complement automatically when extracting the sequence from the fasta file, or will I have to implement an additional step to get that?

Any info on this would be greatly appreciated.

cufflinks gffread gff reverse complement • 4.7k views
ADD COMMENT
0
Entering edit mode

Hello! Any ideas why I'm getting a seg fault when I run this?

gffread BF_annotation.gtf -W -O -E -L -F -w BF_transcripts.fa -g BF_genome.fa
GFF Warning: merging overlapping/adjacent feature segment exon (8436614-8437382) with exon (8436535-8436613) for GFF ID KDM4A on scaffold_1
GFF Warning: merging overlapping/adjacent feature segment exon (319524-319638) with exon (319006-319523) for GFF ID ZFYVE28 on scaffold_199
GFF Warning: merging overlapping/adjacent feature segment exon (1352880-1353047) with exon (1351726-1352879) for GFF ID FLAD1 on scaffold_218
Segmentation fault (core dumped)

and when I run without certain arguments: gffread -w BF_transcripts.fa -W -O -E -g ./BF_genome.fa BF_annotation.gtf

it outputs warnings:

GFF Warning: merging overlapping/adjacent feature segment exon (8436614-8437382) with exon (8436535-8436613) for GFF ID KDM4A on scaffold_1
GFF Warning: merging overlapping/adjacent feature segment exon (319524-319638) with exon (319006-319523) for GFF ID ZFYVE28 on scaffold_199
GFF Warning: merging overlapping/adjacent feature segment exon (1352880-1353047) with exon (1351726-1352879) for GFF ID FLAD1 on scaffold_218
Error (GFaSeqGet): subsequence cannot be larger than 1560
Error getting subseq for BLEC (1..1561)!

I hadn't seen this post before and made a related post: C: Get sequences for all genes in a GTF file

Thanks!

ADD REPLY
0
Entering edit mode

I cleaned the .gtf from some genes (4 of them) and it run without warnings, but I'm not clear this is how I should have proceeded.

ADD REPLY
0
Entering edit mode

In the warning messages, it is just saying that some of your gene exon co-ordinates overlap. For example:

GENE1
___________________________________________
TRANSCRIPTION --> 
                                 _______________________________________________________
                                 TRANSCRIPTION -->               
                                                                                   GENE2

You can see that there is some overlap in the middle. This is the situation that you have on 3 of your scaffolds (1, 199, 218).

This is not necessarily a problem. Even in the human transcriptome CDS FASTA file from GENCODE, many genes have the same sequences because they also overlap. It may be a problem depending on what are your next analysis steps.

ADD REPLY
5
Entering edit mode
3.9 years ago

Sorry the the late response.

Yes, it certainly will read the reverse complement and take strand-specific information into account:

cat mergedAll.gtf 
BacterialGenome.1   Cufflinks   exon    2   500 .   +   .   gene_id "XLOC_000001"; transcript_id "TCONS_00000001"; exon_number "1"; oId "CUFF.10.1"; tss_id "TSS1";
BacterialGenome.1   Cufflinks   exon    2   500 .   -   .   gene_id "XLOC_000628"; transcript_id "TCONS_00000630"; exon_number "1"; oId "CUFF.1.1"; tss_id "TSS630";

.

gffread -w test.fasta -W -O -E -L -F -g BacterialGenomeReference.fasta mergedAll.gtf

.

cat test.fasta 
>TCONS_00000001 gene=XLOC_000001 loc:BacterialGenome.1|2-500|+ exons:2-500 segs:1-499 oId=CUFF.10.1 tss_id=TSS1
CATCACAGACACTACGCATAAGGGGCTATCACCCTCTATGGCCGGACTTTCCATTCCGTTCTGCTTCTTC
TGCAACAATCACGGGCCTGTTCCGCGTTCGCTCGCCACTACTAGCGGAATCTCAATTGATGTCTTTTCCT
CCGGGTACTTAGATGTTTCAGTTCCCCGGGTTCGCCTCATGCCCCTATGTATTCAGAACATGATACCCAT
CGCTGGGTGGGTTGCCCCATTCAGATATCCACGGATCAAAGCCTGCTCGCGGCTCCCCATGGCTTTTCGC
AGCGTGCCACGTCTTTCATCGCCTCCTGGTGCCAAGGCATCCACCGAATGCCCTTATCGCGCTCATTCAC
CACACATGCACAGGAGCCATCCACCCTGGGGCAGACAGTCCCGCACATAAGAAAGTGCAGTTCATCTCAC
GACCACTCTATTCTCTTTACGTCGTCCATGTTCGCTTACGCCACATCGCACAGCACGCAATACCCCGGTT
CCCCGGACC
>TCONS_00000630 gene=XLOC_000628 loc:BacterialGenome.1|2-500|- exons:2-500 segs:1-499 oId=CUFF.1.1 tss_id=TSS630
GGTCCGGGGAACCGGGGTATTGCGTGCTGTGCGATGTGGCGTAAGCGAACATGGACGACGTAAAGAGAAT
AGAGTGGTCGTGAGATGAACTGCACTTTCTTATGTGCGGGACTGTCTGCCCCAGGGTGGATGGCTCCTGT
GCATGTGTGGTGAATGAGCGCGATAAGGGCATTCGGTGGATGCCTTGGCACCAGGAGGCGATGAAAGACG
TGGCACGCTGCGAAAAGCCATGGGGAGCCGCGAGCAGGCTTTGATCCGTGGATATCTGAATGGGGCAACC
CACCCAGCGATGGGTATCATGTTCTGAATACATAGGGGCATGAGGCGAACCCGGGGAACTGAAACATCTA
AGTACCCGGAGGAAAAGACATCAATTGAGATTCCGCTAGTAGTGGCGAGCGAACGCGGAACAGGCCCGTG
ATTGTTGCAGAAGAAGCAGAACGGAATGGAAAGTCCGGCCATAGAGGGTGATAGCCCCTTATGCGTAGTG
TCTGTGATG
ADD COMMENT
0
Entering edit mode

In gffread, if I need only TCONS_00033653 and gene_id "XLOC_028018 ID. what is command I should use?

Actually, I used gffread -w ECvsEL_transcript.fasta -g v.corymbosum_GDV_reftransV1.fasta EC_vs_EL.gtf

But I get only TCONS_00033653 ID

Thank you so much kanjana

ADD REPLY

Login before adding your answer.

Traffic: 2142 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6