SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
DESeq2 without biol replicates sisterdot Bioinformatics 23 02-24-2019 08:01 AM
DESeq2 Simon Anders Bioinformatics 123 07-06-2015 01:45 AM
DESeq V DESeq2 sebastion RNA Sequencing 35 10-16-2014 06:04 AM
DESeq2-test a subset of genes quinne5 RNA Sequencing 2 10-16-2013 05:19 AM
Error Message in nbinomLRT in DESeq2 ToddB Bioinformatics 13 09-05-2013 06:22 AM

Reply
 
Thread Tools
Old 10-16-2013, 08:23 AM   #1
sindrle
Senior Member
 
Location: Norway

Join Date: Aug 2013
Posts: 266
Default DESeq2 SummarizedExperiment help

Hi!

I feel totally stupid, but I have 48 BAMs after alignment with Tophat2 and want to try DESeq2 for diff.exp. analysis.

I dont understand what "SummarizedExperiment" is and how to create it? I have googled to death, but now I have to ask here.

Thank you so much, this is a new field to me.
sindrle is offline   Reply With Quote
Old 10-16-2013, 09:35 AM   #2
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

Have you used htseq-count yet to get counts?
dpryan is offline   Reply With Quote
Old 10-16-2013, 11:14 AM   #3
sindrle
Senior Member
 
Location: Norway

Join Date: Aug 2013
Posts: 266
Default

Hi!
Do I have to? Im stuck here:

http://www.bioconductor.org/packages...zeOverlaps.pdf

Im trying the "2 A First Example", but I dont understand what "features <- GRanges" is (seqnames/IRanges).

Im also not quite sure how to type in "group_id": I have 2 different groups at two timepoints (as you already know).

Thank you very much!

Last edited by sindrle; 10-16-2013 at 11:30 AM.
sindrle is offline   Reply With Quote
Old 10-16-2013, 02:15 PM   #4
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

Well, you don't have to do anything you don't want to do :P However, it's generally a bit easier to just use a quick command line script than to try to deal with GenomicRanges. The later can do it, but has a much steeper learning curve and will just produce the same results (unless you're trying to do something unusual).

Regarding your GRanges question, seqnames would be things like "chr1" (UCSC chromosome naming scheme) or "10" (Ensembl naming scheme). IRanges are genomic bounds, so just the 5' and 3' boundaries of an exon/CDS/whatever. In general, you can get that stuff all made for you by reading in a GTF or GFF file (I've used import.gff from rtracklayer). The "group_id" in the example is more of a place-holder for gene_id later on. In essence, it's probably meant to be the actual feature level that's counted over (though the don't do that in the example), so you don't actually need it that example. In a real case, you'd read in your annotation file into a GRanges object and then
Code:
 split(GRangesObject, mcols(GRangesObject)$label_of_interest)
to get a GRangesList that you would then summarizeOverlaps on. Here, label_of_interest would probably be the gene_id field of the annotation file (though it could be something else, such as sprintf("%s:E%03d", mcols(blah)$gene_id, mcols(blah)$exon), which is what you'd do for DEXSeq counts).
dpryan is offline   Reply With Quote
Old 10-17-2013, 05:26 AM   #5
sindrle
Senior Member
 
Location: Norway

Join Date: Aug 2013
Posts: 266
Default

Thank you so much!

I installed HTseq and have started HTseq-count on all 48 BAMs.

From one BAM, is this a problem?:

no_feature 7984183
ambiguous 521840
too_low_aQual 0
not_aligned 0
alignment_not_unique 13730878


And also, can you compress the HTseq-count output? Its 16gb for each file..

Last edited by sindrle; 10-17-2013 at 05:46 AM.
sindrle is offline   Reply With Quote
Old 10-17-2013, 05:59 AM   #6
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

Ah, you're saving the wrong file, don't use -o. The command to execute is something like:
Code:
samtools view namesorted.bam | htseq-count -m intersection-nonempty -s no -a 10 - GTF_FILE > namesorted.counts
That file will be much much much smaller and is all that you need.
dpryan is offline   Reply With Quote
Old 10-17-2013, 06:01 AM   #7
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

I should add, "namesorted.bam" and "GTF_FILE" above are just place-holder names. You'll want to change "namesorted.counts" to something meaningful so you can keep track of which counts are associated with which sample. Finally, those options for htseq-count are just examples, you can change them to fit your needs.
dpryan is offline   Reply With Quote
Old 10-17-2013, 06:24 AM   #8
sindrle
Senior Member
 
Location: Norway

Join Date: Aug 2013
Posts: 266
Default

Great!
I aborted the run and updated the CMD as you wrote it.

I was using the "union" option, but changed it to "intersection-nonempty". But I dont actually know why.. :P

Thanks again!

Last edited by sindrle; 10-17-2013 at 06:45 AM.
sindrle is offline   Reply With Quote
Old 10-23-2013, 03:14 AM   #9
frymor
Senior Member
 
Location: Germany

Join Date: May 2010
Posts: 149
Default

Quote:
Originally Posted by sindrle View Post
Great!
I aborted the run and updated the CMD as you wrote it.

I was using the "union" option, but changed it to "intersection-nonempty". But I dont actually know why.. :P

Thanks again!
take a look here to see the difference between the union and the intersection-nonempty.

It basically tells htseq-count how to deal with reads matched to two genes, or what to do if the overlap isn't complete.
frymor is offline   Reply With Quote
Old 10-23-2013, 03:51 AM   #10
sindrle
Senior Member
 
Location: Norway

Join Date: Aug 2013
Posts: 266
Default

Thanks I have looked at that one, but what is the implications for differential expression?
sindrle is offline   Reply With Quote
Old 10-23-2013, 04:20 AM   #11
frymor
Senior Member
 
Location: Germany

Join Date: May 2010
Posts: 149
Default

As far as I understand it, depends on the option you take you'll be counting more or less reads.

If I am not mistaken, reads marked as ambiguous are not taken into consideration in the count of read numbers. So if you take union you'll have more ambiguous reads than the intersection options.

BTW, take a look at featureCounts from the sub read package. It works also very good with the data. It finds more results than htseq-counts and works faster.
frymor 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 01:42 AM.


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