Transcript Document

TUXEDO suite:
BOWTIE + TOPHAT + CUFFLINKS + CUMMERBUND
Enrique Blanco
2013
0. SOURCES
1. RAW FORMAT: FASTQ
Example:
@1-ILLUMINA-GA:90515:1:1:28:2023
GGAGAAGTAGCCCGGAATGGTCGAGTAGGGAAGAGAAAGT
+
a^Z]HY]]SZaV]a`VL\OOWNYT_\`BBBBBBBBBBBBB
@1-ILLUMINA-GA:90515:1:1:28:1470
AAGGGCGCGAAGCCCGTTACCAAGCCCGCCAAGGGAACCG
+
QYFQR__TXNUU]U]VGFFOOZTFZT_`]HNVJV\G__BB
@1-ILLUMINA-GA:90515:1:1:28:2017
CGATACGTTCGCGTTCGTGTAACTTTGTGTTTGCAAATCA
+
aa]a^aa_aa``^_aaaaaa_Z`baaa`a```aa^`^a^T
@1-ILLUMINA-GA:90515:1:1:28:931
ACTGCCCCACATCAGCACACCAAACACGTCCACCGCAGCT
+
INMG[TXXZYMFPTXN]NN]NZ\X_Z_ZN[XYBBBBBBBB
…
2. INDEXES FOR READ MAPPING WITH TOPHAT (BOWTIE)
Drosophila_melanogaster/UCSC/dm3/Annotation/Archives/archive-current/Genes/
genes.gtf
refGene.txt
Drosophila_melanogaster/UCSC/dm3/Sequence/BowtieIndex/
genome.1.ebwt genome.2.ebwt genome.3.ebwt
genome.rev.1.ebwt genome.rev.2.ebwt
genome.4.ebwt
Drosophila_melanogaster/UCSC/dm3/Sequence/WholeGenomeFasta/
genome.fa
Copy only those files into a single folder to export it:
export FLY="/usr/local/molbio/indexes/dm3"
genome.fa.fai
3. FROM FASTQ TO READ MAPS IN dm3: RNAseq single-end
2 samples (WID, EID), 2 replicates (R1, R2)
> tophat -p 4 -G $FLY/genes.gtf -o WID_R1_tophat $FLY/genome rawfiles/WID_R1.fq
> tophat -p 4 -G $FLY/genes.gtf -o WID_R2_tophat $FLY/genome rawfiles/WID_R2.fq
> tophat -p 4 -G $FLY/genes.gtf -o EID_R1_tophat $FLY/genome rawfiles/EID_R1.fq
> tophat -p 4 -G $FLY/genes.gtf -o EID_R2_tophat $FLY/genome rawfiles/EID_R2.fq
Number of processes (threads)
Collection of known
genes (RefSeq) to
map reads on them
Indexes for read mapping (dm3)
Output folder
WID_R1_tophat/
accepted_hits.bam
deletions.bed
insertions.bed
junctions.bed
logs/
prep_reads.info
unmapped.bam
mapped read in BAM format (binary)
Input sample
(for pair-end, two files)
4. SAM/BAM format: READ MAPS
SAM example
1-ILLUMINA-GA:90515:1:88:1105:333
0
CGCCAGTATACACGCAGCACCAGGTGCGCAATGAAGCCCC
(01) 1-ILLUMINA-GA:90515:1:88:1105:333
(02) 0
(03) chr2L
(04) 7732
50
40M
*
0
0
(10) CGCCAGTATACACGCAGCACCAGGTGCGCAATGAAGCCCC
a```````a_\a`][^`Z_P]```Z^QDT_``\\\VS\]\
XA:i:0 MD:Z:40 NM:i:0 XS:A:+ NH:i:1
chr2L
7732
50
40M
*
0
a```````a_\a`][^`Z_P]```Z^QDT_``\\\VS\]\
<-----<-----<-----<-----<-----<------
0
XA:i:0
MD:Z:40 NM:i:0
XS:A:+
NAME
FLAG: (TOPHAT) 0 means mapped in FWD; 16 means mapped in RVS
CHROMOSOME
1-based POSITION IN CHROMOSOME
MAPPING quality
CIGAR alignment
<------ SEQUENCE
<------ PHRED quality score
SAM to BAM (compression)
> samtools view -bS -o ec_snp.bam ec_snp.sam
-rw-rw-r-- 1 eblanco eblanco 321M Jun 25 12:43 accepted_hits.bam
-rw-rw-r-- 1 eblanco eblanco 1.4G Jul 2 12:20 accepted_hits.sam
NH:i:1
5. FROM READ MAPS TO GENE EXPRESSION QUANTIFICATION IN RPKMs
> cufflinks -p 4 -G $FLY/genes.gtf -o WID_R1_cufflinks WID_R1_tophat/accepted_hits.bam
> cufflinks -p 4 -G $FLY/genes.gtf -o WID_R2_cufflinks WID_R2_tophat/accepted_hits.bam
> cufflinks -p 4 -G $FLY/genes.gtf -o EID_R1_cufflinks EID_R1_tophat/accepted_hits.bam
> cufflinks -p 4 -G $FLY/genes.gtf -o EID_R2_cufflinks EID_R2_tophat/accepted_hits.bam
Output folder
Read maps (bam)
genes.fpkm_tracking:
tracking_id
cbt
vg
brk
…
X
-
X
-
gene_id gene_short tss_id
cbt
cbt
TSS16012
vg
vg
TSS14634
brk
brk
TSS16005
locus
chr2L:476445-479688
chr2R:8771705-8786899
chrX:7201973-7205165
X
-
X
-
FPKM
40.8326
45.6418
32.7855
FPKM_conf_lo
37.0929
41.3073
28.8919
FPKM_conf_hi
44.6059
48.0142
35.5787
FPKM_status
OK
OK
OK
isoforms.fpkm_tracking:
tracking_id
NM_176558
NM_001260363
NM_170160
…
X
-
X
-
gene_id gene_short tss_id
locus
length
ash2
ash2
TSS10523 chr3R:20475809-20479098 1986
ash2
ash2
TSS16687 chr3R:20476901-20479098 1271
ash2
ash2
TSS13304 chr3R:20477138-20479098 1413
X
-
FPKM
36.7716
4.4918
10.4132
FPKM_conf_lo
31.0929
2.64044
7.12527
FPKM_conf_hi
39.6772
6.4426
13.5855
FPKM_status
OK
OK
OK
6. MERGE SAMPLES + DIFFERENTIAL EXPRESSION ANALYSIS
> cuffmerge -g $FLY/genes.gtf -s $FLY/genome.fa -p 8 -o merged_transcriptomes data/assemblies.txt
> more data/assemblies.txt
./WID_R1_cufflinks/transcripts.gtf
./WID_R2_cufflinks/transcripts.gtf
./EID_R1_cufflinks/transcripts.gtf
./EID_R2_cufflinks/transcripts.gtf
Uniform exon annotation
> cuffdiff -o DE_info -b $FLY/genome.fa -p 8 -L WID,EID -u merged_transcriptomes/merged.gtf
./WID_R1_tophat/accepted_hits.bam,./WID_R2_tophat/accepted_hits.bam
./EID_R1_tophat/accepted_hits.bam,./EID_R2_tophat/accepted_hits.bam
Sample 1: WID, 2 replicates
Output folder:
bias_params.info
cds.count_tracking
cds.diff
cds_exp.diff
cds.fpkm_tracking
cds.read_group_tracking
gene_exp.diff
genes.count_tracking
genes.fpkm_tracking
genes.read_group_tracking
isoform_exp.diff
isoforms.count_tracking
isoforms.fpkm_tracking
isoforms.read_group_tracking
promoters.diff
read_groups.info
run.info
splicing.diff
tss_group_exp.diff
tss_groups.count_tracking
tss_groups.fpkm_tracking
tss_groups.read_group_tracking
var_model.info
Sample 2: EID, 2 replicates
IMPORTANT
FILES!!!
7. CUFFDIFF FILES
gene_exp.diff
test_id
gene_id
gene
locus
sample_1 sample_2 status
value_1 value_2 log2(FC)
XLOC_004259
XLOC_003297
…
XLOC_012036
XLOC_002876
…
XLOC_004259
XLOC_003297
ap
vg
chr2R:1593706-1615030
chr2R:8771705-8786899
WID
WID
EID
EID
OK
OK
122.969
54.965
1.817
0.563
XLOC_012036
XLOC_002876
ey
so
chr4:718314-741799
chr2R:3306854-3322377
WID
WID
EID
EID
OK
OK
0.086
9.326
40.260
44.971
test_stat
-6.07986 -23.4656
-6.60692 -18.6605
8.85814
2.26965
4.81798
12.8272
p_value q_value
significant
5e-05
5e-05
0.000146346
0.000146346
yes
yes
0.0023
5e-05
0.00504939
0.000146346
yes
yes
isoform_exp.diff
test_id
gene_id
gene
locus
sample_1 sample_2 status
TCONS_00032967 XLOC_014395 pcm
chrX:19379912-19385984 WID
EID
OK
…
TCONS_00032972 XLOC_014397 l(1)G0156 chrX:19412460-19415843 WID
EID
OK
…
value_1 value_2 log2(FC) test_stat
20.9259 10.3616 -1.014 -3.01467
2.86149 41.4922
3.858
4.43662
p_value q_value significant
5e-05 0.000287866 yes
5e-05 0.000287866
yes
8. CUMMERBUND (R): Global plots
(requires R v3.0)
source("http://bioconductor.org/biocLite.R")
biocLite()
biocLite('cummeRbund')
library("cummeRbund")
# load the cuff_data
cuff <- readCufflinks()
pdf("../plots/densityGenes.pdf")
dens<-csDensity(genes(cuff))
dens
dev.off()
pdf("../plots/densityGenesRep.pdf")
densRep<-csDensity(genes(cuff),replicates=T)
densRep
dev.off()
pdf("../plots/scatter1Plot.pdf")
s1<-csScatter(genes(cuff),"WID","EID",smooth=T)
s1
dev.off()
9. CUMMERBUND (R): One gene
mygene <- getGene(cuff,"ewg")
genebarplot <- expressionBarplot(mygene)
pdf("../plots/geneBarplotEWG.pdf")
genebarplot
dev.off()
pdf("../plots/isoBarplotEWG.pdf")
isobarplot <- expressionBarplot(isoforms(mygene))
isobarplot
dev.off()
10. CUMMERBUND (R): Groups of genes
mydegenesID <- getSig(cuff,alpha=0.05,level='genes')
or
mygenesNames <-c("ey","ato","Appl","ey","eya","onecut","Optix","run","sev")
mygenes<-getGenes(cuff,mygenesNames)
csHeatmap(mygenes,cluster='both',replicates=F)
11. CUMMERBUND (R): Mixing groups and graphics
mygenesNames <-c("arm","ewg","NrxIV","chn","ppn","siz")
mygenes<-getGenes(cuff,mygenesNames)
b<-expressionBarplot(mygenes)
pdf("../plots/expressionBarplotSET.pdf")
b
dev.off()
c<-expressionBarplot(isoforms(mygenes))
pdf("../plots/ISOexpressionBarplotSET.pdf")
c
dev.off()