Commit 94dbc9c6 authored by lindenb's avatar lindenb
Browse files

cont

parents 6372d9c9 4bea3a4b
# TP Next Generation Sequencing 2018
# 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 :
......@@ -6,14 +6,21 @@ Au cours de cette séance de TP nous allons chercher les mutations de-novo d'un
* extraire les mutations
* annoter ces variations.
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.
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
Les données sont téléchargeables à l'adresse suivante : http://filex.univ-nantes.fr/get?k=PbitjcCjCWnfH1KqwET
Les donnees sont telechargeable a l'URL suivante:
https://uncloud.univ-nantes.fr/index.php/s/oEYQfdJQcReB7ez
Téléchargez et extrayez le fichier zip avec
......@@ -21,7 +28,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.
......@@ -85,9 +92,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"
......@@ -95,7 +102,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"
......@@ -118,13 +125,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
......@@ -136,9 +143,10 @@ A l'aide de `cat`, concatenez ces deux fichiers dans un fichier que vous nommere
$ cat chr22.fa chrM.fa > ref.fa
```
## TP septembre 20220
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:
## TP septembre 2020
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_bwa_20200917/bin:${PATH}"
......@@ -337,10 +345,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-aa\tSM:SAMPLENAME' ref.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**.
......@@ -348,10 +356,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
......@@ -403,7 +411,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:
......@@ -487,6 +495,9 @@ Options:
Triez `child_unsorted.bam` pour produire un fichier `child.bam`.
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
```
samtools sort -T tmp -O BAM -o child.bam child_unsorted.bam
```
......@@ -536,11 +547,13 @@ 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:
```
samtools view child.bam chr22
samtools view child.bam chrM
# la region ci dessous ne donnera rien si il n'y a pas de read dans la region
samtools view child.bam chr22:23707000-23803000
```
......@@ -552,7 +565,7 @@ Samtools doit pouvoir accèder rapidement à des régions du génome de referenc
### Synopsis:
```
samtools faidx reference.fa
samtools faidx ref.fa
```
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.
......@@ -562,7 +575,7 @@ Une fois indexé, on peut rapidement obtenir des sous-régions du genome de ref
### Synopsis:
```
samtools faidx reference.fa chr:start-end
samtools faidx ref.fa chr:start-end
```
En utilisant `samtools faidx` extrayez la sequence `chrM:50-200`
......@@ -576,7 +589,7 @@ La commande `samtools tview` permet de visualiser interactivement les reads dans
### Synopsis:
```
samtools tview sorted.bam reference.fa
samtools tview sorted.bam ref.fa
```
A l'aide de `samtools view` affichez les premiers reads mappés dans `child.bam`, notez la position génomique du 1er read mappé.
......@@ -671,7 +684,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
```
......@@ -702,4 +715,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
"""
}
.PHONY: all zip
REG= "22:22707000-23803000" "MT:3910-4250"
BASE=
all: $(addprefix DATA/,child.R1.fq.gz mother.R1.fq.gz father.R1.fq.gz)
......
......@@ -22,10 +22,10 @@ result.vcf.gz.tbi: result.vcf.gz
tabix --force --preset vcf $<
result.vcf.gz: result.bcf
bcftools call --ploidy GRCh37 --multiallelic-caller --variants-only --output-type z -o $@ --format-fields GQ,GP $<
bcftools call --multiallelic-caller --variants-only --output-type z -o $@ --format-fields GQ,GP $<
result.bcf : $(addsuffix .bam.bai,child father mother) $(addsuffix .bam,child father mother) ref.fa.fai
samtools mpileup -g --uncompressed --output-tags DP --reference ref.fa --output $@ $(filter %.bam,$^)
samtools mpileup -g --uncompressed --output-tags DP --fasta-ref ref.fa --output $@ $(filter %.bam,$^)
$(eval $(call mapreads,father))
......
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