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

A clarity on the Illumina TruSeq Small RNA prep kit manual

In the TruSeq® Small RNALibrary Prep Guide, below the Figure 1, there is a sentence: "The RNA 3' adapter is modified to target microRNAs and other small RNAs that have a 3' hydroxyl group resulting from enzymatic cleavage by Dicer or other RNA processing enzymes." It's right, but could be very misleading if you are not clear of the diverse picture of transcriptome (scroll down for more detail). I want to emphasize that the 3' hydroxyl group (and the 5'-phosphate group) is NOT specific to microRNAs or any small RNAs. And it doesn't necessarily result from enzymatic cleavage by Dicer. Sonic fragmentation can also break the full length mRNA (with 5'-cap and 3'-polyA) into truncated RNA pieces with 5'-phosphate and 3' hydroxyl free ends. I just called Illumina to confirm that the 3' and 5' ligation steps don't guarantee the selection of miRNAs (but rather any RNAs with 5'-phosphate and 3' hydroxyl ends, if more accurately). The last step of gel purification is the key to select (or enrich, if more accurately) miRNAs.

OK. Here is what I learned from my colleagues about the different RNA species in the trancriptome:

There are 4 species in the transcriptome, where the later 3 are intermediates of transcription (or half product of degradation).

• me₇G_ppp-------------------------3' (1) 
•                 p------------- 3' (2) 
•                 OH-------------3' (3) 
•        ppp---------------------3' (4) 

Only group(2) will ligate to 5’adaptor. The 3' end can also have different format, at least two:

• 5' ----------------- AAAAAA (1) 
• 5' ---------- OH (2) 

Also note there are two enzymes used to repair the 5' ends: CIP and TAP. CIP (Calf Intestinal alkaline phosphatase) can remove the 5’ phosphate group of DNA strand. The TAP (Tobacco acid pyrophosphatase) is to remove the 5' cap structure (or 5'-5' triphosphate linkage) and leave a mono-phosphate at the 5' end. So, applying first CIP and then TAP will convert the above (2) and (4) to (3), then convert (1) to (2). That's one way to do capture the 5' cap structure, same purpose as CAGE (but CAGE use the type IIs restriction enzyme MmeI and type III restriction enzyme EcoP15I)

Memory issue for large alignment file when doing de novo assembly

Memory becomes a bottleneck when the sequencing file is bigger and bigger nowadays. This is specially an issue for de novo transcriptome assembly using RNA-seq data from species like human. For Trinity, "a typical configuration is a multi-core server with 256 GB to 1 TB of RAM". "Trinity partitions RNA-Seq data into many independent de Bruijn graphs, ideally one graph per expressed gene, and uses parallel computing to reconstruct transcripts from these graphs, including alternatively spliced isoforms." If Trinity contructs one de Bruijn graph per gene, I don't know why it still needs such a large memory. For Cufflinks (-g option), I already can see it consumes 40G memory for 1/10 of chr1 (given a 4G bam file). Cufflinks constructs a DAG in memory for the given alignment. Hopefully it's one DAG per chromosome, not a DAG for the whole genome. But even though, chr1 has 249 billion base pairs, it still requires a lot of memory... So sad!

Of course, if you have a machine with super large memory, this won't be a problem. But this is not the case usually. For me, most of our nodes have maximally 90G memory.

So, how to have it run?

Here is few tips I can think of:

1. split the alignment into different chromosome, e.g. chr1.bam, chr2,bam etc. and then call de novo on each of files.

2. down-sampling the bam file if it has high coverage. Here is a post about how to down-sample your bam file (https://www.biostars.org/p/4332/). Basically, you can use "samtools view -s" or "sambamba view -s", or GATK's randomSampleFromStream.pl for downsampling.

3. You may also want to remove the PCR artifacts, by "samtools rmdup" or "sambamba markdup".

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".

samtools flagstat won't output the number of reads mapped

I have been using samtools flagstat to get statistics summary of my BAM/SAM file, like many tutorial suggested (e.g. Dave Tang's note).

But then I noticed that flagstat reports the number of mapping (or alignments/hits, whatever you like to call), but not the number of reads mapped. This can be very different if your alignment contains multiple mappers (and this is especially true if you use default setting or without setting "-g 1" in Tophat). To get the number of reads mapped (e.g. in order to get mappability), you can use the following command:

samtools view -cF 0x100 accepted_hits.bam

Note: Even with multiple mappers, only one hit has flag of 0x100 in the output SAM file. I've already explained this in one of my previous posts.

R / Bioconductor for High-Throughput Sequence Analysis

I would like to recommend a recent workshop material on R/Bioconductor from Marc Carlson et al.

http://www.bioconductor.org/help/course-materials/2013/SeattleMay2013/

PDF: IntermediateSequenceAnalysis2013.pdf
R script: IntermediateSequenceAnalysis2013.R

Here is the TOC list if you want to read more (highlight for suggested reading list):

Illumina HiSeq2000 adaptor and sequencing

View the large version of the figure

I've referred the following material when making the figure:

http://support.illumina.com/sequencing/sequencing_kits/pe_dna_sample_prep_kit/training.ilmn

http://www.broadinstitute.org/files/shared/illuminavids/clusterGenSlides.pdf

Just one learning note:

If insert is not long enough (i.e. shorter than the read length), R1 will have contamination from Rd2 SP (e.g. AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC) and R2 will have contamination from reverse complementary of Rd1 SP (e.g. AGATCGGAAGAGCG)

So, basically you just need to check if the shared complementary part of Rd1 and Rd2 SP (which is AGATCGGAAGAGC) occurs in the reads. If yes, simply trim it and its following part (if any).

Note: if you don't understand the "shared complementary part", please refer my previous blog on Illumina adaptor. Here is the link: http://onetipperday.blogspot.com/2013/06/illumina-hiseq2000-adaptor.html

Here is one solution of howto remove adaptor contamination:

1. save the complementary part into a fasta file, e.g. adaptor.fa

>adaptor_complementary_part

AGATCGGAAGAGC

2. run fastq-mcf to remove adaptor

fastq-mcf -o filted -x 10 -l 15 -w 4 -u adaptor.fa input.fq.gz

How to select MACS peaks based on p-value, fold_enrichment and FDR?

In the output peaks.xls file from MACS(v1.4), there are 3 columns: -10log10(pvalue), fold_enrichment, and FDR(%). We can use these three columns to sort/filter peaks. Here are what they means:

The '−10*log10(Pvalue)' column lists the transformed P value of each peak, which makes peak sorting easier. For example, a P value of 1e−5 would be transformed to 50. The 'fold_enrichment' column shows the ratio of the ChIP-seq read count to the local value of lambda within each peak. The 'FDR(%)' column contains the empirical FDR percentage for each peak. For example, the fourth peak in the list has an 'FDR(%)' value of '6.45' and '−10*log10(Pvalue)' value of '83.67'; using the same P value cutoff of 4.4e−09 = 1083.67/−10, the ratio of the number of peaks identified by MACS after and before exchanging control and ChIP-seq samples is 6.45:100. The FDR column is only available when the control sample is available.

As the MACS paper (Zhang et al., Genome Biology, 2008) clearly stated, MACS uses Poisson distribution to measure the distribution of ChIP-seq tag. Poisson distribution is characterized as having only one parameter, λ, for both mean and variance. The origin of using Poisson distribution to measure ChIPseq tag can be referred to a Nature paper (see its supplementary note):

The number of sequence reads required to map a chromatin feature can be estimated from a simple model.
Suppose that the genome is divided into N non-overlapping bins of fixed size, that a fraction f of these bins contain a particular chromatin feature and that one performs ChIP-Seq with an antibody that enriches the sequence in these bins by a factor of e. If one collects a total of R sequence reads, the number of reads in a bin should approximately follow a Poisson distribution with mean eM for bins containing the feature and M for the other bins, where M = R/N(ef+(1-f)). 

Xianjun: For those who have difficulty to understand the formula, I would suggest to just understand it like this: if there is a single Poisson distribution, the λ should be R/N; when there are two parts/distribution, we just divide N into Nf and N(1-f) and give different weight to them. For the unenriched part, you want to weight the enriched bins with weight e (because it's more enriched). That is R/(Nfe+N(1-f))=M; For the enriched part, you want to weight unenriched bins with 1/e, then the λ is R/(Nf+N(1-f)/e)=eM. (Thanks to Shikui and Sowmya for helping understanding Poisson distribution and the formula!)

MACS uses a dynamic λ to compensate the local fluctuations and biases even in control sample: 

λ_local = max(λ_BG, [λ_1k,] λ_5k, λ_10k)

where λ_BG is a uniform estimation for the whole genome, λ_1k, λ_5k and λ_10k are λ estimated from the 1 kb, 5 kb or 10 kb window centered at the peak location in the control sample, or the ChIP-Seq sample when a control sample is not available (in which case λ_1k is not used). λ_local captures the influence of local biases, and is robust against occasional low tag counts at small local regions. Xianjun: That's why you will get different p-value and fold_enrichment for the same peak when using different control samples.  MACS uses λ_local to calculate the p-value of each candidate peak (Xianjun: Once you have λ, you have the Poisson distribution; once you have a distribution and the observed value, which is the observed reads count in the peak region for this case, then you can get p-value based on the distribution. See wikipedia for p-value computation) and removes potential false positives due to local biases (that is, peaks significantly under λ_BG, but not under λ_local). Candidate peaks with p-values below a user-defined threshold p-value (default 10^-5) are called, and the ratio between the ChIP-Seq tag count and λ_local is reported as the fold_enrichment. 

The question is: how to select MACS peaks based on the measurements?

To illustrate the relationship of the three measurements, I plot them as below. Note: the figures below are generated by Excel. Looking not so nice, but only for illustration purpose. 

X-axis: -10log10(p-value); Y-axis: FDR (%)

X-axis: fold_enrichment;  Y-axis: FDR (%)

I don't know why FDR has a functional relationship with the p-value (Neither Tao).

But we can choose a smaller FDR cutoff (e.g. 5%), since it's underestimated due to a small control size in our case. That would correspond to a p-value of 10^-50, significant enough.

Here is the fold_enrichment histogram if we choose FDR cutoff as 10%.

Most peaks remained have a >5 fold enrichment, should be fine as well. 

methods of calling differential region of ChIP-seq

Related papers to read:

Model-based Analysis of ChIP-Seq (MACS)

MACS can also be applied to differential binding between two conditions by treating one of the samples as the control. Since peaks from either sample are likely to be biologically meaningful in this case, we cannot use a sample swap to calculate FDR, and the data quality of each sample needs to be evaluated against a real control.

Curr Protoc Bioinformatics. Author manuscript; available in PMC 2012 June 1.MACS empirically calculates FDR based on the number of peaks from control over ChIP that are called at the same p-value cutoff. Therefore if no control data is available, the FDR column does not exist in the output tabular file. Technically, MACS can also be applied to identify differential peaks between two conditions by treating one of the samples as the control. However, calculated FDR value should be ignored, as peaks from either sample are likely to be biologically meaningful in this case.  

http://www.nature.com/nprot/journal/v7/n9/full/nprot.2012.101.htmlthe warning message 'unbalanced reads between treatment and control' means that the FDR of the resulting peaks will be overestimated when the control sample has more reads and will be underestimated when the ChIP-seq sample is sequenced more deeply.

MAnorm: a robust model for quantitative comparison of ChIP-Seq data sets

ChIPdiff and MACS identified four to six times more target regions associated with significantly increased ChIP-Seq signals for K562 cells compared with those found for H1 ES cells, whereas MAnorm yielded a similar number of cell type-biased peaks in each cell line. To compare the enrichment of cell type-specifically expressed genes in the sets of target genes of the differential binding regions discovered by the three methods, we selected the same number of target genes associated with top differential binding regions identified by each method. The target genes of top differential binding regions identified by MAnorm contained similar numbers of H1 ES cell highly expressed genes but a greater number of K562 cell highly expressed genes compared to those identified by ChIPdiff and MACS (Supplementary Table 1 in Additional file 4), suggesting MAnorm performs better in detecting differentially binding regions than the other two methods. Importantly, the fold changes of differential binding given by ChIPdiff and MACS were based on the total number of reads, which may not be appropriate, as discussed above. Additionally, MAnorm showed even better enrichment of cell type-specifically expressed genes in differential binding region targets than the method developed by Taslim et al. [12] when applied to ChIP-Seq data presented in their study (Supplementary Table 2 in Additional file 4).

1	Supplementary Table 1. Enrichment of cell-type differentially expressed genes in genes near differential binding peaks defined by MAnorm, MACS and ChIPdiff. To compare the enrichment scores, we selected the same number of target genes associated with top differential binding regions identified by MAnorm with those identified by other two methods.
2	Supplementary Table 1A. Enrichment of genes more highly expressed in H1 ES cells (as compared to K562) in genes near H1 ES enriched peaks  (as compared to K562) defined by MAnorm, MACS and ChIPdiff
3	H3K27ac H1 ES-enriched target genes	Number of Genes	Overlap with H1 ES up-regulated genes	Enrichment Score
4	MAnorm (M>1)	2680	1243	2.49
5	ChIPdiff (default)	1467	884	3.24
6	MAnorm (top 1467 genes; same number of genes as identified by ChIPdiff with default settings)	1467	941	3.45
7	MACS (P<1e-6)	1589	993	3.36
8	MAnorm (top 1589 genes; same number of genes as identified by MACS with P<1e-6)	1589	987	3.34
9
10
11	Supplementary Table 1B. Enrichment of K562 higher expressed genes (as compared to H1 ES) in genes near K562 enriched peaks (as compared to H1 ES) defined by MAnorm, MACS and ChIPdiff
12	H3K27ac K562-enriched target genes	Number of Genes	Overlap with K562 up-regulated genes	Enrichment Score
13	MAnorm (M<-1)	2694	895	2.78
14	ChIPdiff (default)	6733	1402	1.74
15	ChIPdiff (confidence threshold=0.9999999999)	2291	697	2.55
16	MAnorm (top 2291 genes; same number of genes as identified by ChIPdiff with confidence threshold=0.9999999999)	2291	820	3.00
17	MACS (P<1e-6)	9346	1600	1.43
18	MACS (P<1e-150)	1556	567	3.05
19	MAnorm (top 1556 genes; same number of genes as identified by MACS with P<1e-6))	1556	644	3.47
20

diffReps: Detecting Differential Chromatin Modification Sites from ChIP-seq Data with Biological Replicates
...shows that diffReps is the most sensitive method among all the methods compared, followed by edgeR, DESeq, ChIPDiff and, lastly, CCAT+DESeq. At each cutoff, diffReps (negative binomial test) typically detects a few thousands more differential sites than the secondly ranked method, edgeR.
The tool is in Perl (https://code.google.com/p/diffreps/)

Related tools in R:

Detecting differential binding of transcription factors with ChIP-seq, Bioinformatics, 28, 121-122, doi: 10.1093/bioinformatics/btr605
DBChIP: http://pages.cs.wisc.edu/~kliang/DBChIP/DBChIP.pdf

DiffBind : differential binding analysis of ChIP-Seq peak datahttp://bioconductor.org/packages/2.12/bioc/vignettes/DiffBind/inst/doc/DiffBind.pdf

I'd like to make an assessment for the tools later. 

compatible command line for different input format (*.fastq, *.fq, or *.fq.gz)

For example, if I want to run bowtie with fastq or zipped fastq files, here is the workout:

(([[ $i =~ "fq|fastq$" ]] && cat $i) || ([[ $i =~ "gz$" ]] && zcat $i)) | bowtie genome -v 1 -M 1 --best --strata -q - -S accepted_hits.sam 2> mapping.summary

In above example, bowtie will take input in format of *.fastq, *.fq, or *.gz and output the best hit (or one of the best hits if multiple hits map equally best) for each read with 1 mismatch allowed, and redirect the log to mapping.summary file.

Illumina HiSeq2000 adaptor

It's made just to help myself understanding the adaptor structure of Illumina sequencing. I will add the part of sequencing (on flow cell) later.

Continue to read the following part:

http://onetipperday.blogspot.com/2013/12/illumina-hiseq2000-adaptor-and.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:

Unstable output from fastq-mcf (adaptor removal)

Below is what I wrote to the author regarding the problem I just found. Erik suggested to use "-t 0" to fix the problem (see info at the end)

================

I am really frustrated with a problem I am facing while using the fastq-mcf for adaptor removal -- the same read in different files with different number of reads will get different trimming result. WHAT'S THE PROBLEM?

Here is the test example:

adaptor:

$ cat adaptor.fa 

>adapter1

CAAGCAGAAGACGGCATACGAGATTCAAGTGTGACTGGAGTTCCTTGGCACCCGAGAATTCCA

>RC_adapter1

TGGAATTCTCGGGTGCCAAGGAACTCCAGTCACACTTGAATCTCGTATGCCGTCTTCTGCTTG

Input file:

$ head test.fa

@HWI-ST560:107:D25EJACXX:3:1301:11693:29455 1:N:0:GCCTAA

CAGTGCAATGTTAAAAGGGCATTTGGAATTCTCGGGTGCCAAGGAACTCC

+

CCCDFFFFHHHHGJJJJJJJJJJJJJJJJJJJJJJIHHJJIJJJIJJIIJ

@HWI-ST560:107:D25EJACXX:3:1101:1420:1954 1:N:0:TTAGGC

ATACAGTCCGACGATCTGGAATTCTCGGGTGCCAAGGAACTCCAGTCACT

+

@@@DDDDDHHG?FDGIIBEHAFGIIHFHIHGIIEHGGGHFHIIGGGGGHI

When I include different number of reads, the result are differnt for the 1st read:

$ head -n10000 test.fa > test1.fq; ~/bin/ea-utils.1.1.2-537/fastq-mcf -o test1.out -l 16 -q 15 -w 4 -x 20 -u -S -P 33 adaptor.fa test1.fq; head test1.out

@HWI-ST560:107:D25EJACXX:3:1301:11693:29455 1:N:0:GCCTAA

CAGTGCAATGTTAAAAGGGCATTTGGAATTCTCGGGTGCCAAGGAACTCC

+

CCCDFFFFHHHHGJJJJJJJJJJJJJJJJJJJJJJIHHJJIJJJIJJIIJ

This is not correct, as the yellow highlight part is from RC_adapter1 obviously. 

But if I set a smaller number of sample, it will be correct. See below for "head -n100":

$ head -n100 test.fa > test1.fq; ~/bin/ea-utils.1.1.2-537/fastq-mcf -o test1.out -l 16 -q 15 -w 4 -x 20 -u -S -P 33 adaptor.fa test1.fq; head test1.out

@HWI-ST560:107:D25EJACXX:3:1301:11693:29455 1:N:0:GCCTAA

CAGTGCAATGTTAAAAGGGCATT

+

CCCDFFFFHHHHGJJJJJJJJJJ

And if I set "-C 0", it will get correct result as well (even for a larger number of sample):

$ head -n10000 test.fa > test1.fq; ~/bin/ea-utils.1.1.2-537/fastq-mcf -o test1.out -l 16 -q 15 -w 4 -x 20 -u -S -P 33 -C 0 adaptor.fa test1.fq; head test1.out

@HWI-ST560:107:D25EJACXX:3:1301:11693:29455 1:N:0:GCCTAA

CAGTGCAATGTTAAAAGGGCATT

+

CCCDFFFFHHHHGJJJJJJJJJJ

This is very confusing! Please help me on it.

=========== reply from Erik========
fastq-mcf will not bother to clip if the adaptor is rare (< .25% presence in sampling).   The reason is that the dangers of clipping (removing valid sequence) can outweigh the benefits (removing invalid sequence) as the number of adaptor sequences get fewer and fewer.

To override this, you set -t 0 (always clip everything)

Extract unmapped / unpaired reads using bowtie2

There are options in Bowtie2 for unmapped reads:

--un <path>	Write unpaired reads that fail to align to file at <path>. These reads correspond to the SAM records with the FLAGS 0x4 bit set and neither the 0x40 nor 0x80 bits set.
--al <path>	Write unpaired reads that align at least once to file at <path>. These reads correspond to the SAM records with the FLAGS 0x4, 0x40, and 0x80 bits unset.

--un is path to save unmapped reads (in FASTQ format), and --al is to save all aligned reads (including multiple mappers).  What's worthy to note is, --un only works when --no-unal option is unset (i.e. if --no-unal is set, --un <file> will be empty).

These two options are both for single-reads (SR). For paired-end(PE) reads, there are another two corresponding options:

--un-conc <path>

Write paired-end reads that fail to align concordantly to file(s) at <path>. These reads correspond to the SAM records with the FLAGS 0x4 bit set and either the 0x40 or 0x80 bit set (depending on whether it's mate #1 or #2). .1 and .2 strings are added to the filename to distinguish which file contains mate #1 and mate #2. If a percent symbol, %, is used in <path>, the percent symbol is replaced with 1 or 2 to make the per-mate filenames. Otherwise, .1 or .2 are added before the final dot in <path> to make the per-mate filenames.

--al-conc <path>

Write paired-end reads that align concordantly at least once to file(s) at <path>. These reads correspond to the SAM records with the FLAGS 0x4 bit unset and either the 0x40 or 0x80 bit set (depending on whether it's mate #1 or #2). .1 and .2 strings are added to the filename to distinguish which file contains mate #1 and mate #2. If a percent symbol, %, is used in <path>, the percent symbol is replaced with 1 or 2 to make the per-mate filenames. Otherwise, .1 or .2 are added before the final dot in <path> to make the per-mate filenames.

As above, --un-conc is for improperly paired reads ('discordantly', e.g. both mates are uniquely mapped, but aren't in the expcted relative orientation, or aren't within the expected disatance range, or both, controlled by --fr/--rf/--ff, -I, -X options), --al-conc is for concordantly paired reads (including multi-mappers). Respectively, --un-conc only works when --no-discordant is unset (this is what I guess, not test yet)

awk script to correct XS:A tag of Tophat output for strand-specific paired-end reads

It's been noticed that the current Tophat (v2.0.3) can assign wrong XS:A tag for strand-specific paired-end library, at least for some reads. Here is such an example:

HWI-ST560:74:D14ELACXX:4:1201:8827:168386 pr1 chr1 10024104 50 50M = 10020756 -3398 TGGTTCTTGAAACTGCTGGTTCAGCATCTGTGTACTAACATCAATCCCGG IJJJJJJJIIIGJJJJJJJJJJJJJJJJJJJJIJJJJHHHHHFFFFFCCC AS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:50 YT:Z:UU NH:i:1 XS:A:+
HWI-ST560:74:D14ELACXX:4:1201:8827:168386 pR2 chr1 10020756 50 36M1628N14M = 10024104 3398 GATGATTTGAAATATGAGACTTCTAAGGCATAATATTGTTTGCAGTGCAC CCCFFFFFHHHHHJJJJJJJJJJJJJJJJJJJJIIIIJJJJIJIJHJGEC AS:i:0 XM:i:0 XO:i:0 XG:i:0 MD:Z:50 NM:i:0 XS:A:- NH:i:1

This is a dUTP protocol where "R2/r1" FLAG for the read pair indicate the reads are a transcript on the + strand, which means both reads should be assigned as XS:A:+. However /2 is assigned to XS:A:-.

I've written an awk script to solve the problem:

#!/bin/awk -f

BEGIN{

    if(save_discrepancy_to_file!="") system("[ -e " save_discrepancy_to_file " ] && rm " save_discrepancy_to_file);

}

{

    if($1 ~ /^@/) print;

    else

    {

        for(i=1;i<=NF;i++) if($i!~/^XS/) printf("%s\t",$i); else XS0=$i;

        XS1=XS0;

        if($2~/^0x/ || $2~/^[0-9]+$/){   # FLAG in HEX or Decimal format

            if(libtype=="fr-firststrand") XS1=((and($2, 0x10) && and($2, 0x40)) || (and($2,0x80) && !and($2,0x10)))?"XS:A:+":"XS:A:-";

            if(libtype=="fr-secondstrand") XS1=((and($2, 0x10) && and($2, 0x80)) || (and($2,0x40) && !and($2,0x10)))?"XS:A:+":"XS:A:-";

        }

        else if($2~/^[:alpha:]/){   # FLAG in string

            if(libtype=="fr-firststrand") XS1=($2~/r.*1/ || ($2~/2/ && $2!~/r/))?"XS:A:+":"XS:A:-";

            if(libtype=="fr-secondstrand") XS1=($2~/r.*2/|| ($2~/1/ && $2!~/r/))?"XS:A:+":"XS:A:-";

        }

        print XS1;

        if(save_discrepancy_to_file!="" && XS1!=XS0) print >> save_discrepancy_to_file;

    }

}

Three ways to trim adaptor/primer sequences for paired-end reads

1. Understanding the adaptors (skip this part if you're familiar with the Illumina adaptor)

Before trimming anything from the reads, let's get clear what the reads content is.

Taking Trufseq reads (from Illumina HiSeq 2000) as example, here is the read file (fastq) looks like:

$ cat r1.fq 
@3VFXHS1:278:D13Y4ACXX:1:1101:1472:2209 1:N:0:CGATGT
CTGGTATTGTCTCTTCCCACACTGAACTCTGGGGAATTCGATGTGTGGCACAGCCCGGCTCAGCCTGCCCGCTGGTGGGAGCCCCTGGGAAGCTGCGGCGC
+
@@CFDDFFGH>CAEH:CGHIJJJJEIHJJHIJJJ?DHIDIJHGEGHJG;FHC9@B(5@6A=EH:B@B@2=>>B?BDCBD<B52<<ABD?<?B1@A9>B###
@3VFXHS1:278:D13Y4ACXX:1:1101:1434:2224 1:N:0:CGATGT
GGCAGAGCCAATCTTCGGACGTGGTGATTGTCTCCTCTAAGTACAAACAGCGCTATGAGTGTCGCCTGCCAGCTGGAGCTATTCACTTCCAGCGTGAAAGG
+
BC@FFFFFHHHGFHIIJJIJGFHICFCGIHGFHFGGCHD@F?B?BGGHJJIG6D@EHEHHEHCD259?AACD@AC59?,(5>A,;>:@C(::(029?8>@A
@3VFXHS1:278:D13Y4ACXX:1:1101:1712:2247 1:N:0:CGATGT
GTACACTTGAACACATTTTTCTAACCTTAGAAAATACCTACAAGGCCTGTTGTCTTGACCCATTACTCAATTGTCCCTGGCATATTATCTGATCTTCACGT
+
CCCFFFFFHHHHGJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJIJGIIJGGGDHHHIHIHHEIJIJJIIGGIIIIFEHHHHFEFFFDEEEDCEEEDEDC?A
@3VFXHS1:278:D13Y4ACXX:1:1101:3318:2215 1:N:0:CGATGT
AGATCGGAAGAGCACACGTCTGAACTCCAGTCACCGATGTATCTCGTATGCCGTCTTCTGCTTGAAAAAAAAAAGAAAAAAAAACAAGCGACAAGGACAGA
+
CCCFFFFFHHHHHJJJJJIJJJJJJJJJIJGHFFHIAHIFGGIIJJIJJFIJIHJIHHHGEGFE>CFFEB###############################
@3VFXHS1:278:D13Y4ACXX:2:1101:5344:2243 1:N:0:ACAGTG
GATCGGAAGAGCACACGTCTGAACTCCAGTCACACAGTGATCTCGTATGCCGTCTTCTGCTTGAAAAAAAAACAACAGAAAAAACAAAGCGCGAACAGTGC
+
CCCFFFFFHHHHHJIJJIJIJHJJJJJJJHIGIGIHJFHFFDIIJGIIIJJJHHFIJJJJJJHCEHHFD################################

$ cat r2.fq 
@3VFXHS1:278:D13Y4ACXX:1:1101:1472:2209 2:N:0:CGATGT
CTGGATTTGAAATCTTTAGCGGAGCGGGAACGCCGGCGCGGAAGGGTCTCTACACAGGGCCCGGTCCGCCCTTGCGCTCTCCTTAATGNNNNNNNNNNCGC
+
@CCF?EFFHHHGHHHIFGI@HGGIEHIGIJGI6@@E>B8>??:DBD++399>ACCDDDD@DBDD58@BBDDDDD@@<@BDDDDC>CDC#############
@3VFXHS1:278:D13Y4ACXX:1:1101:1434:2224 2:N:0:CGATGT
CCCGGGGCCTCCCATTAAGGTCGCACTTGGACCCATTGCCATAGGTCTGGCTGTGGTAGCGTTTAAGACGATGCTGCTTGGAGGCCTTGGCTGTTTCATCA
+
BCCFFFFDHFFHGIHIHGJJIJJGGIIJJJIDHGJIJEIIIIIJJIIJIJHHHEEF;>DDA>BBB@CAABBBDDCDDDDDDAD@@?CDDDCCB?ACCDDC#
@3VFXHS1:278:D13Y4ACXX:1:1101:1712:2247 2:N:0:CGATGT
GTTACTCAGCATTTATTCATGCCTGCTGTGTACGGAAAGGGCAGTTACAAAGGAAAGCCTTGATGATTCTGCTTCCAAGAAACGTGAAGATCAGATAATAT
+
CBCFFFFFHHHHHJIJJIIJJJJJIIJIHIHIJJJIJIJIJJIIIIIIJJGGIIJJJIJJJJHIJJJJHIJHGHHHHFFFFFEDECBDDDDDCCDDDCDEE
@3VFXHS1:278:D13Y4ACXX:1:1101:3318:2215 2:N:0:CGATGT
AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATCATTAAAAAAAAAAAAAAAACAGAAAAAGAAAGAAGGAAAGGNGAGT
+
CCCFFFFFHHHHHJJJJFHHHJJJIJJGGIDGHIEIJJJIIJGIJ5@FHIB?DFFEEEDEEEDD?####################################
@3VFXHS1:278:D13Y4ACXX:2:1101:5344:2243 2:N:0:ACAGTG
GGATCGGGAAAGGGGGGGGGGGGGGAAAAGGGGGGATTTCCGGGGGGGCCGGTTCTTTTAAAAAAAAAAAAAAAGAAAACAGAAACAGAAGATGGACAACA
+
CCCFFFFFHHHHHJJJJFHHHJJJIJJGGIDGHIEIJJJIIJGIJ5@FHIB?DFFEEEDEEEDD?####################################

First of first, it's critical to understand what your reads file contain; do they contain adaptor sequences? do they contain primer sequence? I strongly recommend to read the description file here: http://genomics.med.tufts.edu/documents/protocols/TUCF_Understanding_Illumina_TruSeq_Adapters.pdf, from which we could know that the constructed dsDNA (before binding to the flow cell for sequencing) looks like: 

Where 

Trufseq Universal Adaptor:
5´AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCT3´

--> reverse complementary

5´AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATCATT3´

TruSeq Indexed Adapter:
5´GATCGGAAGAGCACACGTCTGAACTCCAGTCAC‐NNNNNN‐ATCTCGTATGCCGTCTTCTGCTTG3´

The 6-nt "NNNNNN" is barcode for multiplexing. We noticed that the 3´ of Universal adaptor is reverse-complementary to the 5´ of Indexed adaptor (Why? This is to form the Y-shape adaptor. See ZZ's dUTP figure in another post). 

Combining together the "Overrepresented sequences" of FASTQC output file (for example), we could infer that the reads are contaminated by adaptor/primers (for example the red-marked part in the fastq sequences).

So, next step is how to remove the contamination. 

2. Collect adaptor sequences

To remove contamination, we first should collect all possible "contamination" sources. For our case, we collect all used barcode sequences and generated the adaptor file:

http://zlab.umassmed.edu/~dongx/tracks/RNAseq/BU/Demultiplex_Stats.htm

copy content and save to txt file Demultiplex_Stats.htm

grep -v Undetermined barcode_stat.html | awk '{print ">"$2; print "GATCGGAAGAGCACACGTCTGAACTCCAGTCAC"$4"ATCTCGTATGCCGTCTTCTGCTTG";}' > adaptor.fa

The Demultiplex_Stats.htm file, which contains barcode information of each sample, is usually included in the output folder of sequencing. Otherwise, you can consult from the data producer. 

You may also want to append the universal adaptor:

$cat >> adaptor.fa

>Trufseq_Universal_Adaptor

AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCT

(Ctrl+D)

A full list of commonly-used adaptors can be retrieved from the following URL (returned by Google search) by command:

curl -s http://www.omicsoft.com/downloads/ngs/contamination_list/v1.txt 2>&1 | sed "/^\s\+$/d;s/\t\+/\t/g;s/ /_/g;s/,//g;s/'//g" | awk '{print ">"$1; print $2;}' > adaptors_list.fa 

This does not include Trufseq barcodes used for smallRNA (see below), which you have to include yourself:

http://epigenome.usc.edu/docs/resources/core_protocols/TruSeq%20index%20sequences-3.pdf

3. Remove adaptor 

There might be many ways to remove adaptors, but specifically for PE reads (e.g. both reads are removed if one of the pairs is disqualified), I'd like to introduce three ways to do this.

a. FastqMcf (http://code.google.com/p/ea-utils/wiki/FastqMcf)

usage: fastq-mcf [options] <adapters.fa> <reads.fq> [mates1.fq ...]

Options:

    -h      This help

    -o FIL  Output file (stats to stdout)

    -s N.N  Log scale for clip pct to threshold (2.2)

    -t N    % occurance threshold before clipping (0.25)

    -m N    Minimum clip length, overrides scaled auto (1)

    -p N    Maximum adapter difference percentage (10)

    -l N    Minimum remaining sequence length (19)

    -L N    Maximum sequence length (none)

    -k N    sKew percentage-less-than causing trim (2)

    -q N    quality threshold causing trimming (10)

    -w N    window-size for quality trimming (1)

    -f      force output, even if not much will be done

    -F FIL  remove sequences that align to FIL

    -0      Set all trimming parameters to zero

    -U|u    Force disable/enable illumina PF filtering

    -P N    phred-scale (auto)

    -x N    'N' (Bad read) percentage causing trimming (20)

    -R      Don't remove N's from the fronts/ends of reads

    -n      Don't clip, just output what would be done

    -C N    Number of reads to use for subsampling (200k)

    -S FIL  Save clipped reads to file

    -d      Output lots of random debugging stuff

For example, 

fastq-mcf -o c1.fq -o c2.fq -l 16 -q 15 -w 4 -x 10 -u -P 33 adaptor.fa r1.fq r2.fq &>r.log

b. Trimmomatic (http://www.usadellab.org/cms/index.php?page=trimmomatic)

Paired End Mode:
java -classpath <path to trimmomatic jar> org.usadellab.trimmomatic.TrimmomaticPE [-threads <threads>] [-phred33 | -phred64] [-trimlog <logFile>] <input 1> <input 2> <paired output 1> <unpaired output 1> <paired output 2> <unpaired output 2> <step 1> ...

Step options:

• ILLUMINACLIP:<fastaWithAdaptersEtc>:<seed mismatches>:<palindrome clip threshold>:<simple clip threshold>
• fastaWithAdaptersEtc: specifies the path to a fasta file containing all the adapters, PCR sequences etc. The naming of the various sequences within this file determines how they are used. See below.
• seedMismatches: specifies the maximum mismatch count which will still allow a full match to be performed
• palindromeClipThreshold: specifies how accurate the match between the two 'adapter ligated' reads must be for PE palindrome read alignment.
• simpleClipThreshold: specifies how accurate the match between any adapter etc. sequence must be against a read.
• SLIDINGWINDOW:<windowSize>:<requiredQuality>
• windowSize: specifies the number of bases to average across
• requiredQuality: specifies the average quality required.
• LEADING:<quality>
• quality: Specifies the minimum quality required to keep a base.
• TRAILING:<quality>
• quality: Specifies the minimum quality required to keep a base.
• CROP:<length>
• length: The number of bases to keep, from the start of the read.
• HEADCROP:<length>
• length: The number of bases to remove from the start of the read.
• MINLENGTH:<length>
• length: Specifies the minimum length of reads to be kept.

Following the above example:

java -classpath $CLASSPATH/trimmomatic-0.22.jar org.usadellab.trimmomatic.TrimmomaticPE -phred33 -trimlog r.log r1.fq r2.fq t1.fq t1.unpaired.fq t2.fq t2.unpaired.fq LEADING:3 TRAILING:3 ILLUMINACLIP:adaptor.fa:2:40:15 SLIDINGWINDOW:4:15 MINLEN:16 

These two programs give same results:

@3VFXHS1:278:D13Y4ACXX:1:1101:1472:2209 1:N:0:CGATGT

CTGGTATTGTCTCTTCCCACACTGAACTCTGGGGAATTCGATGTGTGGCACAGCCCGGCTCAGCCTGCCCGCTGGTGGGAGCCCCTGGGAAGCTGCGG

+

@@CFDDFFGH>CAEH:CGHIJJJJEIHJJHIJJJ?DHIDIJHGEGHJG;FHC9@B(5@6A=EH:B@B@2=>>B?BDCBD<B52<<ABD?<?B1@A9>B

@3VFXHS1:278:D13Y4ACXX:1:1101:1434:2224 1:N:0:CGATGT

GGCAGAGCCAATCTTCGGACGTGGTGATTGTCTCCTCTAAGTACAAACAGCGCTATGAGTGTCGCCTGCCAGCTGGAGCTATTCACTTCCAGCGTGAAAGG

+

BC@FFFFFHHHGFHIIJJIJGFHICFCGIHGFHFGGCHD@F?B?BGGHJJIG6D@EHEHHEHCD259?AACD@AC59?,(5>A,;>:@C(::(029?8>@A

@3VFXHS1:278:D13Y4ACXX:1:1101:1712:2247 1:N:0:CGATGT

GTACACTTGAACACATTTTTCTAACCTTAGAAAATACCTACAAGGCCTGTTGTCTTGACCCATTACTCAATTGTCCCTGGCATATTATCTGATCTTCACGT

+

CCCFFFFFHHHHGJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJIJGIIJGGGDHHHIHIHHEIJIJJIIGGIIIIFEHHHHFEFFFDEEEDCEEEDEDC?A

c. HTSeq 
(http://www-huber.embl.de/users/anders/HTSeq/doc/sequences.html#HTSeq.SequenceWithQualities.trim_left_end_with_quals)

Cite what Simon posted here, HTSeq has facilities useful for this.

You can write a little Python script like this:

import itertools
import HTSeq
in1 = iter( HTSeq.FastqReader( "mydata_1.fastq" ) )
in2 = iter( HTSeq.FastqReader( "mydata_2.fastq" ) )
out1 = open( "trimmed_1.fastq", "w" )
out2 = open( "trimmed_2.fastq", "w" )
for read1, read2 in itertools.izip( in1, in2 ):
   read1.trim_right_end( "ACGGTC" )
   read2.trim_left_end( "TTCGAC" )
   read1.write_to_fastq_file( out1 )
   read2.write_to_fastq_file( out2 )
out1.close()
out2.close()

I've not figured out how to apply a fasta file of adaptor. 

How to tell which library type to use (fr-firststrand or fr-secondstrand)?

First of all, as a bioinformatian, you should ask the data producer (e.g. the one who prepared the RNAseq library) which protocol they used to generate the data.

Tophat manual page has listed the general strand-specific protocol:

Library Type	Examples	Description
fr-unstranded	Standard Illumina	Reads from the left-most end of the fragment (in transcript coordinates) map to the transcript strand, and the right-most end maps to the opposite strand.
fr-firststrand	dUTP, NSR, NNSR	Same as above except we enforce the rule that the right-most end of the fragment (in transcript coordinates) is the first sequenced (or only sequenced for single-end reads). Equivalently, it is assumed that only the strand generated during first strand synthesis is sequenced.
fr-secondstrand	Ligation, Standard SOLiD	Same as above except we enforce the rule that the left-most end of the fragment (in transcript coordinates) is the first sequenced (or only sequenced for single-end reads). Equivalently, it is assumed that only the strand generated during second strand synthesis is sequenced.

In case you don't know the library-type, you can still figure it out by yourself. Tophat FAQ page provided a solution for that (http://tophat.cbcb.umd.edu/faq.html#library_type). But more simply (comparing to running 1M reads first), you can choose few reads and BLAT to genome and infer the library-type from the mapping result.

Generally, reads from the left-most end of RNA fragment (always from 5´ to 3´) are always mapped to transcript-strand, and (for pair-end sequencing) reads from the right-most end are always mapped to the opposite strand. See the arrows direction in the below schema. This is because the sequencer always read from 5´ to 3´.

Summary of library type protocols (for Tophat/Bowtie)

But regarding to which strand the RNA fragment is synthesized from, this involves different strand-specific protocols. Thanks to the illustration figure (see below) from Zhao Zhang, we could see that for example dUTP method is to only sequence the strand from the first strand synthesis (the original RNA strand is  degradated due to the dUTP incorporated), so the /2 read is from the original RNA strand.

Strand-specific library protocols (Credit: Zhao Zhang)

Taking a real example, first getting some reads (in fasta format) from the paired-end sequencing fastq file using command like:

$ zcat ~/nearline/rnaseq/BU/Jul2012/Sample_3576_H_01.R1.fastq.gz | sed 's/@//g;s/ /_/g' | awk '{if(NR%4==1)print ">"$0;if(NR%4==2) print $0;}' | head

$ zcat ~/nearline/rnaseq/BU/Jul2012/Sample_3576_H_01.R2.fastq.gz | sed 's/@//g;s/ /_/g' | awk '{if(NR%4==1)print ">"$0;if(NR%4==2) print $0;}' | head

Blatting them in UCSC Genome Browser

Below is screenshot for top hits of one pair of reads. They mapped to exons of OS9 genes (the left one is /1 and right one is /2, with opposite direction). We see that /1 mapped to transcript direction, /2 mapped to opposite direction, which means it can only be fr-secondstrand or fr-unstrand (cannot be fr-firststrand).

Continuing to look at other reads in the file, we can find examples like these:

where /2 mapped to transcript strand and /1 mapped to the opposite strand. Combining with the observation from above, we can conclude that this is a fr-unstrand library.