README.md 28.6 KB
Newer Older
Pierre LINDENBAUM's avatar
Pierre LINDENBAUM committed
1
# TP Next Generation Sequencing 2021
lindenb's avatar
1st  
lindenb committed
2

lindenb's avatar
cont    
lindenb committed
3
4
5
6
7
Au cours de cette séance de TP nous allons chercher les mutations de-novo d'un enfant dont l'exome ainsi que celui de ses deux parents ont été séquencés par la technologie Illumina. Il vous faudra :

  * aligner les reads sur un genome de reference
  * extraire les mutations
  * annoter ces variations.
lindenb's avatar
1st  
lindenb committed
8
9


Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
10
Certaines étapes seront effectuées plusieurs fois (père, mère, enfant) je vous conseille donc de noter vos commandes (ou de le mettre dans un shell script, ou dans un fichier *Makefile* pour ceux qui maitrisent ces outils)
lindenb's avatar
1st  
lindenb committed
11

lindenb's avatar
cont    
lindenb committed
12
13
Certains outils ont besoin de place; le cas échéant, faites de la place, videz vos répertoires temporaires.

Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
14
15


Pierre LINDENBAUM's avatar
Pierre LINDENBAUM committed
16
Le TP peut être long pour les débutants: concentrez vous sur la production du fichier VCF final, et cherchez à répondre aux questions de détails plus tard.
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
17

lindenb's avatar
cont    
lindenb committed
18

lindenb's avatar
cont    
lindenb committed
19
20
21
22
## workflow

![workflow.png](workflow.png "workflow")

lindenb's avatar
cont    
lindenb committed
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
## Petit rappel sur la structure d'un Makefile

> Les Makefiles sont des fichiers, généralement appelés **Makefile**, utilisés par le programme make pour exécuter un ensemble d'actions


Un Makefile est un fichier constitué de plusieurs règles de la forme : 


```
cible1 : dependance1 dependance2 dependance3 dependanceN
    commande1
    commande2

cible2 cible3 : dependance1 dependance2 dependance3 dependanceN
    commande1
    commande2

cible4:
    commande1
    commande2

```

Chaque commande est prefixé par une tabulation.

Lorsque l'ensemble des dépendances est analysé et si la cible ne correspond pas à un fichier existant ou si un fichier dépendance est plus récent que la régle, les différentes commandes sont exécutées. 

Raccourcis:

  * `$@` : nom de la cible courante
  * `$^` : toutes les dependance
  * `$<` : la premiere dependance

lindenb's avatar
1st  
lindenb committed
56
57
## Recupération des données FASTQ

lindenb's avatar
cont    
lindenb committed
58
59
60
Les donnees sont telechargeable a l'URL suivante:

https://uncloud.univ-nantes.fr/index.php/s/oEYQfdJQcReB7ez
lindenb's avatar
1st  
lindenb committed
61

Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
62
63
64
65
66
67
Téléchargez et extrayez le fichier zip avec 

```
unzip DATA.zip
```

Solena LE SCOUARNEC's avatar
Solena LE SCOUARNEC committed
68
Ne **JAMAIS** dezipper les fichiers .fq.gz (dans le futur, pour acceder au contenu de ces fichiers utilisez ` gunzip -c `, ou bien de nombreux outils NGS savent lire le format *gzip*).
lindenb's avatar
cont    
lindenb committed
69
70
71
72
73

Dans le repertoire *DATA* se trouvent les fichiers FASTQ contenant les short-reads pour father/mother/child.


```
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
74
75
76
77
78
79
DATA/child.R1.fq.gz
DATA/child.R2.fq.gz
DATA/father.R1.fq.gz
DATA/father.R2.fq.gz
DATA/mother.R1.fq.gz
DATA/mother.R2.fq.gz
lindenb's avatar
cont    
lindenb committed
80
81
82
```


lindenb's avatar
1st  
lindenb committed
83

Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
84
Nous travaillons en **paired-end**: il y a donc un fichier **R1** et **R2** pour chaque individu.
lindenb's avatar
1st  
lindenb committed
85

Pierre LINDENBAUM's avatar
Pierre LINDENBAUM committed
86
* (quand vous aurez votre fichier VCF) Comptez le nombre de lignes R1 chez l'enfant; Vérifiez que c'est un multiple de 4.
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
87
88

```
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
89
gunzip -c DATA/child.R1.fq.gz | wc -l
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
90
91
```

Pierre LINDENBAUM's avatar
Pierre LINDENBAUM committed
92
* (quand vous aurez votre fichier VCF) Comptez le nombre de lignes R2 chez l'enfant; Pourquoi ce même nombre ?
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
93
94

```
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
95
gunzip -c DATA/child.R2.fq.gz | wc -l
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
96
97
98
```


Pierre LINDENBAUM's avatar
Pierre LINDENBAUM committed
99
* (quand vous aurez votre fichier VCF) Comparez les 10 premiers noms de reads de l'enfant entre le fichier R1 et R2, quelle est la différence ?
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
100
101

```
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
102
103
$ gunzip -c DATA/child.R1.fq.gz | paste - - - - | cut -f 1 | head
$ gunzip -c DATA/child.R2.fq.gz | paste - - - - | cut -f 1 | head
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
104
105
```

Pierre LINDENBAUM's avatar
Pierre LINDENBAUM committed
106
*  (quand vous aurez votre fichier VCF) Quelle est la longueur des short-reads R1 chez l'enfant ?
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
107
108

```
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
109
$ gunzip -c DATA/child.R1.fq.gz | paste - - - - | cut -f 2 | head -n 1 | tr -d '\n' | wc -c
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
110
111
```

Pierre LINDENBAUM's avatar
Pierre LINDENBAUM committed
112
* (quand vous aurez votre fichier VCF) Dans les short-reads R1 de l'enfant, extrayez la premiere base de la séquence d'ADN, quelle est la proportion de ces bases (nombre de A , nombre de T, nombre de G, etc... ?). 
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
113
114

```
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
115
$ gunzip -c DATA/child.R1.fq.gz | paste - - - - | cut -f 2 | cut -c 1| sort | uniq -c | sort -n
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
116
117
```

Pierre LINDENBAUM's avatar
Pierre LINDENBAUM committed
118
* (quand vous aurez votre fichier VCF) Quelles sont les bases que l'on trouve dans les séquences d'ADN, n'y a-t-il que A-T-G-C ?
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
119
120

```
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
121
$ gunzip -c DATA/child.R1.fq.gz | paste - - - - | cut -f 2 | tr -d 'ATGC\n'
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
122
123
```

Pierre LINDENBAUM's avatar
Pierre LINDENBAUM committed
124
* (quand vous aurez votre fichier VCF) Quels sont les caractères ascii de qualité que l'on trouve dans les lignes de *qualité* ?
lindenb's avatar
1st  
lindenb committed
125

Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
126
```
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
127
$ gunzip -c DATA/child.R1.fq.gz | paste - - - - | cut -f 4 | grep -o . | sort | uniq
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
128
129
```

lindenb's avatar
1st  
lindenb committed
130
131
## Recupération du génome de reference

Solena LE SCOUARNEC's avatar
Solena LE SCOUARNEC committed
132
Pour aller vite, nous nous contenterons des chromosomes 22 et de l'ADN mitochondrial humains. 
lindenb's avatar
1st  
lindenb committed
133

lindenb's avatar
cont    
lindenb committed
134

Pierre LINDENBAUM's avatar
Pierre LINDENBAUM committed
135
Telechargez les séquences fa.gz  **chr22** et **chrM** de la version **hg19/GRCh37** du genome humain.
lindenb's avatar
1st  
lindenb committed
136

Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
137
```
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
138
139
$ wget -O chr22.fa.gz "http://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/chr22.fa.gz"
$ wget -O chrM.fa.gz "http://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/chrM.fa.gz"
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
140
```
Pierre LINDENBAUM's avatar
Pierre LINDENBAUM committed
141
142


Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
143
dézippez les deux fichiers.
Pierre LINDENBAUM's avatar
Pierre LINDENBAUM committed
144
145

```
Pierre LINDENBAUM's avatar
Pierre LINDENBAUM committed
146
$ gunzip chr22.fa.gz chrM.fa.gz
Pierre LINDENBAUM's avatar
Pierre LINDENBAUM committed
147

Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
148
```
Pierre LINDENBAUM's avatar
Pierre LINDENBAUM committed
149
150


lindenb's avatar
cont    
lindenb committed
151
152
153



Pierre LINDENBAUM's avatar
Pierre LINDENBAUM committed
154
* (quand vous aurez votre fichier VCF) Quelle est la taille du chromosome 22 ? Quelle est la taille du chrM ?
Pierre LINDENBAUM's avatar
Pierre LINDENBAUM committed
155

Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
156
157
158
159
```
$ grep -v '>' chr22.fa | tr -d '\n' | wc -c
$ grep -v '>' chrM.fa | tr -d '\n' | wc -c
```
Pierre LINDENBAUM's avatar
Pierre LINDENBAUM committed
160

Pierre LINDENBAUM's avatar
Pierre LINDENBAUM committed
161
* (quand vous aurez votre fichier VCF) Quelles sont les premières bases du chromosome 22 ? 
Pierre LINDENBAUM's avatar
Pierre LINDENBAUM committed
162

Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
163
164
165
```
$ more chr22.fa
```
Pierre LINDENBAUM's avatar
Pierre LINDENBAUM committed
166

Pierre LINDENBAUM's avatar
Pierre LINDENBAUM committed
167
(quand vous aurez votre fichier VCF) Pourquoi observez vous cela ? Combien de lignes y a-t-il avant que cela change ?
Pierre LINDENBAUM's avatar
Pierre LINDENBAUM committed
168

Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
169
170
171
```
$ cat -n chr22.fa | grep -v '>'  | grep -v N -m1
```
Pierre LINDENBAUM's avatar
Pierre LINDENBAUM committed
172

Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
173
A l'aide de `cat`, concatenez ces deux fichiers dans un fichier que vous nommerez **ref.fa**.
Pierre LINDENBAUM's avatar
Pierre LINDENBAUM committed
174
175

```
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
176
$ cat chr22.fa  chrM.fa > ref.fa
Pierre LINDENBAUM's avatar
Pierre LINDENBAUM committed
177
178
```

lindenb's avatar
cont    
lindenb committed
179
180
181
182
183
184
185
186
187
188
189
190
Par exemple dans un Makefile ça donnera:

```
ref.fa: chr22.fa chrM.fa
    cat $^ > $@

chr22.fa chrM.fa:
    wget -O $@.fa.gz "http://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/$@.fa.gz"
    gunzip $@.fa.gz

```

Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
191

lindenb's avatar
cont    
lindenb committed
192
##  Mapper les Fastq sur le genome de reference
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
193

Solena LE SCOUARNEC's avatar
Solena LE SCOUARNEC committed
194
A priori tous les outils ont été installés. Pas besoin d'installer samtools, bwa etc... comme décrit ci-dessous. En revanche vous devez indiquer au `bash` où se trouvent les outils en faisant:
Pierre LINDENBAUM's avatar
Pierre LINDENBAUM committed
195
196

```
lindenb's avatar
cont    
lindenb committed
197
export PATH="/usr/local/opt/sam_bwa_20200917/bin:${PATH}"
Pierre LINDENBAUM's avatar
Pierre LINDENBAUM committed
198
```
Pierre LINDENBAUM's avatar
Pierre LINDENBAUM committed
199

Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
200
201
202
203
204
on vérifie en tapant:

```
$ which samtools bwa

lindenb's avatar
cont    
lindenb committed
205
206
/usr/local/opt/sam_bwa_20200917/bin/samtools
/usr/local/opt/sam_bwa_20200917bin/bwa
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
207
208
```

lindenb's avatar
cont    
lindenb committed
209
Refaite cette operation a chaque fois 
Pierre LINDENBAUM's avatar
Pierre LINDENBAUM committed
210

lindenb's avatar
1st  
lindenb committed
211
212
## Prise en main de BWA

lindenb's avatar
cont    
lindenb committed
213
*BWA* est l'outil que va nous permettre de mapper les reads sur le génome de reference. *BWA* lui-même contient une série de sous programmes que l'on peut lister en tapant `bwa`
lindenb's avatar
1st  
lindenb committed
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240

```
$ bwa

Program: bwa (alignment via Burrows-Wheeler transformation)

Usage:   bwa <command> [options]

Command: index         index sequences in the FASTA format
         mem           BWA-MEM algorithm
(...)
```

On peut obtenir l'aide de chaque sous-programme en invoquant le nom du sous-programme.


```
$ bwa index

Usage:   bwa index [-a bwtsw|is] [-c] <in.fasta>

Options: -a STR    BWT construction algorithm: bwtsw or is [auto]
         -p STR    prefix of the index [same as fasta name]
         -6        index files named as <in.fasta>.64.* instead of <in.fasta>.* 

```

Pierre LINDENBAUM's avatar
Pierre LINDENBAUM committed
241

lindenb's avatar
1st  
lindenb committed
242
243
## Prise en main de SAMTOOLS

lindenb's avatar
cont    
lindenb committed
244
samtools est le couteau suisse des formats en NGS. Tout comme bwa, samtools lui-même contient une série de sous programmes que l'on peut afficher en tapant `samtools`
lindenb's avatar
1st  
lindenb committed
245

Pierre LINDENBAUM's avatar
Pierre LINDENBAUM committed
246

lindenb's avatar
1st  
lindenb committed
247
248
249
250
```
$ samtools

Program: samtools (Tools for alignments in the SAM format)
lindenb's avatar
cont    
lindenb committed
251
Version: 1.9-250-ga6a160b (using htslib 1.10.2-22-gbfc9f0d)
lindenb's avatar
1st  
lindenb committed
252
253
254

Usage:   samtools <command> [options]

lindenb's avatar
cont    
lindenb committed
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
Commands:
  -- Indexing
     dict           create a sequence dictionary file
     faidx          index/extract FASTA
     fqidx          index/extract FASTQ
     index          index alignment

  -- Editing
     calmd          recalculate MD/NM tags and '=' bases
     fixmate        fix mate information
     reheader       replace BAM header
     targetcut      cut fosmid regions (for fosmid pool only)
     addreplacerg   adds or replaces RG tags
     markdup        mark duplicates

  -- File operations
     collate        shuffle and group alignments by name
     cat            concatenate BAMs
     merge          merge sorted alignments
     mpileup        multi-way pileup
     sort           sort alignment file
     split          splits a file by read group
     quickcheck     quickly check if SAM/BAM/CRAM file appears intact
     fastq          converts a BAM to a FASTQ
     fasta          converts a BAM to a FASTA

  -- Statistics
     bedcov         read depth per BED region
     coverage       alignment depth and percent coverage
     depth          compute the depth
     flagstat       simple stats
     idxstats       BAM index stats
     phase          phase heterozygotes
     stats          generate stats (former bamcheck)

  -- Viewing
     flags          explain BAM flags
     tview          text alignment viewer
     view           SAM<->BAM<->CRAM conversion
     depad          convert padded BAM to unpadded BAM
lindenb's avatar
1st  
lindenb committed
295
296
297
298
299
300
301
302
303
304

(...)
```

On peut également obtenir l'aide de chaque sous-programme en invoquant le nom du sous-programme.


```
$ samtools view

lindenb's avatar
cont    
lindenb committed
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
Usage: samtools view [options] <in.bam>|<in.sam>|<in.cram> [region ...]

Options:
  -b       output BAM
  -C       output CRAM (requires -T)
  -1       use fast BAM compression (implies -b)
  -u       uncompressed BAM output (implies -b)
  -h       include header in SAM output
  -H       print SAM header only (no alignments)
  -c       print only the count of matching records
  -o FILE  output file name [stdout]
  -U FILE  output reads not selected by filters to FILE [null]
  -t FILE  FILE listing reference names and lengths (see long help) [null]
  -X       include customized index file
  -L FILE  only include reads overlapping this BED FILE [null]
  -r STR   only include reads in read group STR [null]
  -R FILE  only include reads with read group listed in FILE [null]
  -d STR:STR
           only include reads with tag STR and associated value STR [null]
  -D STR:FILE
           only include reads with tag STR and associated values listed in
           FILE [null]
  -q INT   only include reads with mapping quality >= INT [0]
  -l STR   only include reads in library STR [null]
  -m INT   only include reads with number of CIGAR operations consuming
           query sequence >= INT [0]
  -f INT   only include reads with all  of the FLAGs in INT present [0]
  -F INT   only include reads with none of the FLAGS in INT present [0]
  -G INT   only EXCLUDE reads with all  of the FLAGs in INT present [0]
  -s FLOAT subsample reads (given INT.FRAC option value, 0.FRAC is the
           fraction of templates/read pairs to keep; INT part sets seed)
  -M       use the multi-region iterator (increases the speed, removes
           duplicates and outputs the reads as they are ordered in the file)
  -x STR   read tag to strip (repeatable) [null]
  -B       collapse the backward CIGAR operation
  -?       print long help, including note about region specification
  -S       ignored (input format is auto-detected)
  --no-PG  do not add a PG line
      --input-fmt-option OPT[=VAL]
               Specify a single input file format option in the form
               of OPTION or OPTION=VALUE
  -O, --output-fmt FORMAT[,OPT[=VAL]]...
               Specify output format (SAM, BAM, CRAM)
      --output-fmt-option OPT[=VAL]
               Specify a single output file format option in the form
               of OPTION or OPTION=VALUE
  -T, --reference FILE
               Reference sequence FASTA FILE [null]
  -@, --threads INT
               Number of additional threads to use [0]
      --write-index
               Automatically index the output files [off]
      --verbosity INT
               Set level of verbosity

lindenb's avatar
1st  
lindenb committed
360
361
362
363
364
365

```


## Construction de l'index du genome de reference

Pierre LINDENBAUM's avatar
Pierre LINDENBAUM committed
366
On utilise `bwa index` pour construire l'index de la reference.
lindenb's avatar
cont    
lindenb committed
367

lindenb's avatar
1st  
lindenb committed
368
```
lindenb's avatar
cont    
lindenb committed
369
370
371
372
373
374
375
376
$ bwa index ref.fa

[bwa_index] Pack FASTA... 1.02 sec
[bwa_index] Construct BWT for the packed sequence...
[BWTIncCreate] textLength=101670074, availableWord=19153804
[BWTIncConstructFromPacked] 10 iterations done. 31594474 characters processed.
[BWTIncConstructFromPacked] 20 iterations done. 58366762 characters processed.
(...)
lindenb's avatar
1st  
lindenb committed
377
378
```

Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
379
380
381
382
383
vérifiez que des fichiers d'index ont été créés

```
ls ref.fa.*
```
lindenb's avatar
1st  
lindenb committed
384
385
386
387
388
389
390
391
392


## Mapping des séquences fastq sur le génome de reference

maintenant que l'index du génome de référence a été créé, nous cherchons a connaitre la position des reads sur ce génome de reference.

le synopsis de la commande bwa pour mapper les reads est :

```
Pierre LINDENBAUM's avatar
cont    
Pierre LINDENBAUM committed
393
bwa mem  -R '@RG\tID:SAMPLENAME-aa\tSM:SAMPLENAME' ref.fa sample.R1.gz sample.R2.gz > sample.sam
lindenb's avatar
1st  
lindenb committed
394
395
```

Solena LE SCOUARNEC's avatar
Solena LE SCOUARNEC committed
396
l'option `-R` permet de créer un **READ-GROUP** pour associer les reads à un échantillon (donc ici, remplacez SAMPLENAME par 'child').
lindenb's avatar
1st  
lindenb committed
397

Pierre LINDENBAUM's avatar
PG    
Pierre LINDENBAUM committed
398

lindenb's avatar
1st  
lindenb committed
399
400
la sortie standard de cette commande est un fichier 'texte' **SAM**.

Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
401
Mappez les reads de l'enfant `child.R1.fq.gz` et `child.R2.fq.gz` avec cette commande et redirigez la sortie standard vers un fichier `child_unsorted.sam`
lindenb's avatar
1st  
lindenb committed
402

Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
403
```
Pierre LINDENBAUM's avatar
aa    
Pierre LINDENBAUM committed
404
bwa mem  -R '@RG\tID:child\tSM:child' ref.fa DATA/child.R1.fq.gz DATA/child.R2.fq.gz > child_unsorted.sam
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
405
406
```

Pierre LINDENBAUM's avatar
Pierre LINDENBAUM committed
407
* (quand vous aurez votre fichier VCF)  Que renvoie la commande suivante ?
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
408
409

```
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
410
file child_unsorted.sam
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
411
412
```

Pierre LINDENBAUM's avatar
Pierre LINDENBAUM committed
413
* (quand vous aurez votre fichier VCF) Combien y-a-t-il de lignes dans ce fichier sam ?
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
414
415

```
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
416
wc -l child_unsorted.sam 
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
417
418
419
```


Pierre LINDENBAUM's avatar
Pierre LINDENBAUM committed
420
* (quand vous aurez votre fichier VCF) Combien y-a-t-il de lignes dans l'en-tête ?
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
421
422

```
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
423
grep  '^@' child_unsorted.sam | wc -l
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
424
425
426
427
```



Pierre LINDENBAUM's avatar
Pierre LINDENBAUM committed
428
* (quand vous aurez votre fichier VCF) Observez les lignes de l'en-tête commençant par `@SQ` , y a-t-il une différence avec le nombre de bases que vous aviez observé dans les séquences Fasta des chr22 et chrM ?
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
429

Pierre LINDENBAUM's avatar
Pierre LINDENBAUM committed
430
* (quand vous aurez votre fichier VCF) Retrouvez la ligne `@RG` du 'read-group' dans l'en-tête, vérifiez que tous les reads dans ce fichier portent une reférence à ce read groupe dans la colonne de méta-données.
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
431

Pierre LINDENBAUM's avatar
Pierre LINDENBAUM committed
432
* (quand vous aurez votre fichier VCF)  Combien y-a-t-il d'alignements ? Comparez au nombre de reads dans les fichiers child.R1.fq.gz et child.R2.fq.gz
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
433

Pierre LINDENBAUM's avatar
Pierre LINDENBAUM committed
434
* (quand vous aurez votre fichier VCF) Identifiez les colonnes d'un fichier SAM:
lindenb's avatar
1st  
lindenb committed
435
436
437
438
439
440
441
442
443
444
445
  * READ NAME
  * FLAG
  * CHROMOSOME
  * POS
  * MAPPING QUALITY
  * CIGAR (réprésentation compacte de l'alignement)
  * MATE CHROMOSOME
  * MATE POS
  * DISTANCE READ-MATE (négatif si mate placé devant)
  * DNA SEQ
  * QUALITIES
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
446
447

```
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
448
grep -v '^@' child_unsorted.sam | head -n 1 | tr "\t" "\n" | cat -n
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
449
450
```

Pierre LINDENBAUM's avatar
Pierre LINDENBAUM committed
451
* (quand vous aurez votre fichier VCF) La deuxième colonne contient le SAM FLAG: des méta-informations sur le mapping des reads. Quel est le SAM Flag trouvé le plus souvent ? Cherchez sa signification dans http://broadinstitute.github.io/picard/explain-flags.html ?
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
452

Pierre LINDENBAUM's avatar
Pierre LINDENBAUM committed
453
* (quand vous aurez votre fichier VCF) Prenez un read au hasard, verifiez que vous retrouvez la même séquence+qualité dans le fichier FASTQ originel (attention la séquence peut-être reverse-complémentée)
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
454

Pierre LINDENBAUM's avatar
Pierre LINDENBAUM committed
455
* (quand vous aurez votre fichier VCF) Y'a t-il des CIGAR String représentant des insertions/deletions ou contenant du clipping ? C'est à dire que la cigar-string est différente de `(nombre)M`
lindenb's avatar
1st  
lindenb committed
456

lindenb's avatar
cont    
lindenb committed
457
## Conversion de SAM en BAM
lindenb's avatar
1st  
lindenb committed
458

Solena LE SCOUARNEC's avatar
Solena LE SCOUARNEC committed
459
pour permettre aux algorithmes d'aller plus vite, nous utilisons samtools pour convertir le fichier en format texte/SAM `child_unsorted.sam` en format binaire/BAM `child_unsorted.bam`
lindenb's avatar
1st  
lindenb committed
460
461


lindenb's avatar
cont    
lindenb committed
462
### Synopsis:
lindenb's avatar
1st  
lindenb committed
463

lindenb's avatar
cont    
lindenb committed
464
```
lindenb's avatar
cont    
lindenb committed
465
samtools view -O BAM -o out.bam in.sam
lindenb's avatar
cont    
lindenb committed
466
467
```

lindenb's avatar
cont    
lindenb committed
468
469

* Option **-O BAM** : en souhaite écrite du binaire (**BAM**)
lindenb's avatar
cont    
lindenb committed
470
471
* Option **-o** : nom du fichier de sortie

Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
472
Generez le fichier `child_unsorted.bam`
lindenb's avatar
cont    
lindenb committed
473

Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
474
```
lindenb's avatar
cont    
lindenb committed
475
samtools view -O BAM -o child_unsorted.bam child_unsorted.sam
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
476
477
```

Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
478
* que renvoie la commande `file child_unsorted.bam` ?
lindenb's avatar
cont    
lindenb committed
479
480
481
482
483
484
485
486
487
* affichez les options de `samtools view` avec

```
samtools view
```

* affichez le corps du BAM

```
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
488
samtools view  child_unsorted.bam
lindenb's avatar
cont    
lindenb committed
489
490
491
492
493
```

* affichez l'en-tete du BAM avec

```
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
494
samtools view -H child_unsorted.bam
lindenb's avatar
cont    
lindenb committed
495
496
497
498
499
```

* affichez le bam complet avec

```
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
500
samtools view -h child_unsorted.bam
lindenb's avatar
cont    
lindenb committed
501
502
```

Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
503
* en vous aidant de  http://broadinstitute.github.io/picard/explain-flags.html et de l'option `-f (flag)` de samtools view, comptez les reads qui étaient à l'origine dans le fichier child.R1.fq.gz (a.k.a 'first in pair')
lindenb's avatar
cont    
lindenb committed
504
505


Solena LE SCOUARNEC's avatar
Solena LE SCOUARNEC committed
506
## Tri du BAM sur la position génomique
lindenb's avatar
cont    
lindenb committed
507

lindenb's avatar
cont    
lindenb committed
508
Pour l'instant les reads ne sont pas triés sur la position génomique (chrom/start). Pour aller plus vite, les algorithmes détectant les mutations ont besoin que les reads soient triés. Pour cela on utilise `samtools sort`:
lindenb's avatar
cont    
lindenb committed
509
510
511
512

### Synopsis:

```
lindenb's avatar
cont    
lindenb committed
513
514
515
Usage: samtools sort [options...] [in.bam]
Options:
  -l INT     Set compression level, from 0 (uncompressed) to 9 (best)
lindenb's avatar
cont    
lindenb committed
516
  -u         Output uncompressed data (equivalent to -l 0)
lindenb's avatar
cont    
lindenb committed
517
  -m INT     Set maximum memory per thread; suffix K/M/G recognized [768M]
lindenb's avatar
cont    
lindenb committed
518
519
520
  -M         Use minimiser for clustering unaligned/unplaced reads
  -K INT     Kmer size to use for minimiser [20]
  -n         Sort by read name (not compatible with samtools index command)
lindenb's avatar
cont    
lindenb committed
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
  -t TAG     Sort by value of TAG. Uses position as secondary index (or read name if -n is set)
  -o FILE    Write final output to FILE rather than standard output
  -T PREFIX  Write temporary files to PREFIX.nnnn.bam
  --no-PG    do not add a PG line
      --input-fmt-option OPT[=VAL]
               Specify a single input file format option in the form
               of OPTION or OPTION=VALUE
  -O, --output-fmt FORMAT[,OPT[=VAL]]...
               Specify output format (SAM, BAM, CRAM)
      --output-fmt-option OPT[=VAL]
               Specify a single output file format option in the form
               of OPTION or OPTION=VALUE
      --reference FILE
               Reference sequence FASTA FILE [null]
  -@, --threads INT
               Number of additional threads to use [0]
      --verbosity INT
               Set level of verbosity
lindenb's avatar
cont    
lindenb committed
539

lindenb's avatar
cont    
lindenb committed
540
541
```

lindenb's avatar
cont    
lindenb committed
542
543
544



Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
545
Triez `child_unsorted.bam` pour produire un fichier `child.bam`.
lindenb's avatar
cont    
lindenb committed
546

Pierre LINDENBAUM's avatar
cont    
Pierre LINDENBAUM committed
547
548
549
Ici l'option  `-T tmp` definit un prefix pour trier sur disque si les donnees ne tiennent pas en memoire: https://en.wikipedia.org/wiki/External_sorting


Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
550
```
lindenb's avatar
cont    
lindenb committed
551
samtools sort -T tmp -O BAM -o child.bam child_unsorted.bam
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
552
553
```

lindenb's avatar
cont    
lindenb committed
554
* En utilisant 'samtools view', vérifiez que le fichier est trié en regardant la première ligne de son en-tête qui doit contenir le mot `coordinate`.
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
555
556

```
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
557
samtools view -H child.bam 
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
558
559
```

lindenb's avatar
cont    
lindenb committed
560
561
* En utilisant 'samtools view', Verifiez que le fichier est trié en observant la position croissante CHROM/POS des reads dans le fichier BAM.

Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
562
```
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
563
samtools view  child.bam | cut -f 3,4
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
564
565
```

lindenb's avatar
cont    
lindenb committed
566
567
568
569
570
571
572
## Petite digression: Le format CRAM.

Prenons le temps de voir le nouveau format CRAM qui est l'équivalent de BAM mais nous ne stockons que les différences avec le genome de reference.

Pour génerer du *CRAM* à partir de BAM nous devons indiquer où se trouve le genome de reference.

```
Pierre LINDENBAUM's avatar
cont    
Pierre LINDENBAUM committed
573
samtools view -O CRAM --reference ref.fa -o child.cram child.bam
lindenb's avatar
cont    
lindenb committed
574
575
576
```

* Comparer la taille des fichiers **CRAM** et **BAM**. Le fichier **CRAM** est-il bien plus petit que le **BAM* ?.
Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
577

lindenb's avatar
cont    
lindenb committed
578
* Pour voir les reads dans un **CRAM** il faut utiliser `samtools view` en joignant le genome de reference.
lindenb's avatar
cont    
lindenb committed
579

lindenb's avatar
cont    
lindenb committed
580
581
582
583
```
samtools view --reference reference.fa child.cram
```

Pierre LINDENBAUM's avatar
Pierre LINDENBAUM committed
584
* (quand vous aurez votre fichier VCF) Y a-t-il le même nombre de reads entre le fichier  `child.cram` et `child.bam` ?
lindenb's avatar
cont    
lindenb committed
585
586
587

## Indexation du BAM

lindenb's avatar
cont    
lindenb committed
588
589
Nous revenons au fichier **BAM**.

Solena LE SCOUARNEC's avatar
Solena LE SCOUARNEC committed
590
L'indexation consiste à créer un fichier d'index **'.bai'** à partir d'un bam trié. Cet index permet d'accéder rapidement aux reads mappés dans une région génomique précise. Ici on se sert de la commande `samtools index`
lindenb's avatar
cont    
lindenb committed
591
592
593
594
595
596
597
598


### Synopsis:

```
samtools index  in.bam
```

Pierre LINDENBAUM's avatar
cont    
Pierre LINDENBAUM committed
599

lindenb's avatar
cont    
lindenb committed
600
Indexez le fichier **child.bam** , vérifiez le fonctionnement de l'index en effectuant des requêtes par région:
lindenb's avatar
cont    
lindenb committed
601
602
603
604

```
samtools view  child.bam chr22
samtools view  child.bam chrM
Pierre LINDENBAUM's avatar
cont    
Pierre LINDENBAUM committed
605
# la region ci dessous ne donnera rien si il n'y a pas de read dans la region
lindenb's avatar
cont    
lindenb committed
606
607
608
609
610
611
612
613
614
615
616
samtools view  child.bam chr22:23707000-23803000
```

## Indexation du genome de reference avec samtools.

Samtools doit pouvoir accèder rapidement à des régions du génome de reference. Pour cela il lui faut préalablement créer un petit fichier d'index pour la séquence de reference:


### Synopsis:

```
Pierre LINDENBAUM's avatar
cont    
Pierre LINDENBAUM committed
617
samtools faidx ref.fa
lindenb's avatar
cont    
lindenb committed
618
619
```
 
lindenb's avatar
cont    
lindenb committed
620
Indexez votre génome de reference. La commande doit generer un fichier **.fai** qui contient les informations nécessaires (offset de la séquence dans le fichier, taille des lignes,.... ) pour accèder à des sous-séquences.
lindenb's avatar
cont    
lindenb committed
621

lindenb's avatar
cont    
lindenb committed
622
Une fois indexé, on peut rapidement obtenir des sous-régions du genome de ref au format fasta.
lindenb's avatar
cont    
lindenb committed
623
624
625
626

### Synopsis:

```
Pierre LINDENBAUM's avatar
cont    
Pierre LINDENBAUM committed
627
samtools faidx ref.fa chr:start-end
lindenb's avatar
cont    
lindenb committed
628
629
630
631
632
633
634
```
 
En utilisant `samtools faidx` extrayez la sequence `chrM:50-200`
 

## Visualisation interactive avec tview 

Solena LE SCOUARNEC's avatar
Solena LE SCOUARNEC committed
635
La commande `samtools tview` permet de visualiser interactivement les reads dans une région donnée. Il faut que le BAM et la séquence de réference soient indexés.
lindenb's avatar
cont    
lindenb committed
636
637
638
639
640


### Synopsis:

```
Pierre LINDENBAUM's avatar
cont    
Pierre LINDENBAUM committed
641
samtools tview  sorted.bam ref.fa
lindenb's avatar
cont    
lindenb committed
642
643
```

Pierre LINDENBAUM's avatar
Pierre LINDENBAUM committed
644

Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
645
646
Ouvrez `child.bam` avec `samtools tview child.bam ref.fa`.

Pierre LINDENBAUM's avatar
Pierre LINDENBAUM committed
647
Appuyez sur la touch 'g' ('goto') et tapez la position `chr:pos` (par exemple `chr22:22353001`) du read précédent. Utilisez les flêches pour naviguer, '?' pour obtenir de l'aide, 'q' pour quitter.
lindenb's avatar
1st  
lindenb committed
648

Pierre LINDENBAUM's avatar
tabix    
Pierre LINDENBAUM committed
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
```
        22352991  22353001  22353011  22353021  22353031  22353041
CAGGTAGGAAAAGGCCTCAGAGACCAGGGTCAGCCACACAGCCTGATTCTGACTCTCGTGTCAAAGATCA
.................M.S.........Y....MM.W..................TS............
CAG TAGGAAAAGGCCTCAGAGACCAGGGTCAGCCACACAGCCTGATTCTGACTCTTGTGTCAAAGATCA
CAGGTAG  aaaggcctcagagaccagggtcagccacacagcctgattctgactcttgtgtcaaagatca
CAGGTAGGAAA GGCCTCAGAGACCAGGGTCAGCCACACAGCCTGATTCTGACTCTTGTGTCAAAGATCA
CAGGTAGGAAA  GCCTCAGAGACCAGGGTCAGCCACACAGCCTGATTCTGACTCTTGTGTCAAAGATCA
CAGGTAGGAAAA  CCTCAGAGACCAGGGTCAGCCACACAGCCTGATTCTGACTCTTGTGTCAAAGATCA
caggtaggaaaagg CTCAGAGACCAGGGTCAGCCACACAGCCTGATTCTGACTCTTGTGTCAAAGATCA
CAGGTAGGAAAAGGCCTCA  GACCAGGGTCAGCCACACAGCCTGATTCTGACTCTTGTGTCAAAGATCA
CAGGTAGGAAAAGGCCTCAGAG   AGGGTCAGCCACACAGCCTGATTCTGACTCTTGTGTCAAAGATCA
CAGGTAGGAAAAGGCCTAACAGAC AGGGTCAGCCACACAGCCTGATTCTGACTCTTGTGTCAAAGATCA
CAGGTAGGAAAAGGCCTCAGAGACCA    CAGCCACACAGCCTGATTCTGACTCTTGTGTCAAAGATCA
CAGGTAGGAAAAGGCCTCAGAGACCA    CAGCCACACAGCCTGATTCTGACTCTTGTGTCAAAGATCA
CAGGTAGGAAAAGGCCTCAGAGACCAGGGT agccacacagcctgattctgactcttgtgtaaaagatca
CAGGTAGGAAAAGGCCTAACAGACCAGGGC  GCCACACAGCCTGATTCTGACTCTTGTGTCAAAGATCA
CAGGTAGGAAAAGGCCTAACAGACCAGGGC  gccacacagcctgattctgactcttgtgtcaaagatca
CAGGTAGGAAAAGGCCTCAGAGACCAGGGTCA ccacacagcctgattctgactcttgtgtcaaagatca
CAGGTAGGAAAAGGCCTCAGAGACCAGGGTCAGCCA   AGCCTGATTCTGACTCTTGTGTCAAAGATCA
CAGGTAGGAAAAGGCCTCAGAGACCAGGGTCAGCCACACAG  TGATTCTGACTCTTGTGTCAAAGATCA
CAGGTAGGAAAAGGCCTCAGAGACCAGGGTCAGCCACACAGCC GATTCTGACTCTTGTGTCAAAGATCA
CAGGTAGGAAAAGGCCTCAGAGACCAGGGTCAGCCACACAGCCT ATTCTGACTCTTGTGTCAAAGATCA
C*GGTAGGAAAAGGCCTCAGAGACCAGGGTCAGCCACACAGCCTG TTCTGACTCTTGTGTCAAAGATCA
CAGGTAGGAAAAGGCCTCAGAGACCAGGGTCAGCCACACAGCCTGA TCTGACTCTTGTGTCAAAGATCA
CAGGTAGGAAAAGGCCTCAGAGACCAGGGTCAGCCACACAGCCTGATTCTGACT  TGTGTCAAAGATCA
CAGGTAGGAAAAGGCCTCAGAGACCAGGGTCAGCCACACAGCCTGATTCTGACTCT GTGTCAAAGATCA
CAGGTAGGAAAAGGCCTCAGAGACCAGGGTCAGCCACACAGCCTGATTCTGACTCTT TGTCAAAGATCA
CAGGTAGGAAAAGGCCTCAGAGACCAGGGTCAGCCACACAGCCTGATTCTGACTCTTGTG CAAAGATCA
CAGGTAGGAAAAGGCCTCAGAGACCAGGGTCAGCCACACAGCCTGATTCTGACTCTTGTGTCAAA ATCA
CAGGTAGGAAAAGGCCTCAGAGACCAGGGTCAGCCACACAGCCTGATTCTGACTCTTGTGTCAAAG TCA
CAGGTAGGAAAAGGCCTCAGAGGCCAGGGTCAGCCACACAGCCTGATTCTGACTCTTGTGTCAAAGATCA
CAGGTAGGAAAAGGCCTCAGAGACCAGGGTCAGCCACACAGCCTGATTCTGACTCTTGTGTCAAAGATCA
CAGGTAGGAAAAGGCCTAACAGACCAGGGCCAGCACCTCAGCCTGATTCTGACTCTTCTGGCAAAGATCC
caggtaggaaaaggcctaacagaccagggccagcacctcagcctgattctgactcttctggcaaagatcc
caggtaggaaaaggcctcagagaccagggtcagccacacagcctgattctgactcttgtgtcaaagatca
CAGGTAGGAAAAGGCCTCAGAGACCAGGGTCAGCCACACAGCCTGATTCTGACTCTTGTGTCAAAGATCA
caggtaggaaaaggcctaacagaccagggccagcacctcagcctgattctgactcttctggcaaagatcc
CAGGTAGGAAAAGGCCTCAGAGACCAGGGTCAGCCACACAGCCTGATTCTGACTCTTGTGTCAAAGATCA
CAGGTAGGAAAAGGCCTCAGAGACCAGGGTCAGCCACACAGCCTGATTCTGACTCTTGTGTCAAAGATCA
CAGGTAGGAAAAGGCCTCAGAGACCAGGGTCAGCCACACAGCCTGATTCTGACTCTTGTGTCAAAGATCA
CAGGTAGGAAAAGGCCTCAGAGACCAGGGTCAGCCACACAGCCTGATTCTGACTCTTGTGTCAAAGATCA
CAGGTAGGAAAAGGCCTAACAGACCAGGGCCAGCACCTCAGCCTGATTCTGACTCTTCTGGCAAAGATCC
```
lindenb's avatar
1st  
lindenb committed
693
694


Solena LE SCOUARNEC's avatar
Solena LE SCOUARNEC committed
695
## Mapping des fichiers parentaux.
lindenb's avatar
cont    
lindenb committed
696
697
698
699
700
701


### Mother 

Utilisez

Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
702
703
  * DATA/mother.R1.fq.gz
  * DATA/mother.R2.fq.gz
lindenb's avatar
cont    
lindenb committed
704
705
706
707
708
709
710

et créez le fichier **mother.bam**

### Father 

Utilisez

Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
711
712
  * DATA/father.R1.fq.gz
  * DATA/father.R2.fq.gz
lindenb's avatar
cont    
lindenb committed
713
714
715
716
717
718
719
720

et créez le fichier **father.bam**


## Calling des mutations.

Le calling des mutation se fait à l'aide de 

lindenb's avatar
cont    
lindenb committed
721
* `bcftools mpileup` qui va générer les bases trouvées à toutes les positions pour les 3 bams. Cet outil génére un fichier binaire de variants (BCF).
lindenb's avatar
cont    
lindenb committed
722
* `bcftools view` va transformer le **BCF** (binaire) en format *VCF* (*Variant Call Format*).
lindenb's avatar
cont    
lindenb committed
723
724
725
726

### Synopsis:

```
lindenb's avatar
cont    
lindenb committed
727

Pierre LINDENBAUM's avatar
cont    
Pierre LINDENBAUM committed
728
$ bcftools mpileup --output-type b --fasta-ref ref.fa --output mutations.bcf  -a 'FORMAT/AD'  -a 'FORMAT/DP'   file1.bam file2.bam ... fileN.bam 
lindenb's avatar
cont    
lindenb committed
729
730
```

lindenb's avatar
cont    
lindenb committed
731
732
* `-a FORMAT/DP` ajouter l'information du nombre de reads sous chaque mutation.
* `-a FORMAT/AD` ajouter l'information du nombre de reads REF/ ALT sous chaque mutation.
lindenb's avatar
cont    
lindenb committed
733
734
735
736

puis:

```
lindenb's avatar
cont    
lindenb committed
737
$ bcftools call  --ploidy GRCh37 -f GQ  --multiallelic-caller --variants-only --output-type z -o result.vcf.gz --format-fields GQ,GP mutations.bcf
lindenb's avatar
cont    
lindenb committed
738
739
740
741
742
743
```



## le fichier VCF

Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
744
Réperez la ligne commençant par "#CHROM" , c'est l'en-tête standard d'un fichier *VCF*
lindenb's avatar
cont    
lindenb committed
745
746
747
748
749
750
751
752
753
754
755
756
757


* 1. #CHROM
* 2.  POS
* 3.  ID
* 4.  REF
* 5.  ALT
* 6.  QUAL
* 7.  FILTER
* 8.  INFO
* 9.  FORMAT
* child, mother, father

Solena LE SCOUARNEC's avatar
Solena LE SCOUARNEC committed
758
Le contenu d'un fichier VCF sera étudié en détail lors de la prochaine séance.
lindenb's avatar
cont    
lindenb committed
759
760
761

## Annotez le fichier VCF

Pierre LINDENBAUM's avatar
gff    
Pierre LINDENBAUM committed
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
Nous allons utiliser bcftools csq pour annoter le VCF:

### Synopsis
```
About: Haplotype-aware consequence caller.
Usage: bcftools csq [options] in.vcf

Required options:
   -f, --fasta-ref <file>          reference file in fasta format
   -g, --gff-annot <file>          gff3 annotation file

CSQ options:
   -b, --brief-predictions         annotate with abbreviated protein-changing predictions
   -c, --custom-tag <string>       use this tag instead of the default BCSQ
   -l, --local-csq                 localized predictions, consider only one VCF record at a time
   -n, --ncsq <int>                maximum number of consequences to consider per site [16]
   -p, --phase <a|m|r|R|s>         how to handle unphased heterozygous genotypes: [r]
                                     a: take GTs as is, create haplotypes regardless of phase (0/1 -> 0|1)
                                     m: merge *all* GTs into a single haplotype (0/1 -> 1, 1/2 -> 1)
                                     r: require phased GTs, throw an error on unphased het GTs
                                     R: create non-reference haplotypes if possible (0/1 -> 1|1, 1/2 -> 1|2)
                                     s: skip unphased hets
Options:
   -e, --exclude <expr>            exclude sites for which the expression is true
       --force                     run even if some sanity checks fail
   -i, --include <expr>            select sites for which the expression is true
       --no-version                do not append version and command line to the header
   -o, --output <file>             write output to a file [standard output]
   -O, --output-type <b|u|z|v|t>   b: compressed BCF, u: uncompressed BCF, z: compressed VCF
                                   v: uncompressed VCF, t: plain tab-delimited text output [v]
   -r, --regions <region>          restrict to comma-separated list of regions
   -R, --regions-file <file>       restrict to regions listed in a file
   -s, --samples <-|list>          samples to include or "-" to apply all variants and ignore samples
   -S, --samples-file <file>       samples to include
   -t, --targets <region>          similar to -r but streams rather than index-jumps
   -T, --targets-file <file>       similar to -R but streams rather than index-jumps
       --threads <int>             use multithreading with <int> worker threads [0]
   -v, --verbose <int>             verbosity level 0-2 [1]

Example:
   bcftools csq -f hs37d5.fa -g Homo_sapiens.GRCh37.82.gff3.gz in.vcf

   # GFF3 annotation files can be downloaded from Ensembl. e.g. for human:
   ftp://ftp.ensembl.org/pub/current_gff3/homo_sapiens/
   ftp://ftp.ensembl.org/pub/grch37/release-84/gff3/homo_sapiens/
```



Téléchargez le fichier d'annotation gff pour le chromosome chr22:


Pierre LINDENBAUM's avatar
Pierre LINDENBAUM committed
814
Pour hg19:
Pierre LINDENBAUM's avatar
gff    
Pierre LINDENBAUM committed
815
816

```
Pierre LINDENBAUM's avatar
Pierre LINDENBAUM committed
817
wget -O annot.gff.gz "http://ftp.ensembl.org/pub/grch37/release-99/gff3/homo_sapiens/Homo_sapiens.GRCh37.87.chromosome.22.gff3.gz"
Pierre LINDENBAUM's avatar
gff    
Pierre LINDENBAUM committed
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
```

Les fichiers GFF  utilisent la notation '22' alors que nos chromosomes utilisent 'chr22'. Nous allons donc ajouter ce prefix avec `sed` pour les lignes qui commencent par '22'

```
gunzip -c annot.gff.gz | sed 's/^22/chr22/' | gzip > annotations.gff.gz
```

et nous annotons le fichier vcf avec le fichier de gènes GFF avec la command `bcftools csq`

```
bcftools csq --phase a --fasta-ref ref.fa --gff-annot annot.gff.gz -O b -o result.annot.vcf.gz result.vcf.gz
```

```

bcftools view result.annot.vcf.gz


(...)
chr22	22712467	.	T	G	730	.	DP=332;VDB=0.784239;SGB=-2.07845;RPB=0.40
8006;MQB=0.451512;MQSB=0.413822;BQB=0.0195238;MQ0F=0.0180723;ICB=0.128205;HOB=0.0555556;AC=5;AN=6;DP4=9,3
0,108,129;MQ=51;BCSQ=missense|IGLV1-47|ENST00000390294|IG_V|+|70S>70R|22712467T>G	GT:PL:DP:AD:GP:GQ
:BCSQ	1/1:255,255,0:107:2,105:269,259,0:127:3	0/1:255,0,255:81:36,45:265,0,250:127:2	1/1:255,245,0:88:
1,87:269,249,0:127:3
(...)
```
lindenb's avatar
1st  
lindenb committed
845

lindenb's avatar
cont    
lindenb committed
846

Pierre Lindenbaum's avatar
cont    
Pierre Lindenbaum committed
847
848
## Makefile

Solena LE SCOUARNEC's avatar
Solena LE SCOUARNEC committed
849
pour ceux qui souhaitent, à la fin, voir à quoi ressemble le fichier Makefile pour ce workflow, un workflow est disponible ici: [workflow.mk](workflow.mk)