annotate FunAnn/annotate.pl @ 0:1886293c2243 draft

Uploaded
author elixir-it
date Thu, 23 Jul 2020 12:49:09 +0000
parents
children 9b8afa0b6da4
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
1 $fss=13468;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
2
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
3 unless (-e "CorGAT")
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
4
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
5 {
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
6 system("wget -i https://raw.githubusercontent.com/matteo14c/CorGAT_galaxy/dev/ann.txt");
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
7 system("gzip -d GCF_009858895.2_ASM985889v3_genomic.fna.gz");
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
8 }
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
9
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
10
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
11
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
12
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
13 $gen_code="genetic_code";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
14 die("need genetic code file in the current folder\n") unless -e "genetic_code";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
15 open(IN,$gen_code);
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
16 while(<IN>)
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
17 {
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
18 ($triplet,$oneL)=(split());
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
19 $code{$triplet}=$oneL;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
20 }
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
21
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
22 $genome="GCF_009858895.2_ASM985889v3_genomic.fna";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
23 die("need reference genome file in the current folder\n") unless -e "GCF_009858895.2_ASM985889v3_genomic.fna";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
24 open(IN,$genome);
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
25 while(<IN>)
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
26 {
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
27 next if $_=~/^>/;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
28 chomp;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
29 $seq.=$_;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
30 }
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
31 $table="annot_table.pl";#"simple_annot_mirror";#"annot_table.pl";#sl5 29728 29768 reg
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
32 die("need detailed annotation file for SARS-CoV-2 in the current folder") unless -e "annot_table.pl";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
33 open(IN,$table);
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
34
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
35
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
36 while(<IN>)
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
37 {
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
38 chomp();
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
39 ($gene,$b1,$b2,$e,$annot,$notes)=(split(/\t/));
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
40 $annot{$b1}{$b2}=[$gene,$e,$annot,$notes];
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
41 #print "$gene $b1 $b2 $e $annot $notes\n";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
42 $len=$b2-$b1+1;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
43 if ($gene ne "nsp12" && $gene ne "orf1ab")
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
44 {
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
45 $seqgene=substr($seq,$b1-1,$len+3); #un codone inizio in più e un codone fine in più
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
46 }else{
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
47 $len1=$fss-$b1+1;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
48 $part1=substr($seq,$b1-1,$len1);
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
49 $len2=$b2-$fss+1;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
50 $part2=substr($seq,$fss-1,$len2+3);
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
51 $seqgene="$part1$part2";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
52
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
53 }
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
54 $annot_seq{$gene}=$seqgene;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
55 $Lgenes{$gene}=length($seqgene);
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
56 ($NSTOP_G,$Tseq,$pos_Stop_R)=translate($seqgene,\%code);
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
57 @seq_res=split('',$Tseq);
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
58 for ($i=0;$i<=$#seq_res;$i++)
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
59 {
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
60 $pos=$i+1;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
61 $res=$seq_res[$i];
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
62 }
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
63 }
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
64 %AF_data=%{read_simple_table("af_data_new.csv")};
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
65 %MFE_data=%{read_simple_table("MFE_annot.csv")};
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
66 %epi_data=%{read_epitopes("epitopes_annot.csv")};
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
67 %hyphy_data=%{read_hyphy("hyphy_novel.csv")};
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
68
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
69
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
70
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
71 $var_File=shift;#""cl7.csv";#"phenetic_indels_sars_cov2.csv";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
72 $out_File=shift;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
73 open(OUT,">$out_File");
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
74 die("no input\n") unless -e $var_File;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
75 open(IN,$var_File);
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
76 $header=<IN>;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
77 @header=(split(/\s+/,$header));
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
78 print OUT "POS\tREF\tALT\tannot\tAF\tEpitopes\tHyphy\tMFE\n";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
79 while(<IN>)
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
80 {
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
81 ($change,@pos)=(split());
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
82 next unless $change=~/\|/;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
83 ($pos,$allele)=(split(/_/,$change))[0,1];
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
84 ($ref,$alt)=(split(/\|/,$allele))[0,1];
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
85 $AF=$AF_data{"$pos$ref$alt"} ? $AF_data{"$pos$ref$alt"} : 0;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
86 next if $alt=~/N/;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
87 $annot_string="";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
88 $contained=0;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
89 $epitope_string="";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
90 $hyphy_string="";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
91 $MFE_string= $MFE_data{"$pos$ref$alt"} ? $MFE_data{"$pos$ref$alt"} : "NA";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
92 #next if $ref eq"." || $alt eq ".";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
93 foreach $b1 (sort{$a<=>$b} keys %annot)
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
94 {
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
95 foreach $b2 (sort{$a<=>$b} keys %{$annot{$b1}})
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
96 {
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
97 if ($pos<=$b2 && $pos>=$b1)
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
98 {
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
99 $contained=1;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
100 $type=$annot{$b1}{$b2}[1];
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
101 $namegene=$annot{$b1}{$b2}[0];
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
102 if ($type eq "cds")
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
103 {
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
104 #print "cds";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
105 @res=annot_CDS($pos,$ref,$alt,$namegene);
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
106 $annot_string.=$res[0];
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
107 if ($namegene ne "orf1ab")
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
108 {
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
109 $hyphy_string.=$res[1];
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
110 $epitope_string.=$res[2];
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
111 }
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
112 }else{
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
113 $rel_pos=$pos-$b1+1;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
114 $annot_string.="$namegene:nc.$ref$rel_pos$alt,NA,NA;";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
115 $epitope_string="NA" if $epitope_string eq "";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
116 $hyphy_string="NA" if $hyphy_string eq "";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
117 }
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
118 }
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
119 }
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
120 }
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
121 $epitope_string=~s/\s+/;/g;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
122 #$epitope_string="EpiT:$epitope_string" unless $epitope_string eq "NA";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
123 print OUT "$pos\t$ref\t$alt\t$annot_string\t$AF\t$epitope_string\t$hyphy_string\t$MFE_string\n" #if $contained==1;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
124 }
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
125
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
126 sub translate
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
127 {
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
128 @orig_seq=split('',$_[0]);
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
129 %gen_code=%{$_[1]};
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
130 $type="";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
131 $NSTOP=0;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
132 $Tseq="";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
133 $pos_Stop=0;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
134 for ($i=0;$i<=$#orig_seq;$i+=3)
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
135 {
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
136 $AA=join('',@orig_seq[$i..$i+2]);
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
137 $res=$gen_code{$AA};
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
138 $pos_Stop=$i if $pos_Stop ==0 && $res eq "*";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
139 $NSTOP++ if $res eq "*";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
140 $Tseq.=$res;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
141 }
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
142 #print "$seq\n";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
143 return($NSTOP,$Tseq,$pos_Stop);
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
144 }
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
145
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
146 sub annot_CDS
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
147 {
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
148 $pos=$_[0];
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
149 $ref=$_[1];
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
150 $alt=$_[2];
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
151 $namegene=$_[3];
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
152 my $hyphy_string="NA";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
153 my $epitopes_string="NA";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
154 #print "$pos $ref $alt $namegene\n";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
155 $pos_inG=$pos-$b1+1;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
156 $pos_inG++ if $pos >$fss && ($annot{$b1}{$b2}[0] eq "nsp12" || $annot{$b1}{$b2}[0] eq "orf1ab");
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
157 $mod=$pos_inG%3;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
158 $rel_pos=int($pos_inG/3);
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
159 $rel_pos++ if $mod !=0;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
160
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
161 if (length($ref)==1 && $ref ne "."&& $alt ne ".")
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
162 {
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
163
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
164 if ($mod ==1)
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
165 {
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
166 $triplet=substr($seq,$pos-1,3);
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
167 @Bs=split('',$triplet);
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
168 die("1\n $triplet b:$Bs[0] r:$ref") unless ($Bs[0] eq $ref);
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
169 $Bs[0]=$alt;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
170 }elsif ($mod==2){
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
171 $triplet=substr($seq,$pos-2,3);
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
172 @Bs=split('',$triplet);
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
173 die("2\n $triplet b:$Bs[1] r:$ref")unless ($Bs[1] eq $ref);
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
174 $Bs[1]=$alt;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
175 }elsif ($mod==0){
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
176 $triplet=substr($seq,$pos-3,3);
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
177 @Bs=split('',$triplet);
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
178 die("3\n $triplet b:$Bs[2] r:$ref")unless ($Bs[2] eq $ref);
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
179 $Bs[2]=$alt;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
180 }
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
181 #print "$pos_inG $relpos $mod @Bs\n";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
182 $Atriplet=join("",@Bs);
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
183 $A1=$code{$triplet};
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
184 $A2=$code{$Atriplet};
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
185 $eff=$A1 eq $A2 ? "synonymous" : "missense";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
186 $eff="start_lost" if $eff eq "missense" && $rel_pos ==1;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
187 $eff="stopGain" if $A2 eq "*" && $A1 ne "*"; #stopG
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
188 $eff="stopLoss" if $A1 eq "*" && $A2 ne "*"; #stopL
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
189 $eff="S" if $A2 eq "*" && $A1 eq "*";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
190 $hyphy_string=$hyphy_data{"$namegene$rel_pos"} if $hyphy_data{"$namegene$rel_pos"};
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
191 $epitopes_string=join(" ",@{$epi_data{"$namegene$rel_pos$A1"}}) if $epi_data{"$namegene$rel_pos$A1"};
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
192 #return("$namegene:c.$pos_inG$ref>$alt,p.$A1$rel_pos$A2,$mod,$eff;",$hyphy_string,$epitopes_string);
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
193 return("$namegene:c.$pos_inG$ref>$alt,p.$A1$rel_pos$A2,$eff;",$hyphy_string,$epitopes_string);
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
194 }else{
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
195 $eff="";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
196 $gene=$annot_seq{$namegene};
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
197 $lgene=length($gene);
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
198 $upstream=substr($gene,0,$pos_inG-1);
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
199 $change=substr($gene,$pos_inG,length($ref));
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
200 $downstream=substr($gene,$pos_inG+length($ref)-1);
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
201 $len=length($alt);
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
202 unless ($ref=~/\./)
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
203 {
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
204 $modSeq="$upstream$alt$downstream";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
205 }else{
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
206 $modSeq="$upstream$alt$change$downstream";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
207 }
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
208 $modSeq=~s/\.//g;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
209
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
210 ($NSTOP_G,$Tseq_R,$pos_Stop_R)=translate($gene,\%code);
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
211 ($NSTOP_R,$Tseq_alt,$pos_Stop_T)=translate($modSeq,\%code);
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
212 $CDS_annot_string=".";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
213 #print "$namegene,aa: sr:$NSTOP_G sa:$NSTOP_R psr:$pos_Stop_R psa:$pos_Stop_T\n";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
214 if ($NSTOP_R>$NSTOP_G && ($ref=~/\./ || $alt=~/\./))
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
215 {
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
216
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
217 $eff="frameshift";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
218 $eff.="Ins" if $ref=~/\./;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
219 $eff.="Del" if $alt=~/\./;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
220 $truncation=$pos_Stop_T/3;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
221 $ref_Seq=substr($Tseq_R,$truncation-1,1);
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
222 return("$namegene:c.$pos_inG$ref>$alt,p.$ref_Seq$truncation*,$eff;",$hyphy_string,$epitopes_string);
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
223 #return("$namegene:$pos_inG,$rel_pos,$mod,Frameshift,$eff;",$hyphy_string,$epitopes_string);#\t$CDS_annot_string");
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
224 }elsif($NSTOP_R==$NSTOP_G && $pos_Stop_T<($pos_Stop_R-$len) ){
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
225 # print "2\n";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
226 $truncation=$pos_Stop_T/3;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
227 $ref_Seq=substr($Tseq_R,$truncation-1,1);
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
228 $eff="Truncating";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
229 $eff.="Ins" if $ref=~/\./;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
230 $eff.="Del" if $alt=~/\./;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
231 return("$namegene:c.$pos_inG$ref>$alt,p.$ref_Seq$truncation*,$eff;",$hyphy_string,$epitopes_string);
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
232 #return("$namegene:$pos_inG,$rel_pos,$mod,Trunc:$truncation,$eff;",$hyphy_string,$epitopes_string);#\t$CDS_annot_string=");
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
233 }else{
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
234 %used_epitopes=();
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
235 $hyphy_string="";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
236 $epitopes_string="";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
237 $pre="";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
238 $Cref=$ref;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
239 $Calt=$alt;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
240 if ($mod==0)
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
241 {
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
242 $pre=substr($gene,$pos_inG-3,2);
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
243 }elsif($mod==2){
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
244 $pre=substr($gene,$pos_inG-2,1);
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
245 }
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
246 $Cref="$pre$Cref";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
247 $Calt="$pre$Calt";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
248 $post=$pos_inG+length($ref)-1;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
249 while (length($Cref) %3!=0)
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
250 {
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
251 $Cref.=substr($gene,$post,1);
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
252 $Calt.=substr($gene,$post,1);
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
253 $post++;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
254 if ($post>length($gene))
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
255 {
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
256 $local_gene=substr($gene,0,length($gene)-3);
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
257 $eff="TruncatingDel";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
258 #print "$post\n";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
259 $copy_pos=$pos_inG;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
260 #print "$pos_inG\n";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
261 #print length($gene)."\n";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
262 #die("muoro");
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
263 $Loc_Ref=substr($local_gene,$copy_pos);
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
264 $Talt="-";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
265 while(length($Loc_Ref)%3!=0)
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
266 {
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
267 $copy_pos--;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
268 $Loc_Ref=substr($local_gene,$copy_pos);
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
269 }
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
270 $rel_pos=int($copy_pos/3);
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
271 $rel_pos++ if $mod !=0;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
272 $Tref=(translate($Loc_Ref,\%code))[1];
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
273 return("$namegene:c.$pos_inG$ref>$alt,p.$Tref$rel_pos$Talt,$eff;",$hyphy_string,$epitopes_string);
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
274 #return("$namegene:$pos_inG,$rel_pos,$mod,$Tref->$Talt,$eff;",$hyphy_string,$epitopes_string)
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
275 }
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
276 }
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
277 if ($alt=~/\./ || $ref=~/\./)
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
278 {
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
279 $Tref="-";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
280 $Talt="-";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
281 $eff="inframe";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
282 if ($ref=~/\./)
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
283 {
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
284 $eff.="Ins";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
285 $Talt=(translate($Calt,\%code))[1];
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
286 }
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
287 if ($alt=~/\./)
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
288 {
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
289 $eff.="Del";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
290 $Tref=(translate($Cref,\%code))[1];
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
291 }
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
292 if ($eff=~/Del/)
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
293 {
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
294 @Tref=split('',$Ttref);
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
295 for ($i=0;$i<=$#Tref;$i++)
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
296 {
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
297 $cur_res=$Tref[$i];
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
298 $cur_pos=$rel_pos+$i;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
299 if ($epi_data{"$namegene$cur_res$cur_pos"})
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
300 {
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
301 $used_epitopes{join(' ',@{$epi_data{"$namegene$cur_res$cur_pos"}})}=1;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
302 }
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
303 if ($hyphy_data{"$namegene$curpos"})
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
304 {
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
305 $hyphy_string.=$hyphy_data{"$namegene$curpos"} . ";";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
306 }
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
307 }
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
308 $Nkeys=keys %used_epitopes;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
309 $epitopes_string.=join(";",keys %used_epitopes) if $Nkeys>=1;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
310 }
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
311 $hyphy_string="NA" if $hyphy_string eq "";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
312 $epitopes_string="NA" if $epitopes_string eq "";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
313 #print "$Talt\n";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
314 return("$namegene:c.$pos_inG$ref>$alt,p.$Tref$rel_pos$Talt,$eff;",$hyphy_string,$epitopes_string);
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
315 #return("$namegene:$pos_inG,$rel_pos,$mod,$Tref->$Talt,$eff;",$hyphy_string,$epitopes_string);#\t$CDS_annot_string");
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
316 }else{
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
317 $Tref=(translate($Cref,\%code))[1];
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
318 $Talt=(translate($Calt,\%code))[1];
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
319 $eff="synonymous" if $Tref eq $Talt;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
320 $eff="missense" if $Tref ne $Talt;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
321 $eff="stopGain" if $Talt=~/\*/;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
322 $eff="stopLoss" if $Tref=~/\*/ && !$Talt=~/\*/;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
323 @Tref=split('',$Ttref);
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
324 for ($i=0;$i<=$#Tref;$i++)
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
325 {
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
326 $cur_res=$Tref[$i];
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
327 $cur_pos=$rel_pos+$i;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
328 if ($epi_data{"$namegene$cur_res$cur_pos"})
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
329 {
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
330 $used_epitopes{join(' ',@{$epi_data{"$namegene$cur_res$cur_pos"}})}=1;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
331 }
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
332 if ($hyphy_data{"$namegene$curpos"})
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
333 {
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
334 $hyphy_string.=$hyphy_data{"$namegene$curpos"} . ";";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
335 }
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
336 }
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
337 $Nkeys=keys %used_epitopes;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
338 $epitopes_string.=join(";",keys %used_epitopes) if $Nkeys>=1;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
339 $hyphy_string="NA" if $hyphy_string eq "";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
340 $epitopes_string="NA" if $epitopes_string eq "";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
341 return("$namegene:c.$pos_inG$ref>$alt,p.$Tref$rel_pos$Talt,$eff;",$hyphy_string,$epitopes_string);
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
342 #return("$namegene:$pos_inG,$rel_pos,$mod,$Tref->$Talt,$eff;",$hyphy_string,$epitopes_string);#\t$CDS_annot_string");
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
343 }
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
344 }
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
345 }
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
346 }
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
347
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
348
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
349 sub read_simple_table
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
350 {
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
351 $file=$_[0];
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
352 die ("$file does not exist") unless -e $file;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
353 my %data=();
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
354 open(IN,$file);
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
355 while(<IN>)
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
356 {
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
357 ($pos,$ref,$alt,$annot)=(split());
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
358 if ($annot=~/;/)
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
359 {
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
360 # $annot=(split(/\;/,$annot))[0];
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
361 }
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
362 $data{"$pos$ref$alt"}=$annot;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
363 }
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
364 return \%data;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
365 }
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
366
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
367 sub read_epitopes
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
368 {
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
369 $file=$_[0];
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
370 die("$file does not exist") unless -e $file;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
371 my %data=();
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
372 open(IN,$file);
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
373 while(<IN>)
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
374 {
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
375 ($gene,$pos,$res,$EPIseq,$num,$HLA)=(split);
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
376 push(@{$data{"$gene$pos$res"}},"$EPIseq,$num,$HLA");
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
377 }
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
378 return \%data;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
379 }
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
380
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
381 sub read_hyphy
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
382 {
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
383 $file=$_[0];
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
384 die("$file does not exist") unless -e $file;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
385 my %data=();
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
386 open(IN,$file);
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
387 %kw=("fel"=>1,"meme"=>1,"kind"=>1);#"betas"=>1);
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
388 while(<IN>)
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
389 {
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
390 chomp();
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
391 ($gene,$pos,@annot)=(split(/\,/));
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
392 foreach $a (@annot)
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
393 {
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
394 ($ref,$key)=(split(/\:/,$a));
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
395 $ref=~s/{//g;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
396 $ref=~s/}//g;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
397 $ref=~s/"//g;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
398 $key=~s/"//g;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
399 #print "$gene $pos $ref $key\n";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
400 if ($kw{$ref})
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
401 {
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
402 $data{"$gene$pos"}.="$ref:$key;";
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
403 }
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
404 }
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
405
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
406 }
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
407 return \%data;
1886293c2243 Uploaded
elixir-it
parents:
diff changeset
408 }