SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
GC content of assembled transcripts DrYak Bioinformatics 2 07-04-2016 09:54 PM
Is blastx reliable enough to annotate transcripts from RNA-seq? mlacencio Bioinformatics 5 02-02-2015 11:34 AM
visualize alternative splicing from cufflink output gtf file lzu Bioinformatics 3 12-04-2014 11:10 AM
how to run tophat without gtf file and annotate the gene mehtaaditya Bioinformatics 0 03-19-2013 08:53 AM
some odd out put in combined GTF RNA-seq Cufflink honey Bioinformatics 0 07-30-2011 08:37 AM

Reply
 
Thread Tools
Old 01-31-2017, 06:52 PM   #1
capricy
Senior Member
 
Location: 63130

Join Date: Apr 2012
Posts: 125
Default annotate cufflink assembled transcripts with reference gtf

Hello, there,

1. I did genome-guided de novo transcripts assembly for my RNAseq data using cufflinks. Here .sam file is from STAR mapping

cufflinks -p 8 /mapping/mapped.sam

2. I then merged the resultant gtf files from the same tissue to have merged.gtf without including reference.gtf

cuffmerge -p 8 gtf.filelist.DeNovo

3. I tried to find the closest gene id for those de novo assembled transcripts

cuffcompare merged.gtf -r reference.gtf

What I've found is that none of my de novo assembled transcripts are mapped to the reference gtf even though some introns are apparently identical between the merged.gtf and reference.gtf

for example:

from the cufflinks merged.gtf, I have

more XLOC_005458.gtf
chr2 Cufflinks exon 25289899 25290661 . . . gene_id "XLOC_005458"; transcript_id "TCONS_00010739"; exon_number "1"; oId "CUFF.5451.1"; tss_i
d "TSS7438";
chr2 Cufflinks exon 25290738 25290883 . . . gene_id "XLOC_005458"; transcript_id "TCONS_00010739"; exon_number "2"; oId "CUFF.5451.1"; tss_i
d "TSS7438";
chr2 Cufflinks exon 25290976 25291190 . . . gene_id "XLOC_005458"; transcript_id "TCONS_00010739"; exon_number "3"; oId "CUFF.5451.1"; tss_i
d "TSS7438";
chr2 Cufflinks exon 25289938 25290082 . . . gene_id "XLOC_005458"; transcript_id "TCONS_00010740"; exon_number "1"; oId "CUFF.5451.2"; tss_i
d "TSS7438";
chr2 Cufflinks exon 25290388 25291177 . . . gene_id "XLOC_005458"; transcript_id "TCONS_00010740"; exon_number "2"; oId "CUFF.5451.2"; tss_i
d "TSS7438";

from the reference.gtf, I have:
2 ensembl_havana CDS 25289989 25290661 . + 0 ccds_id "CCDS15763"; exon_number "1"; gene_biotype "protein_coding"; gene_id "ENSMUSG00000026961"; gene_name "Lrrc26"; gene_source "ensembl_havana"; gene_version "6"; havana_gene "OTTMUSG00000011934"; havana_gene_version "1"; havana_transcript "OTTMUST00000028197"; havana_transcript_version "1"; p_id "P45943"; protein_id "ENSMUSP00000028337"; protein_version "6"; tag "basic"; transcript_biotype "protein_coding"; transcript_id "ENSMUST00000028337"; transcript_name "Lrrc26-001"; transcript_source "ensembl_havana"; transcript_support_level "1"; transcript_version "6"; tss_id "TSS86428";

2 ensembl_havana CDS 25290738 25291057 . + 2 ccds_id "CCDS15763"; exon_number "2"; gene_biotype "protein_coding"; gene_id "ENSMUSG00000026961"; gene_name "Lrrc26"; gene_source "ensembl_havana"; gene_version "6"; havana_gene "OTTMUSG00000011934"; havana_gene_version "1"; havana_transcript "OTTMUST00000028197"; havana_transcript_version "1"; p_id "P45943"; protein_id "ENSMUSP00000028337"; protein_version "6"; tag "basic"; transcript_biotype "protein_coding"; transcript_id "ENSMUST00000028337"; transcript_name "Lrrc26-001"; transcript_source "ensembl_havana"; transcript_support_level "1"; transcript_version "6"; tss_id "TSS86428";

Apparently the same intron (25290661 .. 25290738) exists in both the de novo assemble transcript and the reference. So my question is why the XLOC_005458 from cufflinks output is not mapped to the Lrrc26 in reference.gtf even though they share the same gene region?

Thanks for any inputs!

C.

Last edited by capricy; 02-01-2017 at 07:15 AM.
capricy is offline   Reply With Quote
Old 01-31-2017, 09:21 PM   #2
shunyip
Member
 
Location: New York

Join Date: Oct 2013
Posts: 20
Default

RNA molecules can suffer from degradation. However, introns are identified by splice junctions and are often in the middle of the RNA reads. So, it is more likely for introns to be identified correctly. If you want all genes to be mapped very similar to the reference, you might need higher sequencing depth and/or higher quality data.
shunyip is offline   Reply With Quote
Old 02-01-2017, 03:19 AM   #3
capricy
Senior Member
 
Location: 63130

Join Date: Apr 2012
Posts: 125
Default

Then what is the easy way to annotate those assembled transcripts? I meant, I would like to find the closest reference gene IDs for the transcripts.

Thanks.

C.
capricy is offline   Reply With Quote
Old 02-01-2017, 06:37 AM   #4
shunyip
Member
 
Location: New York

Join Date: Oct 2013
Posts: 20
Default

Hi Capricy,

You can supply your gene annotation (reference.gtf) to Cufflinks during assembly, using the -g argument.
Or you can use bedtools intersect to overlap and combine your merged.gtf and reference.gtf. Here is its document. You need to convert the gtf files into bed files for this method.
http://bedtools.readthedocs.io/en/la...intersect.html

I hope this helps,
shunyip is offline   Reply With Quote
Old 02-01-2017, 07:10 AM   #5
capricy
Senior Member
 
Location: 63130

Join Date: Apr 2012
Posts: 125
Default

According to the cuffcompare document, if I use -r <reference.gtf>, the output should be able to identify the overlapped transfrags. But it did not in my case.

Just wonder if there is something wrong with my steps?

C.
capricy is offline   Reply With Quote
Old 02-01-2017, 07:22 AM   #6
capricy
Senior Member
 
Location: 63130

Join Date: Apr 2012
Posts: 125
Default

I didn't use -g since I only would like to see the de novo assembled transfrags.
capricy is offline   Reply With Quote
Old 02-01-2017, 07:52 AM   #7
shunyip
Member
 
Location: New York

Join Date: Oct 2013
Posts: 20
Default

Then, it would seem that an easy way for you is to use bedtools.

You can convert a gtf file to bed file using:
Code:
cut -f 1,4,5,9 yourfile.gtf > yourfile.bed
This extracts the 1st, 4th, 5th and 9th columns from the gtf files and write them to a new file.

Then, you can use bedtools intersect to overlap the two files.
It seems that the -loj and -wao arguments suit your case well. You can take a look.
shunyip 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 02:49 PM.


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