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:
Merci de votre aide et de vos conseil pour le futur :-)
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);
Partager