J'ai des problèmes avec un script que j'ai écrit. Je l'ai décomposé pour identifier le problème et le voici.
Fichier d'entrée (délimité par des tabulations):
GT:AD:DP:GQ:PL 0/1 0/0 0/0 GT:AD:DP:GQ:PGT:PID:PL 0/1 0/0 0/0 GT:AD:DP:GQ:PGT:PID:PL 0/1 0/0 0/0 GT:AD:DP:GQ:PL 0/1 0/0 0/0 GT:AD:DP:GQ:PGT:PID:PL 0/0 0/0 0/1 GT:AD:DP:GQ:PGT:PID:PL 0/0 0/0 0/1 GT:AD:DP:GQ:PGT:PID:PL 0/0 0/0 0/1 GT:AD:DP:GQ:PGT:PID:PL 0/0 0/0 0/1 GT:AD:DP:GQ:PL ./. 0/0 0/1
Sortie souhaitée (prenez les 3 premiers caractères avant le deux points de chaque échantillon) par ligne et imprimez chaque ligne:
#!/usr/bin/perl use strict; my $inputfile1 = $ARGV[0]; open (FILE1, $inputfile1) or die "Uh oh.. unable to find file $inputfile1"; ##Opens input file my @file1 = <FILE1>; #loads inputfile1 data into array close FILE1; my (@colsplit, @genotypes1, @genotypes2, @genotypes3, @joined); foreach my $line(@file1) { ## process each line, splitting columns and move onto next line @colsplit = split("\t", $line); push (@joined, $colsplit[0]); foreach my $lines(@colsplit) { if ($colsplit[1] =~ m/(^0\/1)/ || $colsplit[1] =~ m/(^0\/0)/ || $colsplit[1]=~ m/(^1\/0)/ || $colsplit[1] =~ m/(^1\/1)/ || $colsplit[1] =~ m/(^.\/.)/) { push (@genotypes1, $1); } if ($colsplit[2] =~ m/(^0\/1)/ || $colsplit[2] =~ m/(^0\/0)/ || $colsplit[2] =~ m/(^1\/0)/ || $colsplit[2] =~ m/(^1\/1)/ || $colsplit[2] =~ m/(^.\/.)/) { push (@genotypes2, $1); } if ($colsplit[3] =~ m/(^0\/1)/ || $colsplit[3] =~ m/(^0\/0)/ || $colsplit[3] =~ m/(^1\/0)/ || $colsplit[3] =~ m/(^1\/1)/ || $colsplit[3] =~ m/(^.\/.)/) { push (@genotypes3, $1); } } } my $i = 0; foreach my $line(@joined) { if ($line =~ m/GT/) { print "$line\t$genotypes1[$i]\t$genotypes2[$i]\t$genotypes3[$i]\n"; $i++; }}
Le code que j'utilise pour faire cette étape ne produit pas le bon 0/0, 0/1, 0/2 codes comme prévu par ligne. Je pense que c'est un problème avec la façon dont j'ai écrit la boucle for, mais je ne suis pas sûr.
GT:AD:DP:GQ:PL 0/1 0/0 0/0 GT:AD:DP:GQ:PGT:PID:PL 0/0 0/0 0/1 GT:AD:DP:GQ:PGT:PID:PL ./. 0/0 0/1 GT:AD:DP:GQ:PL ./. 1/1 ./. GT:AD:DP:GQ:PGT:PID:PL 1/1 0/0 0/0 GT:AD:DP:GQ:PGT:PID:PL 1/1 0/0 0/0 GT:AD:DP:GQ:PGT:PID:PL 1/1 0/0 0/0 GT:AD:DP:GQ:PGT:PID:PL 1/1 0/0 0/0 GT:AD:DP:GQ:PL 1/1 0/1 ./.
Je pense que le problème peut être qu'une fois la première colonne sample1 correspondant, il saute à la deuxième ligne plutôt que d'itérer dans la deuxième colonne sample2. Je ne vois pas comment je gâche ça autrement! Cela me rend fou!
Ma sortie actuelle est:
FORMAT Sample1 Sample2 Sample3 GT:AD:DP:GQ:PL 0/1:17,6:23:85:85,0,370 0/0:51,6:57:17:0,17,1359 0/0:3,0:3:9:0,9,99 GT:AD:DP:GQ:PGT:PID:PL 0/0:3,0:3:0:.:.:0,0,38 0/0:1,0:1:3:.:.:0,3,33 0/1:1,2:3:26:0|1:13813_T_G:81,0,26 GT:AD:DP:GQ:PGT:PID:PL ./.:2,0:2:.:.:.:0,0,0 0/0:1,0:1:3:.:.:0,3,33 0/1:1,2:3:26:0|1:13813_T_G:81,0,26 GT:AD:DP:GQ:PL ./.:0,0:0:.:0,0,0 1/1:0,4:4:12:131,12,0 ./.:0,0:0:.:0,0,0 GT:AD:DP:GQ:PGT:PID:PL 1/1:0,2:2:6:1|1:14590_G_A:90,6,0 0/0:3,0:3:9:.:.:0,9,98 0/0:1,0:1:3:.:.:0,3,30 GT:AD:DP:GQ:PGT:PID:PL 1/1:0,2:2:6:1|1:14590_G_A:90,6,0 0/0:3,0:3:9:.:.:0,9,98 0/0:1,0:1:3:.:.:0,3,30 GT:AD:DP:GQ:PGT:PID:PL 1/1:0,2:2:6:1|1:14590_G_A:90,6,0 0/0:3,0:3:9:.:.:0,9,98 0/0:1,0:1:3:.:.:0,3,30 GT:AD:DP:GQ:PGT:PID:PL 1/1:0,2:2:6:1|1:14590_G_A:90,6,0 0/0:2,0:2:6:.:.:0,6,72 0/0:1,0:1:3:.:.:0,3,30 GT:AD:DP:GQ:PL 1/1:0,7:7:21:186,21,0 0/1:5,4:9:79:79,0,103 ./.:1,0:1:.:0,0,0
Ce n'est clairement pas ce que je veux!
Toute aide serait être reçu avec gratitude.
Ps. Je suis nouveau dans ce domaine, alors allez-y doucement.
3 Réponses :
La question semble claire, conserver les données avant :
dans chaque colonne (sauf dans la première). Mais ensuite, je suis un peu confus par la tentative de code, qui passe inutilement par des modèles spécifiques pour les colonnes.
Voici une interprétation simple de ce qui est décrit
use warnings; use strict; use feature 'say'; my $header_line = <>; # drop the first line while (<>) { # line by line from files given on cmdline, or STDIN chomp; my ($fmt_col, @cols) = split /\t/; # (but sample has spaces, not tabs) s/(.*?):.*/$1/ for @cols; # keep up to (first) : in each field say join "\t", $fmt_col, @cols; # print fields joined by tabs }
Le
code> lit toutes les lignes de tous les fichiers donnés sur la ligne de commande, ou STDIN
; alors soumettez un nom de fichier à traiter sur la ligne de commande lors de l'exécution de ce programme.
Notez que l'exemple publié n'a pas d'onglets, mais plutôt d'espaces; donc le code ci-dessus échouera si on le copie-colle. Changez les espaces en tabulations à tester, ou changez split / \ t /;
en split;
(afin d'utiliser sa valeur par défaut pour fractionner, qui correspond à n'importe quelle quantité d'espace ).
Tous les champs autres que le premier sont modifiés de sorte que seuls les caractères jusqu'au premier :
restent.
Ceci est fait en utilisant le fait que chaque l'élément traité dans la boucle foreach
("topicalizer") est aliasé sur l'élément actuellement traité. Ainsi, lorsque l'expression régulière s ///
change, l'élément correspondant de @cols
est changé. Si cela semble trop cuisant à digérer, veuillez l'écrire gentiment et lentement.
Si quelque chose d'autre est réellement nécessaire à faire, veuillez clarifier.
Merci. Quelques points à noter: 1) Je ne peux pas simplement utiliser les 3 caractères avant la première règle deux-points, car j'ai spécifiquement besoin de retirer les situations où il est 0/0
0/1
1/1
./.
uniquement (à l'exclusion des autres options) puis la partie suivante de mon code prend ces entrées et leur donne un nom selon que ils sont 0/0
etc. C'est pourquoi j'ai utilisé l'expression régulière pour extraire la notation spécifique et ensuite le code traduit cela en conséquence.
@lecb je vois, cela explique cela. J'ajouterai à cette réponse au fur et à mesure que j'y arriverai / (Alors que faire des champs qui n'ont pas ces modèles particuliers? Laissez-les simplement tels quels?)
Je pense avoir résolu le problème. J'ai supprimé foreach mes $ lines (@colsplit) {
et maintenant ça marche! Code ci-dessous:
#!/usr/bin/perl use strict; my $inputfile1 = $ARGV[0]; open (FILE1, $inputfile1) or die "Uh oh.. unable to find file $inputfile1"; ##Opens input file my @file1 = <FILE1>; #loads inputfile1 data into array close FILE1; my (@colsplit, @genotypes1, @genotypes2, @genotypes3, @joined); foreach my $line(@file1) { ## process each line, splitting columns and move onto next line @colsplit = split("\t", $line); push (@joined, $colsplit[0]); if ($colsplit[1] =~ m/(^0\/1)/ || $colsplit[1] =~ m/(^0\/0)/ || $colsplit[1]=~ m/(^1\/0)/ || $colsplit[1] =~ m/(^1\/1)/ || $colsplit[1] =~ m/(^.\/.)/) { push (@genotypes1, $1); } if ($colsplit[2] =~ m/(^0\/1)/ || $colsplit[2] =~ m/(^0\/0)/ || $colsplit[2] =~ m/(^1\/0)/ || $colsplit[2] =~ m/(^1\/1)/ || $colsplit[2] =~ m/(^.\/.)/) { push (@genotypes2, $1); } if ($colsplit[3] =~ m/(^0\/1)/ || $colsplit[3] =~ m/(^0\/0)/ || $colsplit[3] =~ m/(^1\/0)/ || $colsplit[3] =~ m/(^1\/1)/ || $colsplit[3] =~ m/(^.\/.)/) { push (@genotypes3, $1); } } my $i = 0; foreach my $line(@joined) { if ($line =~ m/GT/) { print "$line\t$genotypes1[$i]\t$genotypes2[$i]\t$genotypes3[$i]\n"; $i++; }}
Merci à @zdim en particulier de m'avoir montré une solution plus élégante. E
Une manière possible est d'utiliser le modificateur de continuation ( c code >)
de la correspondance d'expressions régulières. Cela provoque le démarrage d'une opération de correspondance à partir de l'endroit où la dernière a trouvé une correspondance. Comme ça, chaque champ d'une ligne peut être traité.
#!/usr/bin/perl use strict; use warnings; my $fh; $ARGV[0] && open($fh, $ARGV[0]) || die(); # discard first line; <$fh>; foreach my $line (<$fh>) { chomp($line); # capture every non tab character from the beginning of the line; # globally match the pattern repeatedly in the string (g modifier) # keep the current position during repeated matching (c modifier); $line =~ m/^([^\t]*)/gc; print("$1"); # capture every non colon character after a tab character; # globally match the pattern repeatedly in the string (g modifier); # keep the current position during repeated matching (c modifier); while ($line =~ m/\t([^:]*)/gc) { print("\t$1"); } print("\n"); }
Utilisez PyVCF ou une autre bibliothèque pour l'analyse VCF plutôt que de déployer la vôtre.