[piaopiaogg]

--------------------------------------------------------------------------------

Hi Matt,

Really really useful. Thank you

[bioinfosm]

Thanks for the useful resource!

[ttnguyen]

Excellent guide! Looking forward for ChIP-seq analysis also...

[ge_SF]

Very nice guide!

I do have one point- the maketranscriptDb function in GenomicFeatures seems like it is only useful for model organisms (at least, I am not smart enough to make it work for my organism.)

I have managed to use the girafe package to do something similar using just the gff. I can write up something if people think it would be useful.

[MDY]

Gunzip - I unfortunately know of no such guide. Sorry.

JueFish - This command uses regular expressions to rename chromosomes. There are many tutorials on the internet on how to use regular expressions (see http://www.regular-expressions.info/ for example) and the R help for gsub has some useful info too. The point is basically that I have chromosomes named like this:
10.1-129993255
11.1-121843856
12.1-121257530
13.1-120284312
14.1-125194864
15.1-103494974
16.1-98319150
17.1-95272651
18.1-90772031
19.1-61342430
2.1-181748087
3.1-159599783
4.1-155630120
5.1-152537259
6.1-149517037
7.1-152524553
8.1-131738871
9.1-124076172
MT.1-16299
X.1-166650296
Y.1-15902555
1.1-197195432
But my annotation files have them named like this:
chr1
chr10
chr11
chr12
chr13
chr14
chr15
chr16
chr17
chr18
chr19
chr2
chr3
chr4
chr5
chr6
chr7
chr8
chr9
chrM
chrX
chrY
The regular expression just converts the names. You don't have to use a regular expression, it's just a time saving tool.


plassaaw - cufflinks/cuffdiff does not do the same thing as an edgeR analysis. Both are useful in different circumstances (see my response to Azazel's post).

[plassaaw]

Hi,
This seems to be a very good walk through. One thing though, why did you use the R package to check for differential expression as opposed to using cufflinks? In my experience, I have found that cufflinks is very easy to use, since you are producing sam files from bowtie and both programs are developed by UMD.

Just wondering what you thought!

Thanks!

Andrew

[m.nyine]

Hi,

Thanks a lot

This is really helpful for a fresher in bioinformatics like me.

I hope our seniors will take MDY's idea of developing this further.
I know many are working hard to develop such material but there is limited access to it.

Please keep us informed of any developments.

Nyine

[JueFish]

Been checking out the tutorial that was posted on this thread. Can anyone comment what the command:

new_read_chr_names=gsub("(.*)[T]*\\..*","chr\\1",rname(reads))

is doing? I get the making of a new list of chromosome names and that the eaxmple uses gsub to do the substitutions, but I don't under stand what's going on in the first two fields of that command:

"(.*)[T]*\\..*"

and

"chr\\1"

In other words, I have no idea how to make sure that my chromosomal names will match up with ones in the genome (NCBI headers put in a bunch of noise) because of my naivete about syntax. Any thoughts? Thanks.

[gunzip]

Hi Matt,

Really really useful. Thank you.

I was wondering if you, or anyone else reading this, knew of a guide/workflow on novel transcripts/RNA editing, similar to the level of excellence as this guide.

[MDY]

Hi Azazel,

I agree with you that the "differential expression of whole genes taken from UCSC" approach does ignore some important information provided by RNA-seq. However, I do not agree that this means the approach is without value.

The alternative you suggest is to use the data itself to produce an annotation, against which your analysis can proceed. There will be circumstances where this is indeed a superior approach, but this will not always be the case. Firstly, the ability to annotate a gene depends on the level of coverage, which for RNA-seq depends on the level of transcription. By relying on de-novo annotation routines you will ignore, or at the very least bias against, lowly expressed genes. There are certainly situations where the accurate identification of differential expression in lowly expressed genes is vital to the biological question driving the experiment. The example I am familiar with is known-down experiments of polycomb group proteins, where small expression changes in lowly expressed genes are important, but I am sure there are many others. A related point is that smaller experiments may lack the depth for accurate de novo annotation, even for highly expressed genes.

Furthermore, there is obviously extra information that is still ignored by doing differential expression against de-novo annotated genes, such as differential splicing or allele specific expression.

Ultimately, there are many different things that can be done with RNA-seq data and it seems to me that different analysis techniques offer complementary information rather than opposing, mutually exclusive viewpoints. How important different aspects of RNA-seq analysis are will depend on the biological question you are trying to answer.

The point of this guide is not to be a "one size fits all" guide to analyzing RNA-seq, but to provide a step by step introduction to one of the simpler (and possibly more well understood) analysis methods available. It was my original hope (and it still is) that this guide could form the basis for some kind of "RNA-seq analysis wiki", where people with more expertise with other areas of the analysis could add to it. For example, I think a section describing how to create an annotation using a reference genome and some RNA-seq data would be extremely useful, but I don't have the expertise to write it myself.

[Optimistix]

Hi Azazel, you might want to check out tools from the Salzberg lab, in particular cufflinks, cuffdiff and cuffcompare:

http://cufflinks.cbcb.umd.edu/index.html

I'm new to RNASeq data analysis as well, but those three tools do the kind of thing you seem to have in mind.

[Azazel]

Hi Matt,

thanks for putting up this excellent tutorial.

I have one constructive critisism or discussion point though; as I understand it, when checking for differential expression (DE) you only consider reads "overlapping some annotation object, which is usually something like a collection of genes downloaded from the UCSC."

So you suggest checking DE only for something like RefSeq, and taking the number of reads within each RefSeq (or other object) as the expression level.

I think this discards not only much of the information gained by RNA-seq, but also some of the most important information: the most interesting genes are often among the non-annotated genes. Consider for example two cellular states, a very interesting gene might only be expressed in the very unusual state B, and be very highly expressed; while in state A it's not or so lowly expressed that it didn't make it into the annotation. So with this approach a researcher would miss this gene and others like it entirely because it's not in the annotation, although these might be the very genes which explain the biological question at hand.

If I'd use RNA-seq just to identify DE genes which are already annotated in UCSC, I almost might as well have used a tiling array spanning the annotated genes only. (sure RNA-seq is "digital", but the point I'm trying to make is that with UCSC or similar annotation one would ignore 90%+ of the RNA-seq data elsewhere in the genome!)

So I think a better approach would be first to use the RNA-seq data to produce an ad hoc annotation, including information from all sequenced conditions, then check DE against this annotation.

Now the question is of course, what is a very good way to create an annotation, i.e. how to identify the regions spanned by genes, from RNA-seq?

[colindaven]

Excellent introductory guide, thank you!

[KevinLam]

Excellent. This should be a sticky!

[MDY]

flyyuan - I'm not sure what the answer to your first question about biomart. A detailed description of how bowtie decides on a valid match can be found on the bowtie webpage and in particular the manual. You might want to look at this http://bowtie-bio.sourceforge.net/ma...alignment-mode

In brief, in the default mode bowtie will report a read as matching if it has fewer than -n mismatches from the reference in the seed and the sum of the quality scores at ANY mismatching base within the entire read is less than -e.