Transcript J`ai fait séquencer mes petits ARN. Et maintenant

Introduction à l’analyse des données de
séquençage à haut débit en génomique
fonctionnelle.
28 mars 2012, 15:30 – 17:00
J’ai fait séquencer mes petits ARN.
Et Maintenant ?
[email protected]
http://drosophile.org
Les trois principales classes de
petits ARNs chez la drosophile
Hen1
met
+
Produits des snoRNA, tRNA, rRNA.
2S Droso (30nt)
small RNA deep sequencing
(Biases)
20-30nt RNA gel purification
Library “Bar coding”
Que Puis-je Faire avec mes séquences de
petits ARN ?

Annotation

Visualisation

Découverte de loci

Quantification d’expression

Analyse structurale des précurseurs, signatures, …

Mise en évidence de « visiteurs » (virus, …)

…
Informatique  Bioinformatique
Matériel

Un fichier de séquence au format fastq

Un ordinateur avec ~ 8 Mo RAM

Un « Operating System Unix compliant »

Un maniement confortable de cet OS

Quelques logiciels génériques très utiles

Un « vrai » éditeur de texte (TextWrangler, etc..)

R, Gnuplot

…

Une bonne connaissance du web

Le maniement niveau Débutant++ d’un langage de programmation

Perl

Python
Que contient le gros fichier fastq que j’ai téléchargé (et décompressé)
?
* Limite max pour ouvrir un gros fichier texte (~1.2 Go)
 Terminal Unix. Naviguer dans le dossier qui contient le fichier
 Taper la commande more <nom_du_fichier>
lbcd-05:GKG13demo deepseq$ more GKG-13.fastq
@HWIEAS210R_0028:2:1:3019:1114#AGAAGA/1
TNGGAACTTCATACCGTGCTCTCTGTAGGCACCATCAA
+HWIEAS210R_0028:2:1:3019:1114#AGAAGA/1
bBb`bfffffhhhhhhhhhhhhhhhhhhhfhhhhhhgh
@HWIEAS210R_0028:2:1:3925:1114#AGAAGA/1
TNCTTGGACTACATATGGTTGAGGGTTGTACTGTAGGC
+HWIEAS210R_0028:2:1:3925:1114#AGAAGA/1
]B]VWaaaaaagggfggggggcggggegdgfgeggbab
@HWIEAS210R_0028:2:1:6220:1114#AGAAGA/1
TNGGAACTTCATACCGTGCTCTCTGTAGGCACCATCAA
+HWIEAS210R_0028:2:1:6220:1114#AGAAGA/1
aB^^afffffhhhhhhhhhhhhhhhhhhhhhhhchhhh
@HWIEAS210R_0028:2:1:6252:1115#AGAAGA/1
TNCTTGGACTACATATGGTTGAGGGTTGTACTGTAGGC
+HWIEAS210R_0028:2:1:6252:1115#AGAAGA/1
aBa^\ddeeehhhhhhhhhhhhhhhhghhhhhhhefff
@HWIEAS210R_0028:2:1:6534:1114#AGAAGA/1
TNAATGCACTATCTGGTACGACTGTAGGCACCATCAAT
+HWIEAS210R_0028:2:1:6534:1114#AGAAGA/1
aB\^^eeeeegcggfffffffcfffgcgcfffffR^^]
@HWIEAS210R_0028:2:1:8869:1114#AGAAGA/1
GNGGACTGAAGTGGAGCTGTAGGCACCATCAATAGATC
+HWIEAS210R_0028:2:1:8869:1114#AGAAGA/1
aBaaaeeeeehhhhhhhhhhhhfgfhhgfhhhhgga^^
…
…
…
Combien de séquences dans mon fichier ?
 Terminal Unix. Naviguer dans le dossier qui contient le fichier
 Taper la commande wc - l <nom_du_fichier>
lbcd-05:GKG13demo deepseq$ wc -l
GKG-13.fastq
25703828 GKG-13.fastq
>>> 25 703 828 / 4
6 425 957 séquences
Mes séquences contiennent-elles le bon adaptateur ?
Séquence de mon adaptateur: CTGTAGGCACCATCAAT
 Taper la commande cat <nom_du_fichier> | grep CTGTAGG | wc -l
lbcd-05:GKG13demo deepseq$ wc -l
GKG-13.fastq | grep CTGTAGG | wc -l
6355061
6 355 061 sur
6 425 957 séquences
Pas mal
A contrario
lbcd-05:GKG13demo deepseq$ wc -l
308
GKG-13.fastq | grep ATCTCGT| wc -l
Quelle est la qualité de mes séquences ?
http://www.bioinformatics.babraham.ac.uk/projects/fastqc/
Comment retirer l’adaptateur ?
http://hannonlab.cshl.edu/fastx_toolkit/index.html
Séquence de mon adaptateur: CTGTAGGCACCATCAAT
deepseq$ fastq_to_fasta -r –n -i GKG-13.fastq -o GKG-13.fasta
deepseq$ more GKG-13.fasta
>1
AATGGCACTGGAAGAATTCACCTGTAGGCACCATCAAT
>2
TCTCGGTAGAACCTCCACTGTAGGCACCATCAATAGAT
>3
deepseq$ fastx_clipper -a CTGTAGGCACCATCAAT -l 18 -i GKG-13.fasta
-o GKG-13_clipped.fasta
TTTGTGACCGACACTAACGGGTACTGTAGGCACCATCA
>4
TGGAATGTAAAGAAGTATGGAGCTGTAGGCACCATCAA
>5
>18
GTCAGCAACTTGATTCCAGCAATCTGTAGGCACCATCA
AATGGCACTGGAAGAATTCAC
>6
>20
deepseq$ more GKG-13_clipped.fasta
AATGGCACTGGAAGAATTCACGGGCTGTAGGCACCATC
TTTGTGACCGACACTAACGGGTA
>7
>21
TGGAAGACTAGTGATTTTGTTCTGTAGGCACCATCAAT
TGGAATGTAAAGAAGTATGGAG
>8
>22
TGAACACAGCTGGTGGTATCCCTGTAGGCACCATCAAT
GTCAGCAACTTGATTCCAGCAAT
>23
AATGGCACTGGAAGAATTCACGGG
>24
TGGAAGACTAGTGATTTTGTT
>25
TGAACACAGCTGGTGGTATCC
deepseq$ fastq_to_fasta -r -n -i GKG-13.fastq |
>26
TAAGTACTAGTGCCGCAGGA
>27
fastx_clipper -a CTGTAGGCACCATCAAT -l 18 -o GKG-13_clip-pipe.fasta
TGAACACAGCTGGTGGTATC
>28
TAGGAACTTCATACCGTGCTCT
J’utilise fastx_clipper et fastQC pour visualiser la distribution de taille
de mes séquences
deepseq$ fastx_clipper -a CTGTAGGCACCATCAAT -l 0 -i GKG-13.fastq -o GKG-13_clipped.fastq
deepseq$ more GKG-13_clipped.fastq
@HWIEAS210R_0028:2:1:1313:1120#AGAAGA/1
AATGGCACTGGAAGAATTCAC
+HWIEAS210R_0028:2:1:1313:1120#AGAAGA/1
fe\gggd\fgeeeggdaggag
@HWIEAS210R_0028:2:1:1387:1119#AGAAGA/1
TCTCGGTAGAACCTCCA
+HWIEAS210R_0028:2:1:1387:1119#AGAAGA/1
gggggeggfffgggfff
@HWIEAS210R_0028:2:1:1849:1120#AGAAGA/1
TTTGTGACCGACACTAACGGGTA
+HWIEAS210R_0028:2:1:1849:1120#AGAAGA/1
hhhhhhhhhfhgfhhhhgehhha
http://bowtie-bio.sourceforge.net/
Bowtie aligne des reads sur un génome
de référence préalablement préparé
Je télécharge Bowtie, je l’installe, et je lis le manuel
Je télécharge mon génome au format FASTA
Je prépare mon « index » Bowtie
deepseq$ bowtie-build fasta_libraries/dmel-all-chromosome-r5.37.fasta dmel-r5.37
~5 min
deepseq$ ls –laht
-rw-r--r-- 1 deepseq
-rw-r--r-- 1 deepseq
-rw-r--r-- 1 deepseq
-rw-r--r-- 1 deepseq
-rw-r--r-- 1 deepseq
-rw-r--r-- 1 deepseq
staff
staff
staff
staff
staff
staff
49M Mar 24 17:24 dmel-r5.37.rev.1.ebwt
19M Mar 24 17:24 dmel-r5.37.rev.2.ebwt
49M Mar 24 17:20 dmel-r5.37.1.ebwt
19M Mar 24 17:20 dmel-r5.37.2.ebwt
331K Mar 24 17:16 dmel-r5.37.3.ebwt
39M Mar 24 17:16 dmel-r5.37.4.ebwt
J’aligne mes reads avec bowtie
deepseq$ bowtie ~/bin/bowtie/indexes/5.43_Dmel/5.43_Dmel -f GKG-13_clip-pipe.fasta -v 1 -k 1 -p 12 --al
droso_matched_GKG-13.fa --un unmatched_GKG13.fa > GKG13_bowtie_output.tabulated
~/bin/bowtie/indexes/5.43_Dmel/5.43_Dmel
-f GKG-13_clip-pipe.fasta
-v 1
-k 1
-p 12
--al droso_matched_GKG-13.fa
--un unmatched_GKG13.fa
> GKG13_bowtie_output.tabulated
# reads processed: 5997502
# reads with at least one reported alignment: 5045151 (84.12%)
# reads that failed to align: 952351 (15.88%)
Reported 5045151 alignments to 1 output stream(s)
… et je récupère
deepseq$ ls -laht
-rw-r--r-- 1 deepseq staff 351M Mar 24 17:46 GKG13_bowtie_output.tabulated
-rw-r--r-- 1 deepseq staff 156M Mar 24 17:46 droso_matched_GKG-13.fa
-rw-r--r-- 1 deepseq staff 28M Mar 24 17:46 unmatched_GKG13.fa
Un fichier d’alignement
deepseq$ more GKG13_bowtie_output.tabulated
21
+
2L
20487495
TGGAATGTAAAGAAGTATGGAG
30
3L
15836559
GTGAATTCTCCCAGTGCCAAG
25
+
3R
5916902 TGAACACAGCTGGTGGTATCC
23
2L
11953462
CCCGTGAATTCTTCCAGTGCCATT
27
+
3R
5916902 TGAACACAGCTGGTGGTATC
26
3R
9289997 TCCTGCGGCACTAGTACTTA
18
2L
11953465
GTGAATTCTTCCAGTGCCATT
22
3R
8377246 ATTGCTGGAATCAAGTTGCTGAC
20
+
3L
11650036
TTTGTGACCGACACTAACGGGTA
24
+
2R
16493585
TGGAAGACTAGTGATTTTGTT
28
+
3L
10358380
TAGGAACTTCATACCGTGCTCT
35
+
X
18022302
CTTGTGCGTGTGACAGCGGCT
41
3RHet 138608 TGGCGACCGTGACAGGACCCG
42
+
3R
5916902 TGAACACAGCTGGTGGTATCC
Un fichier des séquences alignées
Un fichier des séquences non alignées
deepseq$ more droso_matched_GKG-13.fa
>21
TGGAATGTAAAGAAGTATGGAG
>26
TAAGTACTAGTGCCGCAGGA
>24
TGGAAGACTAGTGATTTTGTT
>23
AATGGCACTGGAAGAATTCACGGG
>27
TGAACACAGCTGGTGGTATC
deepseq$ more unmatched_GKG13.fa
>29
AGGGGGCTATTTCACTACTGGA
>33
CGATGATGACGGTACCCGTAGA
>37
GCTAGTCGGTACTTGAAAC
>59
TGGTTGCAATAGCTTCTGGCGGA
>61
GATGAGTGCTAGATGTAGGGA
Un pipeline d’annotations « génomiques »
Sequence reads (fasta format)
Bowtie
Pre-miRNAs (miRBase)
Unmatched reads
Bowtie
Non coding RNAs
Unmatched reads
Bowtie
Transposons
Matched reads
(fasta)
Read Count
Matched reads
(fasta)
Read Count
hierarchical
Matched reads
(fasta)
Read Count
Matched reads
(fasta)
Read Count
annotation
Unmatched reads
Bowtie
Genes
of
sequence
Unmatched reads
Bowtie
Matched reads
(fasta)
Read Count
Viruses, transgenes, etc…Matched reads
(fasta)
Read Count
Intergenic regions
Unmatched reads
Bowtie
Remaining unmatched sequences
datasets
Je veux visualiser mes reads dans un « Genome Browser »
http://samtools.sourceforge.net/
Un pipeline sommaire pour préparer un fichier de visualisation
deepseq$ bowtie -v 1 -M 1 --best /Users/deepseq/bin/bowtie/indexes/5.37_Dmel -p 12 -f GKG-13_clip-pipe.fasta -S | samtools view
-bS -o GKG-13_clip-pipe.fasta.bam - ; samtools sort GKG-13_clip-pipe.fasta.bam GKG-13_clip-pipe.fasta.bam.sorted ; samtools
index GKG-13_clip-pipe.fasta.bam.sorted.bam

306K GKG-13_clip-pipe.fasta.bam.sorted.bam.bai
42M GKG-13_clip-pipe.fasta.bam.sorted.bam
80M GKG-13_clip-pipe.fasta.bam
Je veux visualiser mes reads dans un « Genome Browser » (2)
J’upload mes fichiers bam et bai sur un serveur accessible
J’indique l’URL du fichier bam à Ensembl (Gbrowse, Modencode, etc..)
Je veux visualiser mes reads dans un « Genome Browser » (3)
Je navigue dans les régions d’intérêt, après avoir indiqué au Browser d’inclure mon « track
6.43 Kb
18,470,000
18,471,000
18,473,000
Chromosome bands
GKG13_demo
Forward strand
18,472,000
18,474,000
18,475,000
18,474,000
18,475,000
36F
162027
227692 features from 'GKG13_demo' omitted
FlyBase feature
mir-125-RM >
pre miRNA
mir-125-RA >
miRNA
mir-100-RM >
let-7-RM >
pre miRNA
pre miRNA
mir-100-RA >
let-7-RA >
miRNA
miRNA
CR43344-RA >
ncRNA
2L_370 >
Contigs
FlyBase feature
< CG10283-RC
protein coding
< CG10283-RB
protein coding
< CG10283-RA
protein coding
< CG10283-RD
protein coding
Unigene EST clus...
Fly cDNA
UniProtKB protein
%GC
18,470,000
18,471,000
Reverse strand
Gene Legend
protein coding
There are currently 34 tracks turned off.
Ensembl Drosophila melanogaster version 66.539 (BDGP5) Chromosome 2L: 18,469,145 - 18,475,575
18,472,000
18,473,000
6.43 Kb
RNA gene
Je veux visualiser mes reads dans un « Genome Browser » (4)
Encore un…
9.71 Kb
15,547,000
15,548,000
15,549,000
15,550,000
Forward strand
15,551,000
Chromosome bands
15,552,000
15,553,000
15,554,000
15,555,000
15,556,000
71E
84
GKG13_demo
65 features from 'GKG13_demo' omitted
UniProtKB protein
Fly cDNA
Unigene EST clus...
FlyBase feature
AGO2-RB >
protein coding
AGO2-RC >
protein coding
AGO2-RD >
protein coding
3L_311 >
Contigs
3L_312 >
FlyBase feature
< CG7739-RA
protein coding
Unigene EST clus...
Fly cDNA
Fly gold cDNA
437 bp
Forward strand
UniProtKB protein
15,553,800
%GC
15,553,850
15,553,900
15,553,950
15,554,000
Chromosome bands
GKG13_demo
15,554,050
15,547,000
15,548,000
15,549,000
15,550,000
Reverse strand
15,551,000
15,552,000
15,553,000
9.71 Kb
protein coding
There are currently 34 tracks turned off.
Ensembl Drosophila melanogaster version 66.539 (BDGP5) Chromosome 3L: 15,546,772 - 15,556,479
41 features from 'GKG13_demo' omitted
Fly cDNA
AGO2-RB >
protein coding
AGO2-RC >
protein coding
AGO2-RD >
protein coding
3L_312 >
Contigs
FlyBase feature
< CG7739-RA
protein coding
Unigene EST clus...
Fly cDNA
Fly gold cDNA
UniProtKB protein
%GC
15,554,150
71E
80
Gene Legend
FlyBase feature
15,554,100
15,554,000
15,555,000
15,556,000
15,554,200
Un “profiler” maison pour les micros ARNs
deepseq$ miRNA_bowtie_profiler.py GKG-13_clip-pipe.fasta ~/bin/bowtie/indexes/dme_miR_r17.1.ebwt
Sequence reads (fasta format)
Bowtie
Pre-miRNAs (miRBase)
Indéxé pour Bowtie
Bowtie Output
Analyse “textuelle”


Cartes des reads par miRNA
Liste de comptage par miR_5p et miR_3p
# bowtie -v 1 -M 1 --best --strata -p 12 --norc --suppress 2,6,7,8 /Users/deepseq/bin/bowtie/indexes/dme_miR_r17 -f GKG-13_clip-pipe.fasta
# reads processed: 5997502
# reads with at least one reported alignment: 3886779 (64.81%)
# reads that failed to align: 2060565 (34.36%)
# reads with alignments sampled due to -M: 50158 (0.84%)
Reported 3886779 alignments to 1 output stream(s)
# Parsing completed in 1 minutes and 36.7 seconds
sizes
offsets
counts
miRNA_bowtie_profiler.py : Cartes des reads, par miR
miRNA_bowtie_profiler.py : Attribution des reads “5p” et “3p”
miRs « 5p »
987 reads
+
*
miRs « 3p »
16003 reads
= 16990, ~ 17009 reads
miRNA_bowtie_profiler.py : Liste de comptage des miRs
..
.
..
.
.
..
.
..
Analyse d’expression différentielle
deepseq$ miRNA_bowtie_profiler.py GKG-13_clip-pipe.fasta ~/bin/bowtie/indexes/dme_miR_r17.1.ebwt
Sequence reads (fasta format)
Bowtie
Pre-miRNAs (miRBase)
Indéxé pour Bowtie
Bowtie Output
Analyse “textuelle”


Cartes des reads par miRNA
Liste de comptage par miR_5p et miR_3p
http://www.r-project.org/
DESeq
Heatplus
edgeR
(Bioconductor)
Profiling des miRNAs durant
la métamorphose de la
drosophile
L3
PF PF+12h
Ecdysone titer
days
Read
count table
Clustering of miRNA read counts after normalization
L3
PF PF+12h
Ecdysone titer
days
 DESeq
 Heatplus
Analyse d’expression différentielle
Larva
« Differential calling » avec le jeu complet de
données
Metamorphosis
Up-regulated
27
Down_regulated
27
« Differential calling » sans replicats
Metamorphosis
Up-regulated
0
Down_regulated
0
Message:
Le Deep Seq n’échappe pas au tests
statistiques
Les réplicats sont nécessaires pour estimer le
bruit biologique
PF
PF+12h
Naive and primed murine pluripotent stem cells have distinct miRNA signatures
1
miR- 302
Cluster
0.8
302/367
miR-290Cluster
290
295
Cluster
17-92
miR17-92
0.6
0.4
0.2
0
ES
EpiSC
M. Cohen-Tannoudji (Institut Pasteur)
A. Jouneau (INRA Jouy en Josas)
E. Heard (Institut Curie)
C. Antoniewski (Institut Pasteur)
ESC1 ESC2 EpiSC3EpiSC1EpiSC2
28/40
Normalized miR read count profiles
0.2
0.6
0.8
0.2
mmu-mir-302a
mmu-mir-291a
mmu-let-7a-1
1.0
1.0
0.8
0.8
% max count
% max count
mmu-let-7a-1
0.4
0.6
0.4
0.0
0.0
0.6
0.8
0.8
mmu-mir-302a
mmu-mir-17
0.4
0.2
0.4
0.6
0.6
0.2
0.2
0.4
0.2
0.4
0.6
0.8
% length
29/40
A lattice of miR read profiles for rapid, visual annotation
30/40
31/40
“Stereo” lattice reveals changes in miR biogenesis between ES and EpiSCs
ESC
EpiSC
% length
32/40
Small RNA signatures
1
−1
0
z−score
2
3
thorax
1
3
5
7
9
11
13
15
17
19
21
overlap (nt)
-100
AUGCUUUCAUGGCAUCCUUAC
AUGCUUUCAUGGCAUCCUUAC
AUGCUUUCAUGGCAUCCUUAC
AUGCUUUCAUGGCAUCCUUAC
AUGCUUUCAUGGCAUCCUUAC
AUGCUUUCAUGGCAUCCUUAC
AUGCUUUCAUGGCAUCCUUAC
AUGCUUUCAUGGCAUCCUUAC
UUUACGAAAGUACCGUAGGAA
|||||||||||||||||||||
123456789.........19...
UUUACGAAAGUACCGUAGGAA
UUUACGAAAGUACCGUAGGAA
UUUACGAAAGUACCGUAGGAA
UUUACGAAAGUACCGUAGGAA
+100
Signature piRNA
19−30nt RNAs
400
0
P-element
−200
−400
0
2000
4000
6000
8000
10000
12000
Coordinates (nt)
Cartographie des ARN de 24-26nt
d’ovaires de drosophile
2
z−score
4
signature
0
Normalized number of reads
200
0
2
4
6
8 10
13
16
19
22
25
28
overlap (nt)
-100
UUGCUUUCAUGGCAUCCUUACCGAUC
AGCUUCUUUACGAAACGAAAGUACCG
|||||||||||||||||||||
12345678910.............
+100
Signature piRNA
19−30nt RNAs
400
0
P-element
−200
−400
0
2000
4000
6000
8000
10000
12000
Coordinates (nt)
Cartographie des ARN de 24-26nt
d’ovaires de drosophile
2
z−score
4
signature
0
Normalized number of reads
200
0
2
4
6
8 10
13
16
19
22
25
28
overlap (nt)
-100
UUGCUUUCAUGGCAUCCUUACCGAUC
AGCUUCUUUACGAAACGAAAGUACCG
|||||||||||||||||||||
12345678910.............
+100