Commit c7f9b32d authored by lindenb's avatar lindenb
Browse files

cont

parents 0de83e6c 044da162
# TP Next Generation Sequencing 2017
# TP Next Generation Sequencing 2019
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.
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 et extraire les mutations.
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)
Certains outils ont besoin de place; le cas échéant, faites de la place, videz vos répertoires temporaires.
La version de l'outil samtools actuellement installée à la fac est une vieille version (1.2). Elle est suffisante pour ce TP mais les commandes avec la nouvelles version de samtools sera légerement différent.
La version de l'outil samtools actuellement installée à la fac est une vieille version (1.2). Elle est suffisante pour ce TP mais les commandes avec la nouvelle version de samtools seront légerement différentes.
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.
## Recupération des données FASTQ
......@@ -20,7 +22,7 @@ Téléchargez et extrayez le fichier zip avec
unzip DATA.zip
```
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*).
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*).
Dans le repertoire *DATA* se trouvent les fichiers FASTQ contenant les short-reads pour father/mother/child.
......@@ -84,9 +86,9 @@ $ gunzip -c DATA/child.R1.fq.gz | paste - - - - | cut -f 4 | grep -o . | sort |
## Recupération du génome de reference
Pour aller vite, nous nous contenterons des chromosomes 22 et et de l'ADN mitochondrial humains.
Pour aller vite, nous nous contenterons des chromosomes 22 et de l'ADN mitochondrial humains.
* Si votre nom de famille commence avant la lettre 'L': Telechargez les séquences fa.gz **chr22** et **chrM** de la version **hg38/GRCh38** du genome humain .
* Si votre nom de famille commence avant la lettre 'L': Telechargez les séquences fa.gz **chr22** et **chrM** de la version **hg38/GRCh38** du genome humain.
```
$ wget -O chr22.fa.gz "http://hgdownload.cse.ucsc.edu/goldenPath/hg38/chromosomes/chr22.fa.gz"
......@@ -94,7 +96,7 @@ $ wget -O chrM.fa.gz "http://hgdownload.cse.ucsc.edu/goldenPath/hg38/chromosomes
```
* Pour les autres, telechargze les séquences fa.gz **chr22** et **chrM** de la version **hg19/GRCh37** du genome humain.
* Pour les autres, telechargez les séquences fa.gz **chr22** et **chrM** de la version **hg19/GRCh37** du genome humain.
```
$ wget -O chr22.fa.gz "http://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/chr22.fa.gz"
......@@ -117,13 +119,13 @@ $ grep -v '>' chr22.fa | tr -d '\n' | wc -c
$ grep -v '>' chrM.fa | tr -d '\n' | wc -c
```
* Quelles sont les premières bases du chromosomes 22 ?
* Quelles sont les premières bases du chromosome 22 ?
```
$ more chr22.fa
```
Pourquoi observez vous cela ? Combien de lignes y-a-t-il avant que cela change ?
Pourquoi observez vous cela ? Combien de lignes y a-t-il avant que cela change ?
```
$ cat -n chr22.fa | grep -v '>' | grep -v N -m1
......@@ -135,9 +137,9 @@ A l'aide de `cat`, concatenez ces deux fichiers dans un fichier que vous nommere
$ cat chr22.fa chrM.fa > ref.fa
```
## TP du 29 septembre 2017
## TP septembre 2018
A priori tous les outils ont-été installé. 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:
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:
```
export PATH="/usr/local/opt/sam_vcf_htslib/bin:${PATH}"
......@@ -284,10 +286,10 @@ maintenant que l'index du génome de référence a été créé, nous cherchons
le synopsis de la commande bwa pour mapper les reads est :
```
bwa mem -R '@RG\tID:SAMPLENAME-aa\tSM:SAMPLENAME' reference.fa sample.R1.gz sample.R2.gz > sample.sam
bwa mem -R '@RG\tID:SAMPLENAME\tSM:SAMPLENAME' reference.fa sample.R1.gz sample.R2.gz > sample.sam
```
l'option `-R` permet de créer un **READ-GROUP** pour associer les reads a un échantillon (donc ici, remplacez SAMPLENAME par 'child').
l'option `-R` permet de créer un **READ-GROUP** pour associer les reads à un échantillon (donc ici, remplacez SAMPLENAME par 'child').
la sortie standard de cette commande est un fichier 'texte' **SAM**.
......@@ -295,10 +297,10 @@ la sortie standard de cette commande est un fichier 'texte' **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.fq.gz DATA/child.R2.fq.gz > child_unsorted.sam
bwa mem -R '@RG\tID:child\tSM:child' ref.fa DATA/child.R1.fq.gz DATA/child.R2.fq.gz > child_unsorted.sam
```
* que renvoie la commande
* Que renvoie la commande suivante ?
```
file child_unsorted.sam
......@@ -350,7 +352,7 @@ grep -v '^@' child_unsorted.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.sam` en format binaire/BAM `child_unsorted.bam`
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`
### Synopsis:
......@@ -445,7 +447,7 @@ L'indexation consiste à créer un fichier d'index **'.bai'** à partir d'un bam
samtools index in.bam
```
Indexez le fichier**child.bam**, vérifiez le fonctionnement de l'index en effectuant des requêtes par région:
Indexez le fichier **child.bam**, vérifiez le fonctionnement de l'index en effectuant des requêtes par région:
```
samtools view child.bam chr22
......@@ -580,7 +582,7 @@ $ samtools mpileup -g --output-tags DP --fasta-ref ref.fa --output mutations.bcf
puis:
```
$ bcftools call --multiallelic-caller --variants-only --output-type z -o result.vcf.gz --format-fields GQ,GP input.bcf
$ bcftools call --multiallelic-caller --variants-only --output-type z -o result.vcf.gz --format-fields GQ,GP mutations.bcf
```
......@@ -611,4 +613,4 @@ Si il vous reste du temps, cherchez dans le fichier VCF les variations de-novo c
## Makefile
pour ceux qui souhaitent , à la fin, voir à quoi ressemble le fichier Makefile pour ce workflow, un workflow est disponible ici: [workflow.mk](workflow.mk)
pour ceux qui souhaitent, à la fin, voir à quoi ressemble le fichier Makefile pour ce workflow, un workflow est disponible ici: [workflow.mk](workflow.mk)
all: bin/nextflow nextflow.nf
bin/nextflow run -resume nextflow.nf --datadir ../../DATA
bin/nextflow:
mkdir -p bin && cd bin && wget -O - -q "https://get.nextflow.io" | bash && ./nextflow
datain = Channel.from("child","father","mother").
map{ S->[S,file("${params.datadir}/"+S+".R1.fq.gz"),file("${params.datadir}/"+S+".R2.fq.gz")] }
process dowmnloadref {
output:
file("ref.fa") into ref
file("ref.fa.*") into ref_idx
script:
"""
export PATH="/usr/local/opt/sam_vcf_htslib/bin:${PATH}"
wget -O - -q "http://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/chr22.fa.gz" | gunzip -c > ref.fa
wget -O - -q "http://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/chrM.fa.gz" | gunzip -c >> ref.fa
samtools faidx ref.fa
bwa index ref.fa
"""
}
process map {
tag "map ${sample}"
input:
set sample,R1,R2 from datain
file ref from ref
file refidx from ref_idx
output:
set sample,file("${sample}.sam") into sams
script:
"""
export PATH="/usr/local/opt/sam_vcf_htslib/bin:${PATH}"
bwa mem -R '@RG\\tID:${sample}\\tSM:${sample}' "${ref}" "${R1}" "${R2}" > "${sample}.sam"
"""
}
process sam2bam {
tag "sam2bam for ${sample}"
input:
set sample,sam from sams
output:
set sample,file("${sample}.unsorted.bam") into unsorted
script:
"""
export PATH="/usr/local/opt/sam_vcf_htslib/bin:${PATH}"
samtools view -b -S "${sam}" > "${sample}.unsorted.bam"
"""
}
process sortbam {
tag "sort ${sample}"
input:
set sample,bam from unsorted
output:
set sample,file("${sample}.bam") into sorted
file("${sample}.bam.bai")
script:
"""
export PATH="/usr/local/opt/sam_vcf_htslib/bin:${PATH}"
samtools sort -T ${sample}.tmp -o "${sample}.bam" ${bam}
samtools index ${sample}.bam
"""
}
process callvariants {
input:
val bams from sorted.map{T->T[1]}.collect()
file ref from ref
file refidx from ref_idx
output:
file("variants.vcf.gz") into vcf
script:
"""
export PATH="/usr/local/opt/sam_vcf_htslib/bin:${PATH}"
samtools mpileup -g --output-tags DP --fasta-ref "${ref}" --output variants.bcf ${bams.join(" ")}
bcftools call --multiallelic-caller --variants-only --output-type z -o variants.vcf.gz --format-fields GQ,GP variants.bcf
"""
}
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