SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Survey: RNA-Seq analysis for Differential Gene/Transcript Expression bodhisattvax Bioinformatics 14 06-12-2013 10:06 AM
differential gene expression analysis between Different strains by RNA-seq qqtwee Bioinformatics 3 07-30-2012 05:19 AM
Differential gene expression analysis with bioreplicates using EdgeR/DESeq iceage Bioinformatics 1 07-07-2012 05:58 PM
Differential gene expression analysis without reference cerebralrust Bioinformatics 7 05-04-2012 03:57 AM

Reply
 
Thread Tools
Old 05-17-2013, 08:55 AM   #1
colaneri
Member
 
Location: Durham

Join Date: Jul 2012
Posts: 30
Default Differential gene expression analysis

I'm trying to establish the best gene expression differential analysis for my purpose: 2 genotypes, 2 experimental situations, 3 biological replicates, 25 million reads per sample (sequenced RNA-seq libraries).
Now I'm using tophat-cufflinks and following the protocol published by Trapnell:
Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks

I'm working with a well annotated model organism "Arabidopsis thaliana"

I have two goals:

First look for diff expression in the already annotated transcriptome: TAIR10
Second: I'm interested in possibility that previously NO annotated genes are differentially expressed between the two genotypes in one of the different experimental conditions.
I have a protocol in mind BUT I will like to be advised for the expertize of this community:

Remember: two genotypes, Two experimental confitions, triplicates MEANS 12 INDEPENDENT LIBRARIES (25 millions reads each)

FIRST PROTOCOL FOR DIFFERENTIAL EXPRESSION:
1) Tophat for each library
2) Merge all the libraries in a single cufflinks (do I need to include the TAIR10.gtf?)
3) Use the final assembly of step two togheter with the 12 acepted hits files from step one in cuffdiff.
4) Use cuffcompare to identify locations of new genes.

How can I automatically extract all those new genes thar are also differentially expressed?

I will appreciate feedbacks for the protocol I have in mind and answer to my questions
colaneri is offline   Reply With Quote
Old 05-17-2013, 10:16 AM   #2
kmcarr
Senior Member
 
Location: USA, Midwest

Join Date: May 2008
Posts: 1,170
Default

What you are describing is exactly what the Cuffmerge program from the Cufflinks suite does. It takes the .gtf files from a set of cufflinks runs and merges them together into a single, non-redundant set of transfrags. You may optionally provides a GTF file contain a set of already known transcripts and the output mapping file will classify the transfrags as known, novel etc. See the class code documentation in the manual.

Now if I may offer a perspective as one who has done a lot of RNA-Seq analysis in Arabidopsis, the A. thaliana genome has been analyzed and annotated to death. The odds of finding a new gene are very, very slim an probably not worth spending time working on unless you have time to waste. I can guarantee you that you will find groups of reads mapping to intergenic space, but will also bet large amounts of money that these will not be new genes. They most likely arise from mis-mapping or spurious transcriptional events.
kmcarr is offline   Reply With Quote
Old 05-20-2013, 11:33 AM   #3
colaneri
Member
 
Location: Durham

Join Date: Jul 2012
Posts: 30
Default no gene/ transcript discovery

Assuming that as you said, the Arabidopsis genome is very well annotated, do I need to run cufflinks?
It is better in any way to perform the analysis just combining tophat, the referece trancriptome and cuffdiff?
Do you suggest to run TopHat with "no gene/transcript discovery"?




Quote:
Originally Posted by kmcarr View Post
What you are describing is exactly what the Cuffmerge program from the Cufflinks suite does. It takes the .gtf files from a set of cufflinks runs and merges them together into a single, non-redundant set of transfrags. You may optionally provides a GTF file contain a set of already known transcripts and the output mapping file will classify the transfrags as known, novel etc. See the class code documentation in the manual.

Now if I may offer a perspective as one who has done a lot of RNA-Seq analysis in Arabidopsis, the A. thaliana genome has been analyzed and annotated to death. The odds of finding a new gene are very, very slim an probably not worth spending time working on unless you have time to waste. I can guarantee you that you will find groups of reads mapping to intergenic space, but will also bet large amounts of money that these will not be new genes. They most likely arise from mis-mapping or spurious transcriptional events.
colaneri is offline   Reply With Quote
Old 05-20-2013, 11:54 AM   #4
kmcarr
Senior Member
 
Location: USA, Midwest

Join Date: May 2008
Posts: 1,170
Default

Quote:
Originally Posted by colaneri View Post
Assuming that as you said, the Arabidopsis genome is very well annotated, do I need to run cufflinks?
No
Quote:
It is better in any way to perform the analysis just combining tophat, the referece trancriptome and cuffdiff?
Yes. It's faster because you are skipping an unneccessary step, and the IDs used for cuffdiff analysis will be the normal TAIR AT IDs instead of cufflinks transfrag IDs (XLOCs) which you would then need to correlate to their TAIR IDs.
Quote:
Do you suggest to run TopHat with "no gene/transcript discovery"?
That's what I normally do.
kmcarr is offline   Reply With Quote
Old 05-20-2013, 03:34 PM   #5
kmcarr
Senior Member
 
Location: USA, Midwest

Join Date: May 2008
Posts: 1,170
Default

Quote:
Originally Posted by colaneri View Post
I want to use tophat in galaxy with the parameter --no-novel-juncs genome

how can I implement the parameter?
Sorry, I don't use Galaxy so can't help you there.
kmcarr is offline   Reply With Quote
Old 05-23-2013, 08:41 AM   #6
dnusol
Senior Member
 
Location: Spain

Join Date: Jul 2009
Posts: 133
Default

I think this was solved elsewhere

http://seqanswers.com/forums/showthread.php?t=16965


HTH
dnusol is offline   Reply With Quote
Old 05-23-2013, 02:51 PM   #7
colaneri
Member
 
Location: Durham

Join Date: Jul 2012
Posts: 30
Default Cuffdiff do not performed promoter preference test

I have run cuffdiff 2 using this command

cuffdiff -p 4 -c 4 --no-update-check /proj/seq/data/TAIR10_Ensembl/Annotation/Archives/archive-2013-03-06-09-54-25/Genes/genes.gtf -o cuffdiff_results ./C_ctrl_rep1_Trim37.tophat_out/accepted_hits.bam ./C_ABA_rep1_Trim40.tophat_out/accepted_hits.bam

even though the work was successfully completed I received the below stated output. My question is why were not test performed for promoter preference or splicing? It is the default option?




Performed 27529 isoform-level transcription difference tests
Performed 25428 tss-level transcription difference tests
Performed 21976 gene-level transcription difference tests
Performed 24788 CDS-level transcription difference tests
Performed 0 splicing tests
Performed 0 promoter preference tests
Performing 0 relative CDS output tests

Writing isoform-level FPKM tracking
Writing TSS group-level FPKM tracking
Writing gene-level FPKM tracking
Writing CDS-level FPKM tracking
Writing isoform-level count tracking
Writing TSS group-level count tracking
Writing gene-level count tracking
Writing CDS-level count tracking
Writing isoform-level read group tracking
Writing TSS group-level read group tracking
Writing gene-level read group tracking
Writing CDS-level read group tracking
Writing read group info
Writing run info
colaneri is offline   Reply With Quote
Old 05-29-2013, 11:24 AM   #8
colaneri
Member
 
Location: Durham

Join Date: Jul 2012
Posts: 30
Default concatenating files before tophat?

Hi
I have a RNA-seq library that has been sequenced multiple times, then I have four fastq files.
Do I need to concatenate them before alignment in tophat?
Can i just list the four files at the end of the tophat command like that?

If my files are fastq1, fastq2, fastq3 and fastq4,

and I do:

tophat -p 4 --segment-length 20 --no-novel-juncs -G /proj/seq/data/TAIR10_Ensembl/Annotation/Archives/archive-2013-03-06-09-54-25/Genes/genes.gtf -o C_ctrl_rep1_THout_6 /proj/seq/data/TAIR10_Ensembl/Sequence/Bowtie2Index/genome fastq1 fastq2 fastq3 fastq4
colaneri is offline   Reply With Quote
Old 05-29-2013, 11:36 AM   #9
colaneri
Member
 
Location: Durham

Join Date: Jul 2012
Posts: 30
Default selecting the approapiate range to trim

Most of my RNA-seq sequenced libraries look like that in the fastqc report (please see the images below)
Do I need to trim the first 10 bases?
Only the first 10 ones?
It is going to improve the results?

[IMG] [/IMG]
colaneri is offline   Reply With Quote
Old 06-10-2013, 02:07 PM   #10
colaneri
Member
 
Location: Durham

Join Date: Jul 2012
Posts: 30
Default Defining replicates and different conditions in Cuffdiff2

Hi peoples
I'm trying to make cuffdiff 2 to compare RNA-seq data from
2 different genotypes in two different conditions and I did 3 biological replicates for each genotype in each condition.
So I have 12 different libraries, I aligned them separately with tophat.
My problem is in running cuffdiff from the command line, I can not get it to work in the way I would like, and I do not know what I'm doing wrong. PLEASE SOME HELP IN HERE GUYS!!!

I did run cuffdiff with this command
cuffdiff -p 8 -c 20 --no-update-check /proj/seq/data/TAIR10_Ensembl/Annotation/Archives/archive-2013-03-06-09-54-25/Genes/genes.gtf -o cuffdiff_ABA_whole_set_results_week \
./C_ctrl_rep1_abaexp.tophat_out/accepted_hits.bam, ./C_ctrl_rep2_abaexp.tophat_out/accepted_hits.bam, ./C_ctrl_rep3_abaexp.tophat_out/accepted_hits.bam \
./C_ABA_rep1_abaexp.tophat_out/accepted_hits.bam, ./C_ABA_rep2_abaexp.tophat_out/accepted_hits.bam, ./C_ABA_rep3_abaexp.tophat_out/accepted_hits.bam \
./B_ctrl_rep1_abaexp.tophat_out/accepted_hits.bam, ./B_ctrl_rep2_abaexp.tophat_out/accepted_hits.bam, ./B_ctrl_rep3_abaexp.tophat_out/accepted_hits.bam \
./B_ABA_rep1_abaexp.tophat_out/accepted_hits.bam, ./B_ABA_rep2_abaexp.tophat_out/accepted_hits.bam, ./B_ABA_rep3_abaexp.tophat_out/accepted_hits.bam


BUT THE RESULT IS THAT ALL THE FILES ARE COMPARED AGAINS THE OTHER, so all samples are considered different instead of 4 groups with triplicates

CAN SOME ONE TELL ME WHAT IS WRONG WITH MY COMMAND LINE?
colaneri is offline   Reply With Quote
Old 06-10-2013, 08:47 PM   #11
kmcarr
Senior Member
 
Location: USA, Midwest

Join Date: May 2008
Posts: 1,170
Default

Quote:
Originally Posted by colaneri View Post
Hi peoples
I'm trying to make cuffdiff 2 to compare RNA-seq data from
2 different genotypes in two different conditions and I did 3 biological replicates for each genotype in each condition.
So I have 12 different libraries, I aligned them separately with tophat.
My problem is in running cuffdiff from the command line, I can not get it to work in the way I would like, and I do not know what I'm doing wrong. PLEASE SOME HELP IN HERE GUYS!!!

I did run cuffdiff with this command
Code:
cuffdiff -p 8 -c 20 --no-update-check /proj/seq/data/TAIR10_Ensembl/Annotation/Archives/archive-2013-03-06-09-54-25/Genes/genes.gtf -o cuffdiff_ABA_whole_set_results_week \
./C_ctrl_rep1_abaexp.tophat_out/accepted_hits.bam, ./C_ctrl_rep2_abaexp.tophat_out/accepted_hits.bam, ./C_ctrl_rep3_abaexp.tophat_out/accepted_hits.bam \
./C_ABA_rep1_abaexp.tophat_out/accepted_hits.bam, ./C_ABA_rep2_abaexp.tophat_out/accepted_hits.bam, ./C_ABA_rep3_abaexp.tophat_out/accepted_hits.bam \
./B_ctrl_rep1_abaexp.tophat_out/accepted_hits.bam, ./B_ctrl_rep2_abaexp.tophat_out/accepted_hits.bam, ./B_ctrl_rep3_abaexp.tophat_out/accepted_hits.bam \
./B_ABA_rep1_abaexp.tophat_out/accepted_hits.bam, ./B_ABA_rep2_abaexp.tophat_out/accepted_hits.bam, ./B_ABA_rep3_abaexp.tophat_out/accepted_hits.bam

BUT THE RESULT IS THAT ALL THE FILES ARE COMPARED AGAINS THE OTHER, so all samples are considered different instead of 4 groups with triplicates

CAN SOME ONE TELL ME WHAT IS WRONG WITH MY COMMAND LINE?
You have spaces after the commas in your command line. The list of BAM files for your bio reps should be separated by commas WITHOUT SPACES, then spaces between the different condition groups. Since you put a space after every BAM file name cuffdiff interpreted them as twelve conditions.

BTW if you are posting blocks of command text or code please use the CODE tag formatting as I have done with your text above. It makes reading lines of code or output much easier and thus easier to spot the errors.
kmcarr is offline   Reply With Quote
Old 06-11-2013, 06:35 AM   #12
colaneri
Member
 
Location: Durham

Join Date: Jul 2012
Posts: 30
Default naming samples in cuffdiff

Than you very much KMCARR!

But the way, when I use the -L option to name the different samples,
do I also have to separate them with commas without spaces?
In term of the names, do I need to use exactly the same name that the one is pointing to the bam file? Or it is just the order of names after the -L option what it matters?
colaneri is offline   Reply With Quote
Old 06-11-2013, 01:11 PM   #13
sdriscoll
I like code
 
Location: San Diego, CA, USA

Join Date: Sep 2009
Posts: 438
Default

Quote:
Originally Posted by colaneri View Post
when I use the -L option to name the different samples, do I also have to separate them with commas without spaces?
Yes.

Quote:
Originally Posted by colaneri View Post
In term of the names, do I need to use exactly the same name that the one is pointing to the bam file? Or it is just the order of names after the -L option what it matters?
It's only the order that matters.
__________________
/* Shawn Driscoll, Gene Expression Laboratory, Pfaff
Salk Institute for Biological Studies, La Jolla, CA, USA */
sdriscoll is offline   Reply With Quote
Old 06-13-2013, 07:49 AM   #14
colaneri
Member
 
Location: Durham

Join Date: Jul 2012
Posts: 30
Default I do not understand this tophat error

I have a fastq file that I used to align sequences with tophat v 1.3 (from a galaxy server) and I have not problem, but when I use the same fastq file to align the sequences with tophat 2 in command line I get this error.

Can you please explain to me why and what does means?

This is the output: (error is highlighted in red at the bottom)

[2013-06-13 00:43:09] Checking for Bowtie
Bowtie version: 2.1.0.0
[2013-06-13 00:43:09] Checking for Samtools
Samtools version: 0.1.19.0
[2013-06-13 00:43:09] Checking for Bowtie index files
[2013-06-13 00:43:09] Checking for reference FASTA file
[2013-06-13 00:43:09] Generating SAM header for /proj/seq/data/TAIR10_Ensembl/Sequence/Bowtie2Index/genome
format: fastq
quality scale: phred33 (default)
[2013-06-13 00:43:12] Reading known junctions from GTF file
[2013-06-13 00:43:16] Preparing reads
[FAILED]
Error running 'prep_reads'
Error: qual length (19) differs from seq length (41) for fastq record !
colaneri is offline   Reply With Quote
Old 06-13-2013, 11:41 PM   #15
sdriscoll
I like code
 
Location: San Diego, CA, USA

Join Date: Sep 2009
Posts: 438
Default

What it means to me is that the qualities string and the read string for at least one of the reads in your fastq file are not the same length. It doesn't explain why it worked on galaxy though.
__________________
/* Shawn Driscoll, Gene Expression Laboratory, Pfaff
Salk Institute for Biological Studies, La Jolla, CA, USA */
sdriscoll is offline   Reply With Quote
Old 06-14-2013, 06:37 AM   #16
colaneri
Member
 
Location: Durham

Join Date: Jul 2012
Posts: 30
Default quality string and read string

Quote:
Originally Posted by sdriscoll View Post
What it means to me is that the qualities string and the read string for at least one of the reads in your fastq file are not the same length. It doesn't explain why it worked on galaxy though.
Thanks for your answer, I found that someway the file was corrupted when copying between my laptop to the server. I rewrite the file and now is working.
colaneri is offline   Reply With Quote
Reply

Tags
cuffcompare, cuffdiff, cufflinks, cuffmerge

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 12:12 PM.


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