Tuxedo tools I (Tophat and Cufflinks) analysis of S. pombe data from Trinity protocols
Bhagirathi Dash
Data, references and resources
Reference: Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat Protoc. 7(3): 562-578.
Machine AWS m4.xlage
Downloads
~/tophat_spombe$ wget http://sourceforge.net/projects/trinityrnaseq/files/misc/TrinityNatureProtocolTutorial.tgz
………………………………
TrinityNatureProtocolTutorial.tgz 100%[==============================================================================================>] 452.06M 19.7MB/s in 24s
2019-04-01 01:52:20 (18.7 MB/s) - 2018TrinityNatureProtocolTutorial.tgz saved [474015394/474015394]
~/tophat_spombe$ wget ftp://ftp.pombase.org/pombe/genome_sequence_and_features/gff3/Schizosaccharomyces_pombe_all_chromosomes.gff3.gz
…………………………..
Schizosaccharomyces_pombe_all_chromosomes 100%[====================================================================================>] 458.59K 1.11MB/s in 0.4s
2019-04-01 01:29:09 (1.11 MB/s) - ‘Schizosaccharomyces_pombe_all_chromosomes.gff3.gz’ saved [469600]
~/tophat_spombe$ wget ftp://ftp.pombase.org/pombe/genome_sequence_and_features/genome_sequence/Schizosaccharomyces_pombe_all_chromosomes.fa.gz
……………………
Schizosaccharomyces_pombe_all_chromosomes.fa. 100%[==============================================================================================>] 3.60M 4.98MB/s in 0.7s
2019-04-01 01:29:55 (4.98 MB/s) - ‘Schizosaccharomyces_pombe_all_chromosomes.fa.gz’ saved [3773743]
~/tophat_spombe$ wget ftp://ftp.ensemblgenomes.org/pub/current/fungi/gtf/schizosaccharomyces_pombe/Schizosaccharomyces_pombe.ASM294v2.42.gtf.gz
…………………
Schizosaccharomyces_pombe.ASM294v2.42.gtf.gz 100%[==============================================================================================>] 785.85K 1.28MB/s in 0.6s
2019-04-01 01:38:59 (1.28 MB/s) - ‘Schizosaccharomyces_pombe.ASM294v2.42.gtf.gz’ saved [804707]
Directory contents
~/tophat_spombe$ ls -l
-rw-rw-r– 1 ubuntu ubuntu 14544225 Apr 1 01:38 Schizosaccharomyces_pombe.ASM294v2.42.gtf
-rw-rw-r– 1 ubuntu ubuntu 12842130 Apr 1 01:29 Schizosaccharomyces_pombe_all_chromosomes.fa
-rw-rw-r– 1 ubuntu ubuntu 3445429 Apr 1 01:29 Schizosaccharomyces_pombe_all_chromosomes.gff3
-rw-rw-r– 1 ubuntu ubuntu 474015394 Feb 9 2013 TrinityNatureProtocolTutorial.tgz
~/tophat_spombe$ tar xvf TrinityNatureProtocolTutorial.tgz
TrinityNatureProtocolTutorial/ TrinityNatureProtocolTutorial/S_pombe_refTrans.fasta TrinityNatureProtocolTutorial/1M_READS_sample/ TrinityNatureProtocolTutorial/1M_READS_sample/Sp.hs.1M.left.fq TrinityNatureProtocolTutorial/1M_READS_sample/Sp.log.1M.right.fq TrinityNatureProtocolTutorial/1M_READS_sample/Sp.ds.1M.right.fq TrinityNatureProtocolTutorial/1M_READS_sample/Sp.ds.1M.left.fq TrinityNatureProtocolTutorial/1M_READS_sample/Sp.log.1M.left.fq TrinityNatureProtocolTutorial/1M_READS_sample/Sp.plat.1M.left.fq TrinityNatureProtocolTutorial/1M_READS_sample/Sp.plat.1M.right.fq TrinityNatureProtocolTutorial/1M_READS_sample/Sp.hs.1M.right.fq TrinityNatureProtocolTutorial/samples_n_reads_described.txt
~/tophat_spombe$ ls -l TrinityNatureProtocolTutorial
drwxrwxr-x 2 ubuntu ubuntu 4096 Feb 6 2013 1M_READS_sample
-rw-rw-r– 1 ubuntu ubuntu 7830629 Feb 6 2013 S_pombe_refTrans.fasta
-rw-rw-r– 1 ubuntu ubuntu 422 Feb 6 2013 samples_n_reads_described.txt
~/tophat_spombe$ ls -l TrinityNatureProtocolTutorial/1M_READS_sample/
-rw-rw-r– 1 ubuntu ubuntu 175846179 Feb 6 2013 Sp.ds.1M.left.fq
-rw-rw-r– 1 ubuntu ubuntu 175846179 Feb 6 2013 Sp.ds.1M.right.fq
-rw-rw-r– 1 ubuntu ubuntu 175736042 Feb 6 2013 Sp.hs.1M.left.fq
-rw-rw-r– 1 ubuntu ubuntu 175736042 Feb 6 2013 Sp.hs.1M.right.fq
-rw-rw-r– 1 ubuntu ubuntu 175741215 Feb 6 2013 Sp.log.1M.left.fq
-rw-rw-r– 1 ubuntu ubuntu 175741215 Feb 6 2013 Sp.log.1M.right.fq
-rw-rw-r– 1 ubuntu ubuntu 175899533 Feb 6 2013 Sp.plat.1M.left.fq
-rw-rw-r– 1 ubuntu ubuntu 175899533 Feb 6 2013 Sp.plat.1M.right.fq
Build a S. pombe genome index
~/tophat_spombe$ bowtie2-build genome_data/Schizosaccharomyces_pombe_all_chromosomes.fa genome
Settings: **Output files: "genome.*.bt2"** …………………..
Input files DNA, FASTA: genome_data/Schizosaccharomyces_pombe_all_chromosomes.fa
…………………………..
bmax according to bmaxDivN setting: 3157744 Using parameters –bmax 2368308 –dcv 1024 Doing ahead-of-time memory usage test Passed! Constructing with these parameters: –bmax 2368308 –dcv 1024
Constructing suffix-array element generator ………………….
Building DifferenceCoverSample ……………… Invoking Larsson-Sadakane on ranks time: 00:00:00
………………..
Calculating bucket sizes Binary sorting into buckets 10% 20% 30% 40% 50% 60% 70% 80% 90% 100% Binary sorting into buckets time: 00:00:00 Splitting and merging Splitting and merging time: 00:00:00 Split 1, merged 7; iterating…
…………………….
Avg bucket size: 1.57887e+06 (target: 2368307)
Converting suffix-array elements to index image
Allocating ftab, absorbFtab
Entering Ebwt loop
Getting block 1 of 8
………………..
……………… Constructing suffix-array element generator ……………
Exited Ebwt loop
fchr[A]: 0
fchr[C]: 4036348
fchr[G]: 6311338
fchr[T]: 8588523
fchr[$]: 12630977
Exiting Ebwt::buildToDisk()
Returning from initFromVector
Wrote 8405204 bytes to primary EBWT file: genome.rev.1.bt2
Wrote 3157752 bytes to secondary EBWT file: genome.rev.2.bt2
Re-opening _in1 and _in2 as input streams
Returning from Ebwt constructor
Headers:
len: 12630977
bwtLen: 12630978
sz: 3157745
bwtSz: 3157745
lineRate: 6
offRate: 4
offMask: 0xfffffff0
ftabChars: 10
eftabLen: 20
eftabSz: 80
ftabLen: 1048577
ftabSz: 4194308
offsLen: 789437
offsSz: 3157748
lineSz: 64
sideSz: 64
sideBwtSz: 48
sideBwtLen: 192
numSides: 65787
numLines: 65787
ebwtTotLen: 4210368
ebwtTotSz: 4210368
color: 0
reverse: 1Total time for backward call to driver() for mirror index: 00:00:04
Additional directories and contents following the run
~/tophat_spombe$ ls -l
-rw-rw-r– 1 ubuntu ubuntu 8405204 Apr 1 01:43 genome.1.bt2
-rw-rw-r– 1 ubuntu ubuntu 3157752 Apr 1 01:43 genome.2.bt2
-rw-rw-r– 1 ubuntu ubuntu 116 Apr 1 01:43 genome.3.bt2
-rw-rw-r– 1 ubuntu ubuntu 3157745 Apr 1 01:43 genome.4.bt2
-rw-rw-r– 1 ubuntu ubuntu 8405204 Apr 1 01:43 genome.rev.1.bt2
-rw-rw-r– 1 ubuntu ubuntu 3157752 Apr 1 01:43 genome.rev.2.bt2
drwxrwxr-x 2 ubuntu ubuntu 4096 Apr 1 01:42 genome_data
Tophat: Map the reads for each sample to the reference genome
~/tophat_spombe$ tophat -I 1000 -i 20 –library-type fr-firststrand -o tophat.Sp_ds.dir genome TrinityNatureProtocolTutorial/1M_READS_sample/Sp.ds.1M.left.fq TrinityNatureProtocolTutorial/1M_READS_sample/Sp.ds.1M.right.fq
[2019-04-01 02:08:20] Beginning TopHat run (v2.1.1)
[2019-04-01 02:08:20] Checking for Bowtie Bowtie version: 2.2.5.0
[2019-04-01 02:08:20] Checking for Bowtie index files (genome)..
[2019-04-01 02:08:20] Checking for reference FASTA file
Warning: Could not find FASTA file genome.fa[2019-04-01 02:08:20] Reconstituting reference FASTA file from Bowtie index Executing: /home/ubuntu/miniconda3/bin/bowtie2-inspect genome > tophat.Sp_ds.dir/tmp/genome.fa
[2019-04-01 02:08:20] Generating SAM header for genome
[2019-04-01 02:08:20] Preparing reads left reads: min. length=68, max. length=68, 997953 kept reads (2047 discarded)
right reads: min. length=68, max. length=68, 998472 kept reads (1528 discarded)[2019-04-01 02:08:40] Mapping left_kept_reads to genome genome with Bowtie2
[2019-04-01 02:09:24] Mapping left_kept_reads_seg1 to genome genome with Bowtie2 (1/2)
[2019-04-01 02:09:26] Mapping left_kept_reads_seg2 to genome genome with Bowtie2 (2/2)
[2019-04-01 02:09:29] Mapping right_kept_reads to genome genome with Bowtie2
[2019-04-01 02:10:13] Mapping right_kept_reads_seg1 to genome genome with Bowtie2 (1/2)
[2019-04-01 02:10:15] Mapping right_kept_reads_seg2 to genome genome with Bowtie2 (2/2)
[2019-04-01 02:10:20] Searching for junctions via segment mapping Coverage-search algorithm is turned on, making this step very slow Please try running TopHat again with the option (–no-coverage-search) if this step takes too much time or memory.
[2019-04-01 02:10:34] Retrieving sequences for splices
[2019-04-01 02:10:35] Indexing splices
Building a SMALL index
[2019-04-01 02:10:35] Mapping left_kept_reads_seg1 to genome segment_juncs with Bowtie2 (1/2)
[2019-04-01 02:10:36] Mapping left_kept_reads_seg2 to genome segment_juncs with Bowtie2 (2/2)
[2019-04-01 02:10:37] Joining segment hits
[2019-04-01 02:10:40] Mapping right_kept_reads_seg1 to genome segment_juncs with Bowtie2 (1/2)
[2019-04-01 02:10:41] Mapping right_kept_reads_seg2 to genome segment_juncs with Bowtie2 (2/2)
[2019-04-01 02:10:43] Joining segment hits
[2019-04-01 02:10:46] Reporting output tracks
[2019-04-01 02:12:02] A summary of the alignment counts can be found in tophat.Sp_ds.dir/align_summary.txt
[2019-04-01 02:12:02] Run complete: 00:03:42 elapsed
Simlilar runs for other paired samples
Additional directories and contents following the run
~/tophat_spombe$ ls -l
drwxrwxr-x 3 ubuntu ubuntu 4096 Apr 1 02:12 tophat.Sp_ds.dir
drwxrwxr-x 3 ubuntu ubuntu 4096 Apr 1 02:17 tophat.Sp_hs.dir
drwxrwxr-x 3 ubuntu ubuntu 4096 Apr 1 02:22 tophat.Sp_log.dir
drwxrwxr-x 3 ubuntu ubuntu 4096 Apr 1 02:28 tophat.Sp_plat.dir
Cufflinks: Assemble expressed genes and transcripts
~/tophat_spombe$ cufflinks -p 8 -o clinks.sp_ds tophat.Sp_ds.dir/accepted_hits.bam
Warning: Could not connect to update server to verify current version. Please check at the Cufflinks website (http://cufflinks.cbcb.umd.edu).
[17:21:23] Inspecting reads and determining fragment length distribution.
Processed 5798 loci. [*************************] 100%
Map Properties:
Normalized Map Mass: 958393.49
Raw Map Mass: 958393.49
Fragment Length Distribution: Empirical (learned)
Estimated Mean: 263.91
Estimated Std Dev: 60.97
[17:21:29] Assembling transcripts and estimating abundances.
Processed 5798 loci. [*************************] 100%
####~/tophat_spombe$ cufflinks -p 8 -o clinks.sp_hs tophat.Sp_hs.dir/accepted_hits.bam ……………..
[17:29:59] Assembling transcripts and estimating abundances.
Processed 6381 loci. [*************************] 100%
~/tophat_spombe$ cufflinks -p 8 -o clinks.sp_log tophat.Sp_log.dir/accepted_hits.bam
………………………
[17:31:01] Assembling transcripts and estimating abundances.
Processed 5889 loci. [*************************] 100%
~/tophat_spombe$ cufflinks -p 8 -o clinks.sp_plat tophat.Sp_plat.dir/accepted_hits.bam
………………….
[17:31:59] Assembling transcripts and estimating abundances.
Processed 5699 loci. [*************************] 100%
Additional directories and contents following the run
~/tophat_spombe$ ls -ldrwxrwxr-x 2 ubuntu ubuntu 4096 Apr 9 17:21 clinks.sp_ds
drwxrwxr-x 2 ubuntu ubuntu 4096 Apr 9 17:29 clinks.sp_hs
drwxrwxr-x 2 ubuntu ubuntu 4096 Apr 9 17:31 clinks.sp_log
drwxrwxr-x 2 ubuntu ubuntu 4096 Apr 9 17:31 clinks.sp_plat
~/tophat_spombe$ ls -l clinks.sp_ds/
-rw-rw-r– 1 ubuntu ubuntu 351199 Apr 9 17:21 genes.fpkm_tracking
-rw-rw-r– 1 ubuntu ubuntu 407649 Apr 9 17:21 isoforms.fpkm_tracking
-rw-rw-r– 1 ubuntu ubuntu 0 Apr 9 17:21 skipped.gtf
-rw-rw-r– 1 ubuntu ubuntu 2432834 Apr 9 17:21 transcripts.gtf
~/tophat_spombe$ ls clinks.sp_ds/transcripts.gtf > assembly.txt
~/tophat_spombe$ ls clinks.sp_hs/transcripts.gtf >> assembly.txt
~/tophat_spombe$ ls clinks.sp_log/transcripts.gtf >> assembly.txt
~/tophat_spombe$ ls clinks.sp_plat/transcripts.gtf >> assembly.txt
~/tophat_spombe$ cat assembly.txt
clinks.sp_ds/transcripts.gtf
clinks.sp_hs/transcripts.gtf
clinks.sp_log/transcripts.gtf
clinks.sp_plat/transcripts.gtf
Cuffmerge: Run Cuffmerge on all assemblies to create a single merged transcriptome annotation
~/tophat_spombe$ cuffmerge -g Schizosaccharomyces_pombe.ASM294v2.42.gtf -s Schizosaccharomyces_pombe_all_chromosomes.fa -p 8 assembly.txt
[Tue Apr 9 17:44:48 2019] Beginning transcriptome assembly merge
[Tue Apr 9 17:44:48 2019] Preparing output location ./merged_asm/
[Tue Apr 9 17:44:49 2019] Converting GTF files to SAM
[17:44:49] Loading reference annotation.
[17:44:49] Loading reference annotation.
[17:44:49] Loading reference annotation.
[17:44:49] Loading reference annotation.
[Tue Apr 9 17:44:49 2019] Quantitating transcripts
Warning: Could not connect to update server to verify current version. Please check at the Cufflinks website (http://cufflinks.cbcb.umd.edu).
Command line:
cufflinks -o ./merged_asm/ -F 0.05 -g Schizosaccharomyces_pombe.ASM294v2.42.gtf -q –overhang-tolerance 200 –library-type=transfrags -A 0.0 –min-frags-per-transfrag 0 –no-5-extend -p 8 ./merged_asm/tmp/mergeSam_file56TGqQ
[bam_header_read] EOF marker is absent. The input is probably truncated.
[bam_header_read] invalid BAM binary header (this is not a BAM file).
File ./merged_asm/tmp/mergeSam_file56TGqQ doesn’t appear to be a valid BAM file, trying SAM…
[17:44:49] Loading reference annotation.
[17:44:49] Inspecting reads and determining fragment length distribution.
Processed 3400 loci.
Map Properties:
Normalized Map Mass: 18656.00
Raw Map Mass: 18656.00
Fragment Length Distribution: Truncated Gaussian (default)
Default Mean: 200
Default Std Dev: 80
[17:44:49] Assembling transcripts and estimating abundances.
Processed 3400 loci.
[Tue Apr 9 17:44:52 2019] Comparing against reference file Schizosaccharomyces_pombe.ASM294v2.42.gtf
Warning: Could not connect to update server to verify current version. Please check at the Cufflinks website (http://cufflinks.cbcb.umd.edu).
No fasta index found for Schizosaccharomyces_pombe_all_chromosomes.fa. Rebuilding, please wait..
Fasta index rebuilt.
Warning: couldn’t find fasta record for ‘AB325691’!
Warning: couldn’t find fasta record for ‘MT’!
Warning: couldn’t find fasta record for ‘MTR’!
………………..
Additional directories and contents following the run
~/tophat_spombe$ ls -lt
drwxrwxr-x 3 ubuntu ubuntu 4096 Apr 9 17:44 merged_asm
-rw-rw-r– 1 ubuntu ubuntu 185 Apr 9 17:44 Schizosaccharomyces_pombe_all_chromosomes.fa.fai
-rw-rw-r– 1 ubuntu ubuntu 119 Apr 9 17:41 assembly.txt
Cuffdiff: Run Cuffdiff by using the merged transcriptome assembly along with the BAM files from TopHat for each replicate
~/tophat_spombe$ cuffdiff -o diff_out -b Schizosaccharomyces_pombe_all_chromosomes.fa -p 8 -L C1,C2 -u merged_asm/merged.gtf tophat.Sp_ds.dir/accepted_hits.bam,tophat.Sp_hs.dir/accepted_hits.bam tophat.Sp_log.dir/accepted_hits.bam,tophat.Sp_plat.dir/accepted_hits.bam
** Here assuming Sp_ds and Sp_hs are replicates. Similarly it is assumed that Sp-log and Sp_plat are replicates.**
Warning: Could not connect to update server to verify current version. Please check at the Cufflinks website (http://cufflinks.cbcb.umd.edu).
[17:53:51] Loading reference annotation and sequence.
Warning: couldn’t find fasta record for ‘AB325691’!
This contig will not be bias corrected.
Warning: couldn’t find fasta record for ‘MT’!
This contig will not be bias corrected.
Warning: couldn’t find fasta record for ‘MTR’!
This contig will not be bias corrected.
[17:53:52] Inspecting maps and determining fragment length distributions.
[17:54:07] Modeling fragment count overdispersion.
[17:54:58] Modeling fragment count overdispersion.
Map Properties: Normalized Map Mass: 933414.68 Raw Map Mass: 922901.69 Number of Multi-Reads: 35945 (with 148785 total hits) Fragment Length Distribution: Empirical (learned) Estimated Mean: 262.47 Estimated Std Dev: 61.71
Map Properties: Normalized Map Mass: 933414.68 Raw Map Mass: 923893.70 Number of Multi-Reads: 64738 (with 379511 total hits) Fragment Length Distribution: Empirical (learned) Estimated Mean: 276.70 Estimated Std Dev: 60.73
Map Properties: Normalized Map Mass: 933414.68 Raw Map Mass: 930065.32 Number of Multi-Reads: 40694 (with 116887 total hits) Fragment Length Distribution: Empirical (learned) Estimated Mean: 263.99 Estimated Std Dev: 60.37
Map Properties: Normalized Map Mass: 933414.68 Raw Map Mass: 927012.86 Number of Multi-Reads: 48144 (with 294370 total hits) Fragment Length Distribution: Empirical (learned) Estimated Mean: 264.10 Estimated Std Dev: 64.51
[17:55:50] Calculating preliminary abundance estimates
Processed 3744 loci. [*************************] 100%
[17:56:26] Learning bias parameters.
[17:56:41] Testing for differential expression and regulation in locus.
Processed 3744 loci. [*************************] 100%
Performed 5100 isoform-level transcription difference tests
Performed 4916 tss-level transcription difference tests
Performed 4431 gene-level transcription difference tests
Performed 4383 CDS-level transcription difference tests
Performed 0 splicing tests
Performed 0 promoter preference tests
Performing 0 relative CDS output tests
Writing isoform-level FPKM tracking
Writing TSS group-level FPKM tracking
Writing gene-level FPKM tracking
Writing CDS-level FPKM tracking
Writing isoform-level count tracking
Writing TSS group-level count tracking
Writing gene-level count tracking
Writing CDS-level count tracking
Writing isoform-level read group tracking
Writing TSS group-level read group tracking
Writing gene-level read group tracking
Writing CDS-level read group tracking
Writing read group info
Writing run info
cummeRbund: Data analysis
library(cummeRbund)Create a CummeRbund database from the Cuffdiff output
cuff_data <- readCufflinks("diff_out/")
cuff_dataCuffSet instance with:
2 samples
6016 genes
7690 isoforms
6997 TSS
5121 CDS
6016 promoters
6997 splicing
4575 relCDSPlot the distribution of expression levels for each sample
csDensity(genes(cuff_data))
Compare the expression of each gene in two conditions with a scatter plot
csScatter(genes(cuff_data), 'C1', 'C2')
Create a volcano plot to inspect differentially expressed genes
csVolcano(genes(cuff_data), 'C1', 'C2')
Genes and transcripts that are differentially expressed between two samples.
gene_diff_data <- diffData(genes(cuff_data))
head(gene_diff_data)gene_diff_data <- diffData(genes(cuff_data))
sig_gene_data <- subset(gene_diff_data, (significant == 'yes'))
head(sig_gene_data)Differentially expressed transcripts or differentially spliced and regulated genes
isoform_diff_data <- diffData(isoforms(cuff_data), 'C1', 'C2')
sig_isoform_data <- subset(isoform_diff_data, (significant == 'yes'))
nrow(sig_isoform_data)[1] 0tss_diff_data <- diffData(TSS(cuff_data), 'C1', 'C2')
sig_tss_data <- subset(tss_diff_data, (significant == 'yes'))
nrow(sig_tss_data)[1] 0cds_diff_data <- diffData(CDS(cuff_data), 'C1', 'C2')
sig_cds_data <- subset(cds_diff_data, (significant == 'yes'))
nrow(sig_cds_data)[1] 0promoter_diff_data <- distValues(promoters(cuff_data))
sig_promoter_data <- subset(promoter_diff_data, (significant == 'yes'))
nrow(sig_promoter_data)[1] 0splicing_diff_data <- distValues(splicing(cuff_data))
sig_splicing_data <- subset(splicing_diff_data, (significant == 'yes'))
nrow(sig_splicing_data)[1] 0relCDS_diff_data <- distValues(relCDS(cuff_data))
sig_relCDS_data <- subset(relCDS_diff_data, (significant == 'yes'))
nrow(sig_relCDS_data)[1] 0