moussa moussaoui

Étape 1 sur 8

Lire un fichier VCF

Sauter les méta-lignes, parser l’en-tête, charger les données

Étape 1 / 8

Lire un fichier VCF

★★★

Contexte

Le VCF (Variant Call Format) est le format universel pour stocker les variants génétiques (SNVs, indels, structuraux). Il est tabulaire mais précédé d’un nombre variable de méta-lignes qui commencent par ##.

les lignes ## décrivent les champs (versions, INFO, FILTER, FORMAT) — à ignorer pour l’analyse,
la ligne #CHROM est l’en-tête tabulaire,
chaque ligne suivante est un variant (8 colonnes obligatoires + colonnes échantillons).

Sauvegardez ce VCF de démonstration sous variants.vcf :

##fileformat=VCFv4.2
##INFO=<ID=DP,Number=1,Type=Integer,Description="Read depth">
##INFO=<ID=AF,Number=1,Type=Float,Description="Allele frequency">
##INFO=<ID=GENE,Number=1,Type=String,Description="Gene symbol">
#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO
chr1	1500	rs1	A	G	35.5	PASS	DP=120;AF=0.45;GENE=BRCA1
chr1	2400	rs2	C	T	42.1	PASS	DP=98;AF=0.51;GENE=BRCA1
chr1	3100	rs3	G	A	18.0	LowQual	DP=22;AF=0.40;GENE=BRCA1
chr1	5200	rs4	T	TAC	50.0	PASS	DP=150;AF=0.48;GENE=TP53
chr2	1100	rs5	G	C	28.0	PASS	DP=85;AF=0.008;GENE=KRAS
chr2	2200	rs6	AGT	A	33.0	PASS	DP=110;AF=0.50;GENE=KRAS
chr3	1500	rs7	C	G	45.0	PASS	DP=180;AF=0.95;GENE=MYC
chr3	3300	rs8	T	C	29.0	PASS	DP=70;AF=0.30;GENE=MYC
chrX	5500	rs9	G	T	15.0	LowQual	DP=18;AF=0.60;GENE=AR
chrX	8800	rs10	A	C	38.5	PASS	DP=130;AF=0.49;GENE=AR

Travail à réaliser

Écrivez lire_vcf(chemin) qui renvoie une liste de dicts, un par variant :

{
  "chrom":  "chr1",
  "pos":    1500,
  "id":     "rs1",
  "ref":    "A",
  "alt":    "G",
  "qual":   35.5,
  "filter": "PASS",
  "info":   "DP=120;AF=0.45;GENE=BRCA1",
}

Affichez ensuite le nombre de variants chargés + la première entrée.

10 variants chargés.
Premier : chr1:1500 A>G  QUAL=35.5  FILTER=PASS

Étape 2 sur 8

Distribution par chromosome

Combien de variants sur chr1, chr2, chrX ?

Étape 2 / 8

Distribution par chromosome

★★★

Contexte

Première inspection : combien de variants par chromosome ? Ce simple comptage révèle déjà beaucoup. Un excès sur un chromosome donné peut indiquer une région mal mappée, un excès en duplications de référence, ou une biologie réelle (chromosome impliqué dans la maladie étudiée).

Travail à réaliser

Écrivez compter_par_chrom(variants) qui renvoie un dict {chrom: nombre}, puis affichez un mini-tableau trié par ordre de fréquence décroissante.

Variants par chromosome :
  chr1  : 4
  chr3  : 2
  chr2  : 2
  chrX  : 2

Étape 3 sur 8

SNV ou indel ?

Classer chaque variant par type

Étape 3 / 8

SNV ou indel ?

★★★

Contexte

Un SNV (Single Nucleotide Variant) change une seule base : REF=A, ALT=G. Un indel ajoute ou supprime des bases : REF=T, ALT=TAC est une insertion ; REF=AGT, ALT=A est une délétion.

Sur le plan clinique, les indels sont souvent plus délétères : un frameshift dans une région codante détruit toute la protéine en aval.

Travail à réaliser

Écrivez type_variant(v) qui renvoie "SNV", "insertion", "deletion", ou "complexe".

Ensuite tableau_types(variants) qui résume :

Types de variants :
  SNV         : 8
  insertion   : 1
  deletion    : 1
  complexe    : 0

Étape 4 sur 8

Ratio Ts/Tv

Transitions vs transversions — un contrôle qualité

Étape 4 / 8

Ratio Ts/Tv

★★★

Contexte

Parmi les SNV, on distingue :

Transitions (Ts) : purine ↔ purine ou pyrimidine ↔ pyrimidine. Soit A↔G, soit C↔T.
Transversions (Tv) : purine ↔ pyrimidine. Toutes les autres (A↔C, A↔T, G↔C, G↔T).

Le rapport Ts/Tv est un contrôle qualité classique : il devrait être autour de 2.0–2.1 pour les exomes humains et autour de 2.0 globalement en génome entier. Un ratio anormalement bas (~1) trahit des faux positifs (erreurs de séquençage), un ratio anormalement haut indique un biais.

Travail à réaliser

Écrivez ratio_ts_tv(variants) qui :

filtre les SNVs (les indels ne comptent pas),
compte transitions et transversions,
renvoie (n_ts, n_tv, ratio).

Transitions   : 5
Transversions : 3
Ts/Tv         : 1.67  (attendu ~2.0 sur exome humain)

Étape 5 sur 8

Parser le champ INFO

Du semicolon-soup à un dict propre

Étape 5 / 8

Parser le champ INFO

★★★

Contexte

Le champ INFO contient les métadonnées importantes du variant : profondeur (DP), fréquence allélique (AF), gène impacté, effet fonctionnel, score CADD… Tout est encodé en cle=valeur séparé par ; :

DP=120;AF=0.45;GENE=BRCA1

Pour faire quoi que ce soit d’utile, il faut convertir ce salmigondis en dict Python.

Travail à réaliser

Écrivez parser_info(s) qui renvoie un dict. Une clé sans = (drapeau booléen) est associée à True.

Puis enrichissez chaque variant avec ce dict via enrichir(variants) (modifie sur place ou renvoie une copie — au choix).

rs1 INFO : {'DP': '120', 'AF': '0.45', 'GENE': 'BRCA1'}
rs1 DP   : 120
rs1 AF   : 0.45
rs1 GENE : BRCA1

Étape 6 sur 8

Filtrage qualité

FILTER, QUAL et DP — garder le bon grain

Étape 6 / 8

Filtrage qualité

★★★

Contexte

Tous les variants détectés ne se valent pas. Un appel à QUAL=15 et DP=18 est probablement un artefact ; à QUAL=42 et DP=120, on est beaucoup plus confiant. Les pipelines standard rejettent :

FILTER ≠ "PASS" (l’appeleur a déjà tagué un problème),
QUAL < 20 (~ 1 chance sur 100 d’erreur),
DP < 30 (profondeur insuffisante).

Travail à réaliser

Écrivez filtrer_qualite(variants, qual_min=20, dp_min=30) qui renvoie deux listes : (garde, rejete). Affichez un mini-rapport :

Filtre : QUAL >= 20 ET DP >= 30 ET FILTER=PASS
Gardés  : 7 / 10
Rejetés : 3
  rs3   chr1:3100   QUAL=18.0  DP=22   FILTER=LowQual
  rs5   chr2:1100   QUAL=28.0  DP=85   FILTER=PASS   -> mais AF basse
  rs9   chrX:5500   QUAL=15.0  DP=18   FILTER=LowQual

(la justification des rejets dépend de votre implémentation ; le format ci-dessus est indicatif).

Étape 7 sur 8

Distribution des allèles (AF)

Rare, peu fréquent, commun

Étape 7 / 8

Distribution des allèles (AF)

★★★

Contexte

La fréquence allélique (AF) — proportion d’allèles alternatifs dans la population — classe les variants en :

Rares : AF < 0.01 (1 %). Souvent les plus intéressants en médecine personnalisée.
Peu fréquents : 0.01 ≤ AF < 0.05.
Communs : AF ≥ 0.05. Polymorphismes courants (SNPs neutres pour la plupart).

Cette stratification est la première étape de toute analyse de gène candidat.

Travail à réaliser

Écrivez classer_af(variants) qui renvoie un Counter {"rare": n, "peu_frequent": n, "commun": n}. Puis affichez :

Distribution AF :
  rare         (AF < 0.01)  : 1
  peu_frequent (0.01-0.05)   : 0
  commun       (AF >= 0.05)  : 9

Variants rares (à investiguer en priorité) :
  rs5  chr2:1100  G>C  GENE=KRAS  AF=0.008

Étape 8 sur 8

Pipeline CLI : vcfqc.py

Rapport complet en une commande

Étape 8 / 8

Pipeline CLI : vcfqc.py

★★★

Contexte

Sept briques posées. Reste à les enchaîner dans un script utilisable. L’objectif n’est pas de remplacer vcftools ou bcftools stats, mais de comprendre comment ils sont construits — et d’avoir un outil de QC à vous que vous pouvez modifier.

Travail à réaliser

Rassemblez tout dans vcfqc.py. main() doit produire :

résumé global : N variants, par chromosome,
types : SNV / insertion / délétion / complexe,
Ts/Tv,
après filtrage qualité : N gardés, N rejetés,
distribution AF + liste des variants rares.

$ python vcfqc.py variants.vcf

Analyse de variants (VCF)