CpG o/e script

**wanaga3166** · 31/01/2014, 14h39

Bonjour à la communauté de bioinformaticien,

Je dois calculer le ratio CpG o/e de plusieurs espèces, à partir d'un script qu'on m'a donné et que je dois modifier.

Du coup, je m'exécute et je fais donc une rapide initiation au Perl (en autonome). En effet, j'ai rajouté en bas de ce script la formule pour calculer le ratio CpG o/e (ce qui paraissait le plus simple), à partir de la même formule que celle décrite dans cet article (http://www.biomedcentral.com/1471-2164/11/483).

Comme toujours, il y a un problème. Or, quand je lance le script, on me répond qu'il est impossible de le lancer et à cause de: division by zéro.

J'ai cherché comment résoudre ce problème. Mais je suis bloqué (je rappelle que je début dans la programmation en Perl :-) ).

Voici le code:

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
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
 
 
#!/usr/bin/perl -w
 
####################################
#
# getfastanamesandlength.pl -f fasta-file -o output_file
#
# reads concatanated fasta files, writes name, length, CpGsand GpCs into TAB file #for each fasta sequence
#####################################
 
use diagnostics;
use strict;
use Carp;
use FileHandle;
use File::Path;
use File::Basename;
use Getopt::Std;
 
my @methylnam         = ();
my @methylseq         = ();
my @methylnam1         = ();
my @methylseq1         = ();
my @methylnam1_sorted = ();
my @temp = ();
my $minimum_length    = 0; #minimal acceptable sequence length
my $acceptable_proc_n = 100;   #minimal acceptable percentage of non ATGC
my $seq_length        = 0;
my $proc_n            = 0;
my $FILE              = FileHandle->new();
my $seq = "";
my $num = 0;
my $i;
my $duplicates_nr = 0;
my $outfile = "";
my $j;
my $seqnumber = 0;
my $seqname = "";
my $CG_number = 0;
my $C_number = 0;
my $G_number = 0;
my $GC_number = 0;
 
use vars qw( $opt_f $opt_o);
# Command line parsing 
#
getopts('f:o:');
#
# Read input file and split into fasta sequences on the fly
#
#replace Mac and PC CRLF by UNIX LF
system ("perl -pi -e 's/\r/\n/g' $opt_f");
open($FILE,$opt_f) || croak(sprintf("Cannot open file \"%s\"",$opt_f));
while ( <$FILE> ) {
    chomp;
    s/\ /_/g;
    if ( /^>(\S+)/) 
         {
         push(@methylnam,$1);
         push(@methylseq,"");
         } 
       else 
         {
	     $methylseq[-1] .= $_;
         }
       }
close($FILE) || croak(sprintf("Cannot close open file \"%s\"",$opt_f));
 
$outfile = $opt_o;
#print header
open($FILE,">".$outfile) || croak(sprintf("Cannot open file \"%s\"",$outfile));
print $FILE "#name\tlength\tCpGs\tGpCs\tCs\tGs\tCpG o\/e\n";
close ($FILE);
 
open($FILE,">>".$outfile) || croak(sprintf("Cannot open file \"%s\"",$outfile));
 
for $i (0 .. $#methylnam)
         {
         @temp = ();
         @temp=split('\|',$methylnam[$i]);
         $seqname=$temp[1];         
         print $FILE (($methylnam[$i])."\t".(length($methylseq[$i]))."\t".($methylseq[$i]=~ s/CG/CG/gi)."\t".($methylseq[$i]=~ s/GC/GC/gi)."\t".($methylseq[$i]=~ s/C/C/gi)."\t".($methylseq[$i]=~ s/G/G/gi)."\t".((($methylseq[$i]=~ s/CG/CG/gi)/(($methylseq[$i]=~ s/C/C/gi)*($methylseq[$i]=~ s/G/G/gi)))*(((length($methylseq[$i]))^2)/(length($methylseq[$i]))-1))."\n");         
         }
close ($FILE);

Merci de votre aide et de vos conseil pour le futur :-)

**wanaga3166** · 03/02/2014, 14h53

Pour éviter ce illegal division by zero, j'ai voulu rajoute des conditions dans la boucles for ci dessous:

Code :

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

1
2
3
4
5
6
7
8
 
for $i (0 .. $#methylnam)
         {
         @temp = ();
         @temp=split('\|',$methylnam[$i]);
         $seqname=$temp[1];         
         print $FILE (($methylnam[$i])."\t".(length($methylseq[$i]))."\t".($methylseq[$i]=~ s/CG/CG/gi)."\t".($methylseq[$i]=~ s/GC/GC/gi)."\t".($methylseq[$i]=~ s/C/C/gi)."\t".($methylseq[$i]=~ s/G/G/gi)."\t".((($methylseq[$i]=~ s/CG/CG/gi)/(($methylseq[$i]=~ s/C/C/gi)*($methylseq[$i]=~ s/G/G/gi)))*(((length($methylseq[$i]))^2)/(length($methylseq[$i]))-1))."\n");         
         }

Les conditions seraient que C*G soient >0. Mais je n'y arrive pas. Des conseils s'il vous plaît. Merci beaucoup.

**wanaga3166** · 03/02/2014, 15h18

ou utiliser

Code :

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

foo

?

**Gardyen** · 03/02/2014, 15h40

peux-tu donner quelques lignes du fichier à traiter ? ainsi que l'erreur exacte ?

tu auras une division par zéro s'il n'y a pas de C ou de G dans ta séquence, ou si ton $methylseq[$i] est vide.

As tu vérifié que tes fichiers intermédiaires ont le bon format ?

**wanaga3166** · 03/02/2014, 16h10

Envoyé par Gardyen

peux-tu donner quelques lignes du fichier à traiter ? ainsi que l'erreur exacte ?

Voici un example d'EST de la fourmi de feu (Solenopsis invicta), sur laquelle je fais les tests. En realité, ce scripts devra tourner sur plusieurs milliers d'EST récuperés sur ESTdb et de plusieurs milliers d'espèces.

L'erreur exacte est:

Illegal division by zero at ./scriptcpg.pl line 80 (#2)
(F) You tried to divide a number by 0. Either something was wrong in
your logic, or you need to put a conditional in to guard against
meaningless input.

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
 
>gi|165907545|gb|FD483965.1|FD483965 FAF01116 Fire Ant cDNA library Solenopsis invicta cDNA 5', mRNA sequence
CCCACTTTTGAGAAGTCAAACATATTAATCCTACATTGATTAACCTGGATCATATTCACT
AGAGCAATTTTTTTCTAATTTTCCCGATTCTCTTTCTCCTTTTCTCTTGAGTGTTAAATG
AATATACCAGAATTTCTTTCAGACTTTAAGAATCCAAAAATTCTTTTTGAAAAAACATAT
CTTTTTACAGTGACTTAGGTAATTTAACCAAGCAGTGAAAAATAAGAGTTTACGGAGGGT
TAAATGTCTTTATTATAACAGGAAAATTCTGCTAAAATATACTTTATAAATTTAATTTAT
TTGTAGCTCTACTTTATGATTGCATTTTAAATCTAATCATATTGGCATAGTTTATAGAAC
TAAAAAATGTTTTATAATTGATATATGGAAATATTATTGAACATCAATGGCACTTGTTTT
TTTATGTTATATCCCAAGATTATAAAATGCTAGCAACTTATCAACGGTATGCCAGCAAAC
TCGATCATGACTGTAGAACATATTCTGTCTAATGTTTTTTTTTTAAGAGCCACTTATCAA
TTCCAGTTAATCTAATTAATTGATTAGNCTATTAGAATTGACATAGTTGTCATTTTACAT
AACAGACAGATGCAGTTGGTCATTCCAATGGAAAATTCAATTAGATTAACCAGAGTTGAT
GAAATGGCCCTAAGCTTGGTTTCTATTCATTACAATTCAACTATAATATAAGAACAGAAA
TG

tu auras une division par zéro s'il n'y a pas de C ou de G dans ta séquence, ou si ton $methylseq[$i] est vide.

Il y a bien des C et des G dans mes séquences. J'ai voulu faire une boucle for pour dire que C*G doit être supérieur à 0. Quant à $methylseq[$i], je ne pense pas qu'il soit vide. En effet, lorsque je modifie la formule (CpG o/e= CG/(C*G) * (l²/(l-1)), en échangeant le / par *, le script fonctionne normalement et il n'y a pas d'illigal division by zero.

As tu vérifié que tes fichiers intermédiaires ont le bon format ?

Les fichiers ont le bon format normalement (fasta). Par contre après trimming, les fichiers sont automatiquement renommés comme dans l'exemple suivant:

fire_ant.fasta_out.fst

Mais, c'est toujours du fasta.

**Gardyen** · 03/02/2014, 17h36

intéressant, j'ai copié ton script et ton fichier d'entrée et je n'ai eu aucun problème pour lancer le script

as-tu une erreur si tu testes uniquement avec ta séquence ?

**wanaga3166** · 04/02/2014, 09h57

Envoyé par Gardyen

intéressant, j'ai copié ton script et ton fichier d'entrée et je n'ai eu aucun problème pour lancer le script

Envoyé par Gardyen

as-tu une erreur si tu testes uniquement avec ta séquence ?

En réalité, lorsque je teste, il y a un fichier Fire ant avec plusieurs séquences fasta comme celle que j'ai donné à titre d'exemple. Ensuite, lorsque je teste uniquement avec cette sequence comme toi, cela marche.

Mais comment faire pour plusieurs séquences ? Où est le problème ?

**Gardyen** · 04/02/2014, 10h32

j'ai testé en dupliquant la séquence et ça fonctionne correctement.
je pense donc que tu devrais vérifier tes séquences, en affichant le nombre de C, de G et sa longueur avant la ligne lançant l'erreur, quelque chose comme

Code :

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

print $methylnam[$i]." ".(length($methylseq[$i]))." ".($methylseq[$i]=~ s/CG/CG/gi)." ".($methylseq[$i]=~ s/GC/GC/gi)." ".($methylseq[$i]=~ s/C/C/gi)." ".($methylseq[$i]=~ s/G/G/gi)." "

CpG o/e script

Bioinformatique Perl

Discussions similaires

Partager

Partager