SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
One gene dominates RNAseq lynchde RNA Sequencing 1 05-14-2015 03:40 PM
Calculating total number of exons in the human copy of a gene rifah23 General 7 10-27-2014 02:21 PM
Getting gene coordinate not exons adrian Bioinformatics 3 08-07-2014 06:16 AM
Using DEXSeq normalized data to cluster candidate gene's exons? mwn84 Bioinformatics 0 10-01-2013 07:34 AM

Reply
 
Thread Tools
Old 03-11-2017, 07:36 AM   #1
sindrle
Senior Member
 
Location: Norway

Join Date: Aug 2013
Posts: 266
Default Analysing only exons from ONE gene using RNAseq

I know that DEXseq and/or Cuffdiff etc. may be used of exon-levels analyses.

However, I am only interested in ONE gene, because a novel exon is discovered in this gene, by other researches. There are one "old" known transcript using exon-1a and the rest of the gene, and this "new" transcript using exon-1b and the rest if the gene.

Now, I want to see if there is a different usage of exon-1a vs. exon-1b between lean and obese people.

My approach is:
1. Make a new custom GTF-annotation file with the novel exon. This GTF-file thus contains info about the "old" and the "new" transcript. Since it does not exists in any known repertoires yet.
2. Using tophat2, only provide chromosome 7 (where the gene is located) and this GTF as annotations.
3. Count the reads uniquely mapped exon-1a and exon-1b. And also sum all reads on the gene.
4. Normalize the exon expression to total gene expression and perhaps to library size?
5. Simply use linear regression i.e. exon~sum.of.reads+library.size+obese

Is this "valid"? Am I overlooking something major? I know I ignore some dispersion estimates etc. and that I ofc could run the whole DEXseq pipeline. But I want to save time and get intuitive answers. DEXseq is hard to understand and hard to explain to coworkers.

Thanks!

Last edited by sindrle; 03-11-2017 at 08:01 AM.
sindrle is offline   Reply With Quote
Old 03-12-2017, 03:46 AM   #2
gringer
David Eccles (gringer)
 
Location: Wellington, New Zealand

Join Date: May 2011
Posts: 836
Default

If you're looking at a single gene and have the full transcript sequences, why not just map to the transcripts and skip the complicated intron/exon searching?

If you're using Tophat2 for exon searches, you can switch to mapping to the transcripts with Bowtie2 without much change in the mapping pipeline.
gringer is offline   Reply With Quote
Old 03-12-2017, 06:57 AM   #3
sindrle
Senior Member
 
Location: Norway

Join Date: Aug 2013
Posts: 266
Default

Thats interesting!

Yes, I have the full transcript sequences. What do you think is the pros and cons between these approaches?
sindrle is offline   Reply With Quote
Old 03-12-2017, 02:24 PM   #4
gringer
David Eccles (gringer)
 
Location: Wellington, New Zealand

Join Date: May 2011
Posts: 836
Default

If the gene (or a subsequence of the gene) is not genome-unique, then non-specific matches to the gene of interest may happen.

I can't really think of any other cons, except for the whole throwing out 99.9% of your data thing, but that may be an advantage in this case because the needle is in the remainder.
gringer is offline   Reply With Quote
Old 03-13-2017, 09:03 AM   #5
sindrle
Senior Member
 
Location: Norway

Join Date: Aug 2013
Posts: 266
Default

Ok, so I got the transcript-sequences for the two transcripts of interest. Can you please explain how I should use them as reference for Bowtie2 (file format etc.)? An example with code would be very nice!

I probably can figure it out myself, but I would greatly appreciate your input, also to check that Im doing it correctly.

Thank you!

Last edited by sindrle; 03-13-2017 at 12:17 PM.
sindrle is offline   Reply With Quote
Old 03-13-2017, 11:18 AM   #6
gringer
David Eccles (gringer)
 
Location: Wellington, New Zealand

Join Date: May 2011
Posts: 836
Default

This should produce a count table for each transcript. It assumes that the most recent version of SAMtools is installed:

Code:
cat transcript1.fa transcript2.fa > allTranscripts_myGene.fa # combine transcripts into a single file
bowtie2-build allTranscripts_myGene.fa allTranscripts_myGene.fa # make bowtie2 index for transcripts
bowtie2 -x allTranscripts_myGene.fa -1 left_reads.fq.gz -2 right_reads.fq.gz | samtools sort -O BAM > reads_vs_myGene.bam # map reads to transcripts
samtools index reads_vs_myGene.bam # generate index file for mapped sequences
samtools idxstats reads_vs_myGene.bam # show mapping counts for transcripts
Note that any shared transcript areas will be randomly distributed among the two transcripts. A better idea of the differential expression may be obtained by looking at the reads in a pileup program like Tablet or IGV and only counting the reads for places where the transcripts differ. This will also give an indication of whether or not unspliced sequences are also in the RNASeq reads, which will also affect count results.
gringer is offline   Reply With Quote
Old 03-13-2017, 12:55 PM   #7
sindrle
Senior Member
 
Location: Norway

Join Date: Aug 2013
Posts: 266
Default

Thank you very much!

The two transcripts each have 2 unique sections (exon-a1 for one, and exon-1b for the other). I guess I can just define where these sequences are, in a GTF-file?, and use HTseq or similar to count the reads.
sindrle is offline   Reply With Quote
Old 03-13-2017, 01:01 PM   #8
gringer
David Eccles (gringer)
 
Location: Wellington, New Zealand

Join Date: May 2011
Posts: 836
Default

Yes, that will probably work. You need to be careful that you're not excluding reads that hit only partially to the unique regions (i.e. check the HTSeq-count options); any hits to the unique region of the transcripts should be counted as an isoform-specific hit.
gringer is offline   Reply With Quote
Old 03-13-2017, 01:12 PM   #9
sindrle
Senior Member
 
Location: Norway

Join Date: Aug 2013
Posts: 266
Default

Thats great!

Any input on how to normalise the counts? i.e. "counts/library size"?
sindrle is offline   Reply With Quote
Reply

Tags
cuffdiff, dexseq, exon, rnaseq

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 04:35 PM.


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