Note from GEUVADIS paper

Note from GEUVADIS papers:

Lappalainen et al. Nature 2013 : Transcriptome and genome sequencing uncovers functional variation in humans, http://dx.doi.org/10.1038/nature12531

 

‘t Hoen et al. Nature Biotechnology 2013: Reproducibility of high-throughput mRNA and small RNA sequencing across laboratories, http://dx.doi.org/10.1038/nbt.2702

1. how to detect sample outlier?

a. before alignment: distance of k-mer profile
b. after alignment: Spearman rank correlation between samples --> D-statistics (i.e.  the median correlation of one sample against all the other samples)
c. gender mismatch: XIST vs. chrY
d. ASE bias rate among heterozygous sites

2. eQTL

a. exon/gene quantification
b. filter out lowly expressed ones (e.g. 0 in >50% samples)
c. for each group, normalize with PEER, adding mean
    c1. use subset (??e.g. chr20, or chr20- 22??) using K=????0,1,,3,,57,10,13,15,20 for each dataset
    c2. run eQTL and number of genes (eGenes) for each K.
    c3. get the optimal K = K(with most number of eQTL genes)
    c4. run PEER on 20,000 exons to get covairantes for the final normalization
    c5. final PEER normalization using all dataset, residual + mean as final quantification
d. transform the final quantification to standard normal distribution (by ?)
e. eQTL using Matrix-eQTL: linear regression of quantification ~ genotypes + genotype_covariates

3. Differential expression analysis

a. TMM normalization (from edgeR)
b. filter: genes with more than 5 counts per million in at least 1 sample were analyzed in pairwise population comparisons
c. tweeDEseq (good for large samples), significance: FDR < 0.05 and log2 fold change greater than 2

Cufflinks mask option (-M/--mask-file) works when ...

Obviously I am not the only one who had questions on the "-M/--mask-file" mask GTF option in Cufflinks:

http://seqanswers.com/forums/showthread.php?t=8190
https://www.biostars.org/p/110289/
http://seqanswers.com/forums/showthread.php?t=29975

And too bad that no one from the Texedo group ever threw a piece of clue!

Here are few tips I found necessary to share in order to have it work:

1. The mask GTF file should have all 9 fields in required format. For example, the strand column should be '+', '-', or '.', not anything else. GTF/GFF file can be extracted from GENCODE (http://www.gencodegenes.org/) or downloaded from UCSC Table browser. It can be also converted from a bed file using Kent's bedToGenePred --> genePredToGtf. But be aware that that the bed file should have at least 6 columns (i.e. including strand column), otherwise the converted GTF file will have a "^@" in the strand column, which results in an invalid GTF.

For example, if you want to exclude all reads mapped to human mitochondrial genome,  you can use
echo "chrM 0 16571 mt 0 ." | bedToGenePred stdin stdout | genePredToGtf file stdin chrM.gtf

2. "-M" option also works for de novo assembly (cufflinks -g).

3. Using "-M" option should theoretically increase the FPKM value (comparing to no mask). So, if you observed opposite tread, there must be something wrong.

4. If you expect a lot of reads from the mask regions (e.g. chrM, rRNAs), you can substract the masked reads from your bam file before feeding to cufflinks, for example using "samtools view -L retained_region.bed".

Natural barcode and rRNA genes

Ribosomal RNA (rRNA) is part of ribosome, where proteins are synthesized from mRNAs. Ribosome is composed by two major types of rRNAs (large subunit -LSU and short subunit - SSU) and >50 ribosomal proteins.

Credit: http://en.wikipedia.org/wiki/File:Ribosome_mRNA_translation_en.svg

The ribosome in eukaryotes is called 80S ribosome. S is the unit of sedimentation coefficients defined by Svedberg. Its large subunit is called 60S and small subunit is called 40S. In most eukaryotes the small ribosomal subunit contains 18S rRNA, while the large subunit contains three rRNA species: the 5S, 5.8S and 28S.

In eukaryotes, the rRNAs are generally encoded by many copied of rRNA genes in the genome. Mammalian cells have 2 mitochondrial (12S and 16S) rRNA molecules and 4 types of cytoplasmic rRNA (the 28S, 5.8S, 18S, and 5S subunits). The 28S, 5.8S, and 18S rRNAs are encoded by a single transcription unit (45S) separated by 2 internally transcribed spacers (ITS) and padded by external transcribed sequences (ETS) at both ends. The 45S rDNA organized into 5 clusters (each has 30-40 repeats) on chromosomes 13, 14, 15, 21, and 22. The region between the clusters is called intergenic spacer (IGS). These are transcribed by RNA polymerase I. (source from wikipedia). See below for details.

Credit: Xianjun Dong

One example of 45S rRNA in Genbank: http://www.ncbi.nlm.nih.gov/gene/100861532

What interested me is that the ITS sequence can be used as barcode to identify different species because of its high degree of variation and ease of amplification.

5S occurs in tandem arrays (~200-300 true 5S genes and many dispersed pseudogenes), the largest one on the chromosome 1q41-42. 5S rRNA is transcribed by RNA polymerase III. Here is one 5S rRNA record in Genbank: http://www.ncbi.nlm.nih.gov/gene/100873410

All animal mitochondria genomes contain 37 genes, 2 rRNA, 22 tRNA and 13 mRNA.The 2 rRNA are 12S rRNA and 16S rRNA, which are encoded by gene MTRNR1 and MT-RNR2, respectively. Research have shown that MT-RNR2 is associated with Alzheimer's disease and MTRNR1 is associated with hearing loss.

Due to the high number of copies of rRNA genes, we usually see big porportion of rRNA from the total RNA. That will bias the RNAseq result. So we need to remove the rRNAs from the library before sequencing (e.g. ribominus) and after sequencing (e.g. cufflinks -M option). 

I have an early post on howto get rRNA gene annotation: http://onetipperday.blogspot.com/2012/08/how-to-get-trnarrnamitochondrial-gene.html

RNA sample quality assessment

RIN: http://www.biomedcentral.com/1471-2199/7/3

The RQN is equivalent to the broadly accepted RIN and is based on a proprietary algorithm that uses three areas of the electropherogram: (A) the area before the 18S peak, (B) the total area of the 18S and 28S peaks, and (C) the ratio of the 28S and 18S peaks. (Reference: http://www.genengnews.com/gen-articles/simplifying-ht-rna-quality-quantity-analysis/4688/)

RIN is highly correlated with RQN:

faster way to get number of reads species

Let's call the unique reads sequences as 'species' in this topic.  So, an easy way to get the number of reads species for a RNAseq (for example) sam file can be

cut -f10 alignment.sam | sort -u | wc -l

But it's not the faster one if the sam file is large. We can improve it by not sorting it (which is unnecessary for this task). Here is it:

cut -f10 alignment.sam | awk '{r[$1]++;}END{for(i in r)j++; print "number of species:", j;}'
or

awk '{r[$10]++;}END{for(i in r)j++; print "number of species:", j;}' alignment.sam

The point is, you don't really need to "sort" the key, but instead to "group" the lines by the key. For example, when using groupBy in bedtools, it requires the input is pre-sorted by the key which is to be grouped. To do that in a quick way (esp. important for large file), we can use array of array for awk or hash of array in perl. Here is the hash of array solution in perl for "sort -k10" of a large file:

perl -e 'while (<>) {$l=$_; @a=split("\t", $l); push(@{$HoA{$a[17]}}, $l);}{foreach $i (sort keys %HoA) {print join("", @{$HoA{$i}});}}'

btw, I found this awk array of array does not work:

awk '{r[$10][length(r[$10])+1]=$0;}END{for(i in r) for (j in r[i]) print r[i][j];}' alignment.sam

how RNAseq quantification tools deal with multi-mapping reads?

Just to give an overview of how the commonly used RNAseq tools deal with multiple-mapping reads.

1. Cufflinks (http://cufflinks.cbcb.umd.edu/howitworks.html#hmul)
By default, Cufflinks will uniformly divide each multi-mapped read to all of the positions it maps to. In other words, a read mapping to 10 positions will count as 10% of a read at each position. If multi-mapped read correction is enabled (-u/--multi-read-correct), Cufflinks will improve its estimation in a 'rescue' manner. That is, Cufflinks will first calculate initial abundance estimates for all transcripts using the uniform dividing scheme. Cufflinks will then re-estimate the abundances dividing each multi-mapped read probabilistically based on the initial abundance estimation of the genes it maps to, the inferred fragment length, and fragment bias (if bias correction is enabled).

2. HTSeq-count (http://www-huber.embl.de/users/anders/HTSeq/doc/count.html)

HTSeq-count can be quite conservative in the sense that it discards reads that cannot be unambiguously be assigned to a gene. This can cause values to be lower than reported by other tools. (This is done on purpose: For differential expression analysis, you need to discount ambiguous reads to avoid that differential signal on one gene shows up in another one that overlaps.) (http://seqanswers.com/forums/showthread.php?t=9129)  Note that, HTSeq-count also reports alignment_not_unique: reads with more than one reported alignment. These reads are recognized from the NH optional SAM field tag. (If the aligner does not set this field, multiply aligned reads will be counted multiple times.)

3. RSEM (http://www.biomedcentral.com/1471-2105/12/323)
It's like the above Cufflinks rescue model, but using EM method to estimate. It also provides a way to generate normalized or weighted bigwig.

4. bedtools multicov

bedtools gives you raw hits which includes every 1 hit for every intersection of every alignment with any features in the GTF no matter how many times it aligned or how many features it hit. So, you have to deal with the apportion yourself, either uniformly or probabilistically.

Don't trust Cufflinks FPKM for short genes

Here is what Cole (author of Cufflinks) commented to the observation of "very high RPKM values from Cufflink":

This issue has been discussed elsewhere on this board. As Nicholas points out, RNA-Seq really isn't reliable for very short transcripts. The reason is that all the fragments that map to these transcripts come from the "tail" of the distribution of library fragment lengths. That is, fragments that map to microRNAs are much, much shorter than most fragments in the library - by design in the RNA-Seq protocol, which size selects away very short inserts. Thus, Cufflinks infers that even though relatively few fragments actually mapped to the microRNAs, there were probably TONS of individual microRNA molecules in the transcriptome before all of the various size selection parts of the protocol kicked in. Cufflinks accordingly increases the FPKM of these short transcripts to compensate for the bias against short fragments in the library.

This compensation was designed to improve accuracy for transcripts that are in the 500bp-1kb range - for longer transcripts, the "edge effects" due to library fragment size aren't much of an issue. However, I wouldn't trust FPKM values for transcripts shorter than your average fragment length. There's really just not enough data in most standard RNA-Seq libraries to say much about small RNA abundance.
I should also point out that other methods use this same bias correction technique (RSEM for example). As far as I'm aware, the "count-based" methods don't, but that doesn't mean they shouldn't. Most of those methods are strictly for differential analysis, where any edge effects are assumed to be affecting each condition the same way. That may or may not be the case in your data.
In any case, the quick answer to this problem is to simply remove or ignore transcripts shorter than around 300bp from your GTF. In a future version, we will be flagging these transcripts as too short for reliable quantification where appropriate.

Here is the full discussion: http://seqanswers.com/forums/archive/index.php/t-17404.html

adaptor remover/clipper with mismatch allowed

$ head -n12 RNAseq/mouse_adult_wt_smallRNAseq_76nt_strand_ox.fastq 
@HWUSI-EAS1533_0026_FC:6:1:1204:14021#0/1
TTCCTATCAATAAATCTCATGTGANNGCAGCTCGTANNCCNTCTNCNNNTNNNNANNNNANANCNNNNNANNNNNN
+HWUSI-EAS1533_0026_FC:6:1:1204:14021#0/1
gggaegggggfdggfcaedbBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB
@HWUSI-EAS1533_0026_FC:6:1:1204:19652#0/1
TGGTATAAACTCTTCTAAAAGACTNNAGCTCGTATGNNGCNTTCNGNNNGNNNNANNNNANANANNNNNANNNNNN
+HWUSI-EAS1533_0026_FC:6:1:1204:19652#0/1
c^^[^a`ddYaddadf_YafSOTZDDYZK\M`W_`BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB
@HWUSI-EAS1533_0026_FC:6:1:1205:19274#0/1
TCAGAAATGTCCTGTGAAACTCGCNNAATTTCGTATGCCGCCTTCTGNNTGNAAANNNAAAACAANNNNANNNNNN
+HWUSI-EAS1533_0026_FC:6:1:1205:19274#0/1
fV^V]b]`ZSTXdRbc]dac`BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB

The basic question is to trim/clip the part after the adaptor sequence (e.g. TCGTATGCCGTCTTCTGCTTG in above example) in the FASTQ reads. It's easy if there is no mismatch. For case with mismatch, here is few solutions I got from google:

1. FAR (http://sourceforge.net/apps/mediawiki/theflexibleadap/index.php?title=Main_Page), which uses Needleman-Wunsch alignment (global alignment) to align the adaptor sequence(s) to the read sequences and find the first best alignment.
   
   far -s stdin.fastq -t stdout.fq -f fastq-sanger -as TCGTATGCCGTCTTCTGCTTG --cut-off 5 --min-overlap 15 --phred-pre-trim 20 --min-readlength 20 --trim-end right --adaptive-overlap yes --max-uncalled 30 --nr-threads 8 --log-level ALL
2.  Scythe only checks for 3'-end contaminants, up to the adapter's length into the 3'-end. For reads with contamination in any position, the program TagDust (http://genome.gsc.riken.jp/osc/english/dataresource/) is recommended. Scythe has the advantages of allowing fuzzier matching and being base quality-aware, while TagDust has the advantages of very fast matching (but allowing few mismatches, and not considering quality) and FDR. TagDust also removes contaminated reads entirely, while Scythe trims off contaminants.

A possible pipeline would run FASTQ reads through Scythe, then TagDust, then a quality-based trimmer, and finally through a read quality statistics program such as qrqc (http://bioconductor.org/packages/devel/bioc/html/qrqc.html) or FASTqc (http://www.bioinformatics.bbsrc.ac.uk/projects/fastqc/).
3. trimLRPatterns() function in BioConductor. Here is an example.
4. Another very useful toolkit is fastx_toolkit, which includes a set of tools for handling fasta and fastq format. fastx_clipper  is for clipping adaptor, but from what I checked it does not allow mismatch (e.g. N is counted as 'match' in fastx_clipper, but in FAR it's mismatch).  Running the following code on above reads only clips the last read. The first one should also be clipped. Don't know why. (Just get reply from Gordon that fastx_clipper does allow mismatch and also count 'N' as mismatch. Just the first two sequences have high error rate (=number_of_mismatch/adaptor_alignment_length, e.g. 33% in this case), so it's not counted as an adaptor. 

$ head -n12 RNAseq/mouse_adult_wt_smallRNAseq_76nt_strand_ox.fastq | fastx_clipper -a TCGTATGCCGTCTTCTGCTTG -c -n -i - 

@HWUSI-EAS1533_0026_FC:6:1:1205:19274#0/1

TCAGAAATGTCCTGTGAAACTCGCNNAATT

+HWUSI-EAS1533_0026_FC:6:1:1205:19274#0/1

fV^V]b]`ZSTXdRbc]dac`BBBBBBBBB