Discussion :

[karaudrey88]

Bonjour
Au cours de mon stage j'ai réalisé un pipeline permettant d'analyser des variations génomiques. Pour chaque région génomique et chaque échantillon, des fichiers txt et des graphiques sont générés permettant cette analyse. Plusieurs analyses sont lancées en parallèle. pour faciliter l'analyse des résultats, il faudrait que je fasse un script perl permettant de faire une synthèse des résultats par région génomique étudiée. Malheureusement je ne sais pas comment m'y prendre.Quelqu'un pourrait il m'aiguiller?
Merci

[Gardyen]

hum que veux-tu dire par faire une synthèse ? tu peux préciser un peu plus ?
puisque tes résultats sont liés à une position, tu peux l'utiliser pour retrouver les infos concernées non ?

[karaudrey88]

Par exemple pour chaque échantillon et chaque région génomique donnés, des fichiers txt contenant la région génomique, les positions de la variation structurale et ses caractéristiques sont générés en parallèle.
Ce que j’appelle faire une synthèse par région génomique, c'est regrouper les variations structurales détectées dans cette région dans un même fichier quelque soit les échantillons.
Et je ne peut pas faire un copier coller avec les commandes perl ou linux car comme les échantillons sont traités en parallèle je ne peut opier plusieurs choses en même temps dans le même fichier. De plus il faut que rajoute le nom de l’échantillon dans les caractéristiques (le nom de l'échantillon est dans le nom du fichier txt).

[Gardyen]

de ce que j'ai compris, chaque fichier de résultat est spécifique d'un échantillon et d'une position ? par contre il y a plusieurs types de fichiers pour un même échantillon et une même position ?

si tu peux savoir quand un job se termine, alors tu peux ensuite transférer les infos vers un support centralisé, comme un fichier spécifique pour une position.
plusieurs jobs en parallèle, plusieurs fichiers de résultat, mais 1 seul "contrôleur" qui exécutera une simple copie une fois le job terminé.

[karaudrey88]

Mon pipeline se réalise en plusieurs étapes, et pour cela je contrôle le déroulement des jobs. Une fois que tous les jobs d'une étape sont terminés je lance la seconde étape. Donc je peux savoir quand mes jobs se terminent. En revanche un simple copier coller des résultats dans un même fichier ne fonctionne pas à l'intérieur du pipeline. Il fonctionne seulement quand il est réalisé "a la main" après que le pipeline ait fini de tourner. De plus, une simple copie ne me suffit pas, j'ai besoin de récupérer le nom du fichier a chaque copie d'une variation structurale dans le fichier centrale afin de savoir a quel échantillon elle appartient.

[Gardyen]

comment démarre-tu le pipeline ?
y a-t-il un contrôleur qui gère les jobs ?
peux-tu montrer des exemples de fichiers de résultats et de fichiers finaux ?

je n'arrive pas à me faire une idée de ce que tu veux au final

[Lolo78]

Tu pourrais donner une idée de tes données en entrée (et de ce que tu veux en sortie)?

Là, pour l'instant, on n'a pas une idée claire de ce que tu veux.

[karaudrey88]

Mon pipeline commence par la définition de toutes les options a mettre en ligne de commandes et gestion des erreurs allant avec. Ensuite commencent les étapes d'analyse.
Je gère mes jobs de la façon suivante :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
if ($submit eq "yes") {
	print STDERR "Launching individual bam file jobs";
}
my (%pileup_jobs_id,$ref_bam_pileup_job,@prefix_list);
foreach my $bam (@bamfiles) {
	#create a directory for each bam
	my @nn=split(/\./,$bam); 
	my $prefix=$nn[0];
	my $command = "mkdir $prefix";
	qx{$command};
 
	push(@prefix_list,$prefix);
	#create qsub for each chr
 
	foreach my $chr (@chrs) {
		$command = "mkdir $prefix/$chr";
		qx{$command};
		#generate qsub file and launch it
		my $qsuboutfile = "$bam\_$chr\_bam2cnv.qsub";
		open (QSUB,">$prefix/$chr/$qsuboutfile")or die "\nError : $qsuboutfile file could not be created: $!\n";
		print QSUB qq(#!/bin/bash
#\$ -S /bin/bash
#\$ -V
#\$ -M $users{$user}
#\$ -m e
#\$ -o $qsuboutfile.out
#\$ -e $qsuboutfile.err
#\$ -cwd
#\$ -q long.q
#\$ -pe thread 5
 
# command(s):
hostname >&2
 
XXX : commandes a effectuer en bash
);
		close QSUB;
 
		#launch qsub job if asked
		my $jobid;
		if ($submit eq "yes") {
			my $dependencies = "";
			if (defined $hold_jid) {
				$dependencies = "-hold_jid $hold_jid";
			}
			my $instruction = "cd $prefix/$chr; qsub $dependencies $qsuboutfile | awk '{print \$3}' ";
			$jobid=qx{$instruction};
			chomp $jobid;
			if ($bam ne $refbam) {
				$pileup_jobs_id{$bam} = $jobid;
			} else {
				$ref_bam_pileup_job =  $jobid;
			}
			print STDERR ".";
		} else {
			$pileup_jobs_id{$bam} = "";
		}
	}
}

Sinon voici l'exemple d'un fichier de sortie généré par le pipeline

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
XXXX	        12900	13076	0.828	      1.45e-07	2.290
XXXX 	14610	14760	1.412       1.44e-10	4.10
XXXX	        15890	16002	1.049       7.10e-12	2.85

Voici trois variations structurales provenant d'un échantillon Y sur une région génomique XXXX.

Et j'aimerai obtenir un fichier du type :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
XXXX	      Y1       12900	13076	0.828	      1.45e-07	2.290
XXXX	      Y8       12900	13076	0.828	      1.45e-07	2.290
XXXX       Y2       12950        13076       0.589       3.02e-03      1.960
XXXX       Y5       12950        13076       0.589       3.02e-03      1.960

En fait je veux regrouper toutes les variations structurales provenant de différents échantillons mais de la même région génomique dans un même fichier en ajoutant la provenance de la variation détectée..

[Lolo78]

Si je comprends bien, tu veux regrouper toutes les lignes commençant par le même mot.

Est-ce ça?

[karaudrey88]

C'est une partie de ce que j'aimerai bien faire effectivement.

[Lolo78]

Je ne comprends toujours pas.

Je regarde tes données en sortie du pipeline et celles que tu voudrais avoir. Et je suppose que les données en sortie du pipeline doivent être triturées pour obtenir ce que tu voudrais avoir au final. Déjà, est-ce que c'est ça?

Entre les deux séries de données, le premier enregistrement est reconnaissable, mais je ne sais pas d'où vient le "Y1" que tu as ajouté, il ne semble venir de nulle part. Quant aux trois lignes suivantes, je ne vois pas du tout d'où elles viennent.

Je précise que je te parle en termes de données (en termes informatiques, donc), parce que je n'y connais rien en génétique: la dernière fois que j'ai fait de la biologie en fac, la génétique en était à ses tous premiers balbutiements, tout cela a énormément changé. Et comme je n'ai jamais eu l'occasion de travailler dans ce domaine, j'en suis resté à ce niveau proche de 0 (j'ai même probablement oublié le peu que je savais à l'époque).

Bref il faudrait quelles sont données en entrée et comment tu les veux en sortie. Pour l'instant, je ne suis même pas sûr que ce soit clair dans ton esprit.

[karaudrey88]

Effectivement, les données doivent être un peu "triturées". Les Y sont les différents échantillons d'ou proviennent les variations structurales.
Les données d'entrées (variations) sont classées en fonction de l'échantillon et de la région génomique étudier. La région génomique est présente dans les données ce qui correspond a XXXX. En revanche les Y (échantillons) sont présents dans le nom du fichier. Donc pour avoir le nom de l'échantillon il faut récupérer le nom du fichier d'ou les variations sont issues.

[Lolo78]

OK, j'ai compris maintenant.

Voici une solution possible: tu lis toutes tes données en entrée en prenant soin de les enrichir avec le nom de l'échantillon, et stokes toutes ces lignes dans un hash dont la clé est la région génomique (xxxx).

A la fin, quand tu as lus tous tes fichiers, tu parcours la hash pour l'imprimer dans un fichier.

[Lolo78]

Quelque chose comme cela (non testé):

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
my %hash;
$nom_echantillon = "Y1";
while (my $line = <FILI_IN>) {
     chomp $line;
     my ($region, @reste_de_la_ligne) = split /\s/, $line;
     my $newline = join "\t", $region, $nom_echantillon, @reste_de_la_ligne, "\n";
     if (defined $hash{$region}) {
          $hash{$region} .= $newline;
     }
     else {
          $hash{$region} = $newline;
     }
}

A la fin de ce processus, quand tu as lu tous tes fichiers, $hash{XXXX} va contenir un enregistrement multi-ligne regroupant toutes les lignes relatives à XXXX.

Il suffit alors de balayer le hash et d'imprimer les valeurs:

Par exemple, si l'ordre dans lequel les régions apparaissent n'a pas d'importance:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

print foreach values %hash;

Si tu veux les régions par ordre alphabétique (ou plutôt "ASCIIbétique"):

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
foreach my $region (sort keys %hash) {
     print OUTFILE $hash{$region};
}

[karaudrey88]

Merci
Effectivement c'est une idée qui ne m'était pas venu à l'esprit, je vas la tester.
En tout cas merci de votre patience et de votre aide.