Commit c22b5cb5 authored by Pierre Lindenbaum's avatar Pierre Lindenbaum
Browse files

cont

parent 6f2eb08e
......@@ -8,7 +8,7 @@ Certains outils ont besoin de place; le cas échéant, faites de la place, videz
## Recupération des données FASTQ
Les données sont téléchargeables à l'adresse suivante : http://filex.univ-nantes.fr/get?k=EWZwzJg76wdtYkKcTC7
Les données sont téléchargeables à l'adresse suivante : http://filex.univ-nantes.fr/get?k=moyINN3cCurgIm1BFwb
Téléchargez et extrayez le fichier zip avec
......@@ -22,67 +22,60 @@ Dans le repertoire *DATA* se trouvent les fichiers FASTQ contenant les short-rea
```
DATA/child.R2.ab.fq.gz
DATA/child.R2.aa.fq.gz
DATA/father.R1.aa.fq.gz
DATA/father.R2.aa.fq.gz
DATA/mother.R2.aa.fq.gz
DATA/mother.R1.ab.fq.gz
DATA/mother.R2.ab.fq.gz
DATA/father.R1.ab.fq.gz
DATA/father.R2.ab.fq.gz
DATA/child.R1.aa.fq.gz
DATA/child.R1.ab.fq.gz
DATA/mother.R1.aa.fq.gz
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
```
Nous travaillons en **paired-end**: il y a donc un fichier **R1** et **R2** pour chaque individu. Les séquenceurs illumina
découpent les FASTQs en morceaux. Il y a donc des series **aa** et **bb** pour chaque individu.
Nous travaillons en **paired-end**: il y a donc un fichier **R1** et **R2** pour chaque individu.
* Comptez le nombre de lignes R1.aa chez l'enfant; Vérifiez que c'est un multiple de 4.
* Comptez le nombre de lignes R1 chez l'enfant; Vérifiez que c'est un multiple de 4.
```
gunzip -c DATA/child.R1.aa.fq.gz | wc -l
gunzip -c DATA/child.R1.fq.gz | wc -l
```
* Comptez le nombre de lignes R2.aa chez l'enfant; Pourquoi ce même nombre ?
* Comptez le nombre de lignes R2 chez l'enfant; Pourquoi ce même nombre ?
```
gunzip -c DATA/child.R2.aa.fq.gz | wc -l
gunzip -c DATA/child.R2.fq.gz | wc -l
```
* Comparez les 10 premiers noms de reads de l'enfant entre le fichier R1 et R2, quelle est la différence ?
```
$ gunzip -c DATA/child.R1.aa.fq.gz | paste - - - - | cut -f 1 | head
$ gunzip -c DATA/child.R2.aa.fq.gz | paste - - - - | cut -f 1 | head
$ gunzip -c DATA/child.R1.fq.gz | paste - - - - | cut -f 1 | head
$ gunzip -c DATA/child.R2.fq.gz | paste - - - - | cut -f 1 | head
```
* Quelle est la longueur des short-reads R1 chez l'enfant ?
```
$ gunzip -c DATA/child.R1.aa.fq.gz | paste - - - - | cut -f 2 | head -n 1 | tr -d '\n' | wc -c
$ gunzip -c DATA/child.R1.fq.gz | paste - - - - | cut -f 2 | head -n 1 | tr -d '\n' | wc -c
```
* 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... ?).
```
$ gunzip -c DATA/child.R1.aa.fq.gz | paste - - - - | cut -f 2 | cut -c 1| sort | uniq -c | sort -n
$ gunzip -c DATA/child.R1.fq.gz | paste - - - - | cut -f 2 | cut -c 1| sort | uniq -c | sort -n
```
* Quelles sont les bases que l'on trouve dans les séquences d'ADN, n'y a-t-il que A-T-G-C ?
```
$ gunzip -c DATA/child.R1.aa.fq.gz | paste - - - - | cut -f 2 | tr -d 'ATGC\n'
$ gunzip -c DATA/child.R1.fq.gz | paste - - - - | cut -f 2 | tr -d 'ATGC\n'
```
* Quels sont les caractères ascii de qualité que l'on trouve dans les lignes de *qualité* ?
```
$ gunzip -c DATA/child.R1.aa.fq.gz | paste - - - - | cut -f 4 | grep -o . | sort | uniq
$ gunzip -c DATA/child.R1.fq.gz | paste - - - - | cut -f 4 | grep -o . | sort | uniq
```
## Recupération du génome de reference
......@@ -295,29 +288,29 @@ l'option `-R` permet de créer un **READ-GROUP** pour associer les reads a un é
la sortie standard de cette commande est un fichier 'texte' **SAM**.
Mappez les reads de l'enfant `child.R1.aa.fq.gz` et `child.R2.aa.fq.gz` avec cette commande et redirigez la sortie standard vers un fichier `child_unsorted.aa.sam`
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`
```
bwa mem -R '@RG\tID:child-aa\tSM:child' ref.fa DATA/child.R1.aa.fq.gz DATA/child.R2.aa.fq.gz > child_unsorted.aa.sam
bwa mem -R '@RG\tID:child-aa\tSM:child' ref.fa DATA/child.R1.fq.gz DATA/child.R2.fq.gz > child_unsorted.sam
```
* que renvoie la commande
```
file child_unsorted.aa.sam
file child_unsorted.sam
```
* Combien y-a-t-il de lignes dans ce fichier sam ?
```
wc -l child_unsorted.aa.sam
wc -l child_unsorted.sam
```
* Combien y-a-t-il de lignes dans l'en-tête ?
```
grep '^@' child_unsorted.aa.sam | wc -l
grep '^@' child_unsorted.sam | wc -l
```
......@@ -326,7 +319,7 @@ grep '^@' child_unsorted.aa.sam | wc -l
* 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.
* Combien y-a-t-il d'alignements ? Comparez au nombre de reads dans les fichiers child.R1.aa.fq.gz et child.R2.aa.fq.gz
* Combien y-a-t-il d'alignements ? Comparez au nombre de reads dans les fichiers child.R1.fq.gz et child.R2.fq.gz
* Identifiez les colonnes d'un fichier SAM:
* READ NAME
......@@ -342,7 +335,7 @@ grep '^@' child_unsorted.aa.sam | wc -l
* QUALITIES
```
grep -v '^@' child_unsorted.aa.sam | head -n 1 | tr "\t" "\n" | cat -n
grep -v '^@' child_unsorted.sam | head -n 1 | tr "\t" "\n" | cat -n
```
* 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 ?
......@@ -353,7 +346,7 @@ grep -v '^@' child_unsorted.aa.sam | head -n 1 | tr "\t" "\n" | cat -n
## Conversion de SAM en BAM
pour permettre aux algorithmes d'aller plus vite, nous utilisons samtools pour convertir le fichier en format texte/SAM `child.aa.sam` en format binaire/BAM `child_unsorted.aa.bam`
pour permettre aux algorithmes d'aller plus vite, nous utilisons samtools pour convertir le fichier en format texte/SAM `child.sam` en format binaire/BAM `child_unsorted.bam`
### Synopsis:
......@@ -366,13 +359,13 @@ samtools view -bS -o out.bam in.sam
* Option **-b** : la sortie est du binaire (**BAM**)
* Option **-o** : nom du fichier de sortie
Generez le fichier `child_unsorted.aa.bam`
Generez le fichier `child_unsorted.bam`
```
samtools view -bS -o child_unsorted.aa.bam child_unsorted.aa.sam
samtools view -bS -o child_unsorted.bam child_unsorted.sam
```
* que renvoie la commande `file child_unsorted.aa.bam` ?
* que renvoie la commande `file child_unsorted.bam` ?
* affichez les options de `samtools view` avec
```
......@@ -382,22 +375,22 @@ samtools view
* affichez le corps du BAM
```
samtools view child_unsorted.aa.bam
samtools view child_unsorted.bam
```
* affichez l'en-tete du BAM avec
```
samtools view -H child_unsorted.aa.bam
samtools view -H child_unsorted.bam
```
* affichez le bam complet avec
```
samtools view -h child_unsorted.aa.bam
samtools view -h child_unsorted.bam
```
* 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.aa.fq.gz (a.k.a 'first in pair')
* 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')
## Tri du BAM sur la position génomique
......@@ -416,48 +409,26 @@ samtools sort -o out.bam -T out in-unsorted.bam
Triez `child_unsorted.aa.bam` pour produire un fichier `child_sorted.aa.bam`.
Triez `child_unsorted.bam` pour produire un fichier `child.bam`.
```
samtools sort -T tmp -o child_sorted.aa.bam child_unsorted.aa.bam
samtools sort -T tmp -o child.bam child_unsorted.bam
```
* 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`.
```
samtools view -H child_sorted.aa.bam
samtools view -H child.bam
```
* En utilisant 'samtools view', Verifiez que le fichier est trié en observant la position croissante CHROM/POS des reads dans le fichier BAM.
```
samtools view child_sorted.aa.bam | cut -f 3,4
samtools view child.bam | cut -f 3,4
```
## Mapping de la deuxième paire child.R1.ab.fq.gz et child.R2.ab.fq.gz
de la même manière que vous aviez mappé les fichiers child.R1.aa.fq.gz et child.R2.aa.fq.gz, mappez les fichiers **child.R1.ab.fq.gz** et **child.R2.ab.fq.gz** de manière à produire un fichier BAM **child_sorted.ab.bam**
Note: Dans un environnement parallélisé (cluster...), le mapping de **child_sorted.aa.bam** et **child_sorted.ab.bam** pourrait être fait en parallele...
## Fusion des BAM pour le même échantillon.
Maintenant que nous avons **child_sorted.aa.bam** et **child_sorted.ab.bam**, nous pouvons les fusionner en un seul fichier BAM **child.bam** pour cela on utilise la command `samtools merge`
### Synopsis:
```
samtools merge -f out.bam in1.bam in2.bam in3.bam ....
# option `-f` : 'forcer' la création du fichier.
```
mergez les fichiers **child_sorted.aa.bam** et **child_sorted.ab.bam** pour produire **child.bam**. A l'aide de `samtools view`, vérifiez que le fichier est toujours trié et que le nombre d'alignements total est conservé.
```
samtools merge -f child.bam child_sorted.aa.bam child_sorted.ab.bam
```
## Indexation du BAM
......@@ -572,10 +543,8 @@ CAGGTAGGAAAAGGCCTAACAGACCAGGGCCAGCACCTCAGCCTGATTCTGACTCTTCTGGCAAAGATCC
Utilisez
* DATA/mother.R1.aa.fq.gz
* DATA/mother.R1.ab.fq.gz
* DATA/mother.R2.aa.fq.gz
* DATA/mother.R2.ab.fq.gz
* DATA/mother.R1.fq.gz
* DATA/mother.R2.fq.gz
et créez le fichier **mother.bam**
......@@ -583,10 +552,8 @@ et créez le fichier **mother.bam**
Utilisez
* DATA/father.R1.aa.fq.gz
* DATA/father.R1.ab.fq.gz
* DATA/father.R2.aa.fq.gz
* DATA/father.R2.ab.fq.gz
* DATA/father.R1.fq.gz
* DATA/father.R2.fq.gz
et créez le fichier **father.bam**
......
java.exe=/cm/shared/apps/java/jdk1.7.0_60/bin/java
java.exe=java
BAMS= /commun/data/projects/plateforme/NTS-015_ENS_Bezieau_HUGODIMS/20150604/align/Samples/10-01-BL-E/BAM/20150604_10-01-BL-E_final.bam \
/commun/data/projects/plateforme/NTS-015_ENS_Bezieau_HUGODIMS/20150604/align/Samples/10-01-BL-M/BAM/20150604_10-01-BL-M_final.bam \
/commun/data/projects/plateforme/NTS-015_ENS_Bezieau_HUGODIMS/20150604/align/Samples/10-01-BL-P/BAM/20150604_10-01-BL-P_final.bam
......@@ -8,57 +8,34 @@ BAMS= /commun/data/projects/plateforme/NTS-015_ENS_Bezieau_HUGODIMS/20150604/ali
REG= 22:22707000-23803000 MT:3910-4250
SPLIT=48000
all: $(addprefix DATA/,child.R1.aa.fq.gz mother.R1.aa.fq.gz father.R1.aa.fq.gz)
all: $(addprefix DATA/,child.R1.fq.gz mother.R1.fq.gz father.R1.fq.gz)
zip: DATA.zip
DATA.zip: $(addprefix DATA/,child.R1.aa.fq.gz mother.R1.aa.fq.gz father.R1.aa.fq.gz)
DATA.zip: $(addprefix DATA/,child.R1.fq.gz mother.R1.fq.gz father.R1.fq.gz)
rm -f $@
zip $@ -r DATA
DATA/child.R1.aa.fq.gz :
DATA/child.R1.fq.gz :
mkdir -p DATA
samtools view -b -o jeter.bam /commun/data/projects/plateforme/NTS-015_ENS_Bezieau_HUGODIMS/20150604/align/Samples/10-01-BL-E/BAM/20150604_10-01-BL-E_final.bam ${REG}
${java.exe} -jar /commun/data/packages/picard/picard-tools-1.126/picard.jar \
SamToFastq I=jeter.bam F=child.R1.fastq.gz F2=child.R2.fastq.gz FU=jeter.fq.gz VALIDATION_STRINGENCY=LENIENT
gunzip -c child.R1.fastq.gz | split -l ${SPLIT} - child.R1.
gunzip -c child.R2.fastq.gz | split -l ${SPLIT} - child.R2.
mv child.R1.aa child.R1.aa.fq
mv child.R1.ab child.R1.ab.fq
mv child.R2.aa child.R2.aa.fq
mv child.R2.ab child.R2.ab.fq
gzip -f child.R1.aa.fq child.R1.ab.fq child.R2.aa.fq child.R2.ab.fq
rm -f jeter.fq.gz jeter.bam child.R1.fastq.gz child.R2.fastq.gz
mkdir -p DATA
mv $(addsuffix .gz,child.R1.aa.fq child.R1.ab.fq child.R2.aa.fq child.R2.ab.fq) DATA/
SamToFastq I=jeter.bam F=DATA/child.R1.fq.gz F2=DATA/child.R2.fq.gz FU=jeter.fq.gz VALIDATION_STRINGENCY=LENIENT
rm -f jeter.fq.gz jeter.bam
DATA/mother.R1.aa.fq.gz :
DATA/mother.R1.fq.gz :
mkdir -p DATA
samtools view -b -o jeter.bam /commun/data/projects/plateforme/NTS-015_ENS_Bezieau_HUGODIMS/20150604/align/Samples/10-01-BL-M/BAM/20150604_10-01-BL-M_final.bam ${REG}
${java.exe} -jar /commun/data/packages/picard/picard-tools-1.126/picard.jar \
SamToFastq I=jeter.bam F=mother.R1.fastq.gz F2=mother.R2.fastq.gz FU=jeter.fq.gz VALIDATION_STRINGENCY=LENIENT
gunzip -c mother.R1.fastq.gz | split -l ${SPLIT} - mother.R1.
gunzip -c mother.R2.fastq.gz | split -l ${SPLIT} - mother.R2.
mv mother.R1.aa mother.R1.aa.fq
mv mother.R1.ab mother.R1.ab.fq
mv mother.R2.aa mother.R2.aa.fq
mv mother.R2.ab mother.R2.ab.fq
gzip -f mother.R1.aa.fq mother.R1.ab.fq mother.R2.aa.fq mother.R2.ab.fq
rm -f jeter.fq.gz jeter.bam mother.R1.fastq.gz mother.R2.fastq.gz
mkdir -p DATA
mv $(addsuffix .gz,mother.R1.aa.fq mother.R1.ab.fq mother.R2.aa.fq mother.R2.ab.fq) DATA/
SamToFastq I=jeter.bam F=DATA/mother.R1.fq.gz F2=DATA/mother.R2.fq.gz FU=jeter.fq.gz VALIDATION_STRINGENCY=LENIENT
rm -f jeter.fq.gz jeter.bam
DATA/father.R1.aa.fq.gz :
DATA/father.R1.fq.gz :
mkdir -p DATA
samtools view -b -o jeter.bam /commun/data/projects/plateforme/NTS-015_ENS_Bezieau_HUGODIMS/20150604/align/Samples/10-01-BL-P/BAM/20150604_10-01-BL-P_final.bam ${REG}
${java.exe} -jar /commun/data/packages/picard/picard-tools-1.126/picard.jar \
SamToFastq I=jeter.bam F=father.R1.fastq.gz F2=father.R2.fastq.gz FU=jeter.fq.gz VALIDATION_STRINGENCY=LENIENT
gunzip -c father.R1.fastq.gz | split -l ${SPLIT} - father.R1.
gunzip -c father.R2.fastq.gz | split -l ${SPLIT} - father.R2.
mv father.R1.aa father.R1.aa.fq
mv father.R1.ab father.R1.ab.fq
mv father.R2.aa father.R2.aa.fq
mv father.R2.ab father.R2.ab.fq
gzip -f father.R1.aa.fq father.R1.ab.fq father.R2.aa.fq father.R2.ab.fq
rm -f jeter.fq.gz jeter.bam father.R1.fastq.gz father.R2.fastq.gz
mkdir -p DATA
mv $(addsuffix .gz,father.R1.aa.fq father.R1.ab.fq father.R2.aa.fq father.R2.ab.fq) DATA/
SamToFastq I=jeter.bam F=DATA/father.R1.fq.gz F2=DATA/father.R2.fq.gz FU=jeter.fq.gz VALIDATION_STRINGENCY=LENIENT
rm -f jeter.fq.gz jeter.bam
1000G_omni2.5.b37.vcf 1000G_phase1.snps.high_confidence.b37.vcf Mills_and_1000G_gold_standard.indels.b37.vcf :
cat /commun/data/pubdb/broadinstitute.org/bundle/1.5/b37/$@ | grep -E '(^##INFO|^##fileformat|^##FORMAT|^##FILTER|##contig=<ID=22|##contig=<ID=MT|^22|^#CHROM)' > $@
......
define mapreads2
define mapreads
$(1).bam.bai : $(1).bam
samtools index $$<
$(1).$(2).sorted.bam : $(1).$(2).unsorted.bam
$(1).bam : $(1).unsorted.bam
samtools sort -o $$@ -T $$(basename $$@) $$<
$(1).$(2).unsorted.bam :$(1).$(2).sam
$(1).unsorted.bam :$(1).sam
samtools view -Sb -o $$@ $$<
$(1).$(2).sam : DATA/$(1).R1.$(2).fq.gz DATA/$(1).R2.$(2).fq.gz ref.fa.bwt
$(1).sam : DATA/$(1).R1.fq.gz DATA/$(1).R2.fq.gz ref.fa.bwt
bwa mem -R '@RG\tID:$$(basename $$@)\tSM:$(1)' ref.fa $$(filter %.fq.gz,$$^) > $$@
endef
define mapreads
$(1).bam.bai : $(1).bam
samtools index $$<
$(1).bam: $(1).aa.sorted.bam $(1).ab.sorted.bam
samtools merge -f $$@ $$^
$(call mapreads2,$(1),aa)
$(call mapreads2,$(1),ab)
endef
......
workflow.png

246 KB | W: | H:

workflow.png

141 KB | W: | H:

workflow.png
workflow.png
workflow.png
workflow.png
  • 2-up
  • Swipe
  • Onion skin
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment