Transcriptomics – towards RNASeq – part III

1. Transcriptomics – towardsRNASeq – part III Federico M. Giorgi – federico.giorgi@gmail.com Analisi del Genoma e BioinformaticaCorso di LaureaSpecialistica in Biotecnologiedellepiante e deglianimali

2. Overview of the course Transcription and Transcriptomics Day 1 23/04/2012 Room β3 Transcriptomics methods: Microarrays Exercises on Microarray analysis RNASeq Day 2 02/05/2012 Room 35 Further Applications of RNASeq Day 3 07/05/2012 Room β3 RNASeq exercises

3. Differential Expression in RNASeq The pipeline for performing Differential Expression of all genes between samples can be summarized into five steps: Map reads onto a reference genome or transcriptome Summarize:i.e.merge counts at a gene-, exon-, isoform- level Normalizeread counts (optional) Calculatefold change,p-values, ranking, differential expression Interpretethe results (Systems Biology)

4. Differential Expression in RNASeq The pipeline for performing Differential Expression of all genes between samples can be summarized into five steps: Map reads onto a reference genome or transcriptome Filter reads Tophat (Spliced exon first aligner) Summarize:i.e.merge counts at a gene-, exon-, isoform- level Normalizeread counts (optional) Cufflinks Calculatefold change,p-values, ranking, differential expression Interpretethe results (Systems Biology)

5. Differential Expression in RNASeq In real cases, raw reads are filtered prior to mapping: Generate reads via deep sequencing (e.g. Illumina)

6. Differential Expression in RNASeq In real cases, raw reads are filtered prior to mapping: Generate reads via deep sequencing (e.g. Illumina) Filter reads

7. Differential Expression in RNASeq In real cases, raw reads are filtered prior to mapping: • Remove contaminant reads • - Bowtie aligner vs. NCBI-BLAST sequence database • - rNA (Vezzi et al., 2001) Generate reads via deep sequencing (e.g. Illumina) Filter reads

8. Differential Expression in RNASeq In real cases, raw reads are filtered prior to mapping: • Remove contaminant reads • - Bowtie aligner vs. NCBI-BLAST sequence database • - rNA (Vezzi et al., 2001) • Trim reads for quality • - Trimmomatic • - rNA Generate reads via deep sequencing (e.g. Illumina) Filter reads

9. Differential Expression in RNASeq In real cases, raw reads are filtered prior to mapping: • Remove contaminant reads • - Bowtie aligner vs. NCBI-BLAST sequence database • - rNA (Vezzi et al., 2001) • Trim reads for quality • - Trimmomatic • - rNA • Correct reads by kmer-frequency* • - Soap corrector • - Allpaths-LG corrector (within Trinity) Generate reads via deep sequencing (e.g. Illumina) Filter reads *Seldom used in quantitative analysis Mostly used in Transcriptome reconstruction based on unnormalized libraries

10. Differential Expression in RNASeq In real cases, raw reads are filtered prior to mapping: • Remove contaminant reads • - Bowtie aligner vs. NCBI-BLAST sequence database • - rNA (Vezzi et al., 2001) • Trim reads for quality • - Trimmomatic • - rNA • Correct reads by kmer-frequency • - Soap corrector • - Allpaths-LG corrector (within Trinity) Generate reads via deep sequencing (e.g. Illumina) Filter reads Map reads onto a reference genome or transcriptome Tophat

11. A real example – contaminants in tomato genomic sequences (Solanumlycopersicum) 80% Completely new transcripts, or mistakes? 6% 4% Environmental samples, no species assigned <0.5% Unknown species Possible true contaminants <0.1% Chloroplasts Other «contaminants»: Mitochondria BLAST: assembled genomic contigsvs. NCBI nr/nt

12. Today: RNASeq excercises I see and I forget I read and I remember I do and I understand Confucius 551 B.C. – 479 B.C.

13. Exercise of today Purpose: Differential Expression Analysis Tools: Tophat/Cufflinks software Alternative approaches «Red» way «Green» way Completely relies on the availableannotation of the Transcriptome Use reads to identify novel junctions and/or transcripts

14. Workflow Raw reads format: fastQ Trimming (rNA) Trimmed reads format: fastQ Mapping (TopHat) Mappings format: SAM Use a transcript database to quantify transcript

15. Workflow Raw reads format: fastQ Trimming (rNA) Trimmed reads format: fastQ Mapping (TopHat) Mappings format: SAM Use a transcript database to quantify transcript «Red» way

16. Workflow Raw reads format: fastQ «Green» way Trimming (rNA) Use reads to identify novel junctions and/or transcripts Trimmed reads format: fastQ Mapping (TopHat) Mappings format: SAM Use a transcript database to quantify transcript «Red» way

17. Workflow Raw reads format: fastQ «Green» way Trimming (rNA) Use reads to identify novel junctions and/or transcripts Generate an updated transcript database Trimmed reads format: fastQ Mapping (TopHat) Mappings format: SAM Use a transcript database to quantify transcript «Red» way

18. Experimental design Simple, treated vs. control with replicates: Control A Treated A vs. Control B Treated B Every sample has two fastQ files, because we are gonna work on paired reads This time, we don’t have a real dataset, like with microarrays, because a full RNASeq pipeline study may take days on a cluster of computers

19. Prepare the directory structure Open the terminal cd /home/ngs/data_crunching/RNAseq mkdir cuffcompare_dir mkdir cuffmerge_dir mkdir cuffdiff_red mkdir cuffdiff_green mkdir /home/ngs/data_crunching/sequences/rna/trimmed cd /home/ngs/data_crunching/sequences/rna

20. Look at the reads How long are these reads? cd /home/ngs/data_crunching/sequences/rna

21. Look at the reads How long are these reads? cd /home/ngs/data_crunching/sequences/rna head controlA_R1.fastq

22. Look at the reads How long are these reads? cd /home/ngs/data_crunching/sequences/rna head controlA_R1.fastq How many reads do we have?

23. Look at the reads How long are these reads? cd /home/ngs/data_crunching/sequences/rna head controlA_R1.fastq How many reads do we have? wc -l controlA_R1.fastq grep "@" controlA_R1.fastq| wc -l (Common RNASeq fastQ file: 16,337,998 reads, >2Gbytes)

24. Trim the reads We will use many for loops to save time: for alldir in controlAcontrolBtreatedAtreatedB do echo "$alldir" is a cool directory done

25. Trim the reads for alldir in controlAcontrolBtreatedAtreatedB do rNA --filter-for-assembly --query1 "$alldir"_R1.fastq --query2 "$alldir"_R2.fastq --threads 1 --min-phred-value-CLC 20 --min-mean-phred-quality 20 --min-size 50 --output trimmed/"$alldir" done

26. Index the reference GENOME Go to the genome directory: cd /home/ngs/data_crunching/reference/ Create a GENOME index for the aligner (TopHat): bowtie-build Populus_minusculus.faPopulus_minusculus

27. Trapnell suite of software Trapnell, Nature Protocols, 2012

28. Trapnell suite of software Trapnell, Nature Protocols, 2012

29. Tophat • Aligns RNA-Seq reads to a genome • Able to identify exon-exon splice junctions • Built on the short read aligner Bowtie • Split reads alignment (for reads spanning two exons) • http://tophat.cbcb.umd.edu/manual.html • Or google it!! Read Reference Exon Exon Intron

30. Tophat alignment Change directory, point the annotation to the right file: cd /home/ngs/data_crunching/RNAseq annotation=/home/ngs/data_crunching/genes/Populus_minusculus.gff3 head $annotation reference=/home/ngs/data_crunching/reference/Populus_minusculus head "$reference".fa Align every read pair against the reference: for alldir in controlAcontrolBtreatedAtreatedB do mkdir "$alldir" mkdir "$alldir"/tophat tophat -o "$alldir"/tophat -G $annotation --max-multihits 10 --initial-read-mismatches 1 --segment-mismatches 0 --segment-length 25 $reference ../sequences/rna/trimmed/"$alldir"_1.fastq ../sequences/rna/trimmed/"$alldir"_2.fastq done

31. Some Tophat options -G: provide TopHat with a file of annotation. TopHat will initially map on the transcriptome and then map the remaining reads to the genome (exon-first spliced approach) -g/--max-multihits: ignore reads with more than this alignments. --initial-read-mismatches: max number of mismatches for read alignment --segment-mismatches: max number of mismatches for split segments --segment-length: length of each segment in which reads are splitted for alignment to the genome

32. First approach: red way «Red» way Completely relies on the availableannotation of the Transcriptome

33. Run Cuffdiff Cuffdiff is one of the Cufflinks sub-packages, for simple Differential Expression Analysis cd /home/ngs/data_crunching/RNAseq annotation=/home/ngs/data_crunching/genes/Populus_minusculus.gff3 cuffdiff-o cuffdiff_red/ -L Control,Treated$annotation controlA/tophat/accepted_hits.bam,controlB/tophat/accepted_hits.bamtreatedA/tophat/accepted_hits.bam,treatedB/tophat/accepted_hits.bam

34. Analyze results (1) cd cuffdiff_red head gene_exp.diff • There is more than one file called *.diff • Let’s open the one named “isoform_exp.diff” with LibreOffice Calc • A status “NOTEST” means that for the gene you don’t have enough data. We only can work on genes with “OK”

35. Analyze results (2) • How many significantly changed genes do we have? • Sort by q-value (corrected p-value) • The “significant” column has the value “yes” if the p-valueof the differential expression is <0.05 after correction for multiple testing • This is the set on which you should focus your interest for downstream analysis such as MapMan

36. Analyze results (3) Cuffdiff/Cufflinks doesn’t provide raw read counts, but calculates RPKMs RPKM (Read per kilobase of exon model per million mapped reads): accounts for both library size and gene lengtheffect To be even more precise, these are called FPKMs, to account for paired reads (Fragments per kilobase per million mapped pairs of reads)

37. Second approach: green way «Green» way Use reads to identify novel junctions and/or transcripts

38. The green way • The first step of the green route is still mapping of reads using TopHat (which in turn relies on Bowtie). We can use the same parameters • Now, we will exploit Cufflinks to predict new transcripts • Cufflinks takes a *.bam format file as input. Authors’ suggest to use the output of TopHat as Cufflink input

39. Run Cufflinks (1) First, we need to index the ALIGNMENTS Cufflinks requires indexed alignments... cd /home/ngs/data_crunching/RNAseq for alldir in controlAcontrolBtreatedAtreatedB do samtoolsindex "$alldir"/tophat/accepted_hits.bam done

40. Run Cufflinks (2) And then, we run Cufflinks itself: annotation=/home/ngs/data_crunching/genes/Populus_minusculus.gff3 for alldir in controlAcontrolBtreatedAtreatedB do cufflinks -g $annotation -o "$alldir"/cufflinks "$alldir"/tophat/accepted_hits.bam done You can change this –g with –G. -g: use the annotation to guide assembly (green way) -G: use the annotation and do not perform assembly (if you use the gff you supplied to tophat, then using –G will be the same as following the red way)

41. Merge assemblies (1) Now we are comparing not only different mappings, but also different assemblies, each produced by Cufflinks on a different set of reads In order to compare counts from different cufflinks assemblies/mappings, we need a further merging step Create a manifest file, that points cuffmerge to the assemblies produced by the two separate Cufflinks runs: cd /home/ngs/data_crunching/RNAseq nano manifest.txt controlA/cufflinks/transcripts.gtf controlB/cufflinks/transcripts.gtf treatedA/cufflinks/transcripts.gtf treatedB/cufflinks/transcripts.gtf Fill the file with four lines Ctrl + o (save) ENTER Ctrl + x (exit)

42. Merge assemblies (2) Simply run cuffmerge: cuffmerge -o cuffmerge_dir -g $annotation manifest.txt

43. Run Cuffdiff Cuffdiff will now work as in the red way: it will just compare comparable gene counts: cd /home/ngs/data_crunching/RNAseq newannotation=cuffmerge_dir/merged.gtf cuffdiff-o cuffdiff_green/ -L Control,Treated $newannotationcontrolA/tophat/accepted_hits.bam,controlB/tophat/accepted_hits.bamtreatedA/tophat/accepted_hits.bam,treatedB/tophat/accepted_hits.bam This is an annotation generated by cuffmerge: it’s an extension of the default one

44. Analyze results cd cuffdiff_green head gene_exp.diff How many significantly changed genes do we have now? Is this identical to the red way? You can notice that NEW genomic areas previously not known to be transcribed were found by Cufflinks in the green way

45. Conclusions • You now have a list of genes which are differentially expressed between two conditions • You can assess the effects of this treatment by checking who these genes are • In this case, we used a (dummy) dataset from Populus minuscula, a not very studied plant, therefore before using Mapman we must assign a functionto the (known and newly found by Cufflinks) gene sequences BLAST Mercator http://mapman.gabipd.org/web/guest/app/mercator

46. Final slide Fabio Marroni And thanks to who provided the data!

47. CummeRbund CummeRbund is an R package to load and visualize Cuffdiff outputs Install CummeRbund: Here, select «BioC Software» R setRepositories() install.packages("cummeRbund") library("cummeRbund") Load the cuffdiff "green way" output: cuff<-readCufflinks(dir="/home/ngs/data_crunching/RNAseq/cuffdiff_green/") cuff

48. CummeRbund CummeRbund is an R package to load and visualize Cuffdiff outputs Some plots describing the dataset: disp<-dispersionPlot(genes(cuff)) disp

49. CummeRbund CummeRbund is an R package to load and visualize Cuffdiff outputs Some plots describing the dataset: brep<-csBoxplot(genes(cuff),replicates=T) brep

50. CummeRbund CummeRbund is an R package to load and visualize Cuffdiff outputs Some plots describing the dataset: dens<-csDensity (genes(cuff),replicates=T) dens