Commentaires 1
c'est quoi cette ligne ? elle sert à quoi ?
je prefererais que ça soit un argument de ton script R.
c'est quoi ALT1 vs ALT2. Pour les dialleliques ? comment utiliser cette information ? comment feras tu pour les multialleliques > 2?
Pour ce problème (connu) des multialleliques, tu peux considerer que dès qu'un allele n'est pas REF, il s'agit du même allele rare ou bien tu peux ajouter une normalisation bcftools norm
pour convertir les multialleliques en dialleliques (les deux solutions sont problèmatiques)
C'est un bon commencement pour lire les VCFs mais rapidement vas re-inventer un parser de VCF, il faut vraiment que tu utises une API existante pour lire le VCF. Je ne connais pas celles en R.
Globalement tu as la bonne idée pour scanner le vcf (scanner les genotypes, compter les genotypes qui portent l'allele rare).
commentaires en anglais STP. si le script est publié ça t'évitera dès maintenant à tout reprendre tes scripts...
par exemple là ton script ne marcherait pas si il y avait un genotype 10/11
même idée ici si le génotype est 2/3