|
0
|
1 #!/usr/bin/env perl
|
|
|
2
|
|
|
3 use strict;
|
|
|
4 use warnings;
|
|
|
5
|
|
|
6 # reports the variant in transcripts separately only if the AA change is different, or distance from splicing site is different
|
|
|
7 @ARGV == 2 or @ARGV == 3 or die "Usage: $0 <hgvs input.txt> <output.txt> [ignore splice distance diff]\n";
|
|
|
8
|
|
|
9 my %ftype_rank = ( protein_coding => 100,
|
|
|
10 processed_transcript => 90,
|
|
|
11 antisense => 80,
|
|
|
12 retained_intron => 70,
|
|
|
13 lincRNA => 60,
|
|
|
14 nonsense_mediated_decay => 50,
|
|
|
15 misc_enrichment_kit_target => 0);
|
|
|
16
|
|
|
17 my %mtype_rank = ( nonsense => 100,
|
|
|
18 frameshift => 99,
|
|
|
19 nonstop => 90,
|
|
|
20 missense => 80,
|
|
|
21 silent => 50,
|
|
|
22 "non-coding" => 40);
|
|
|
23 open(IN, $ARGV[0])
|
|
|
24 or die "Cannot open $ARGV[0] for reading: $!\n";
|
|
|
25 open(OUT, ">$ARGV[1]")
|
|
|
26 or die "Cannot open $ARGV[1] for writing: $!\n";
|
|
|
27 my $succinct = (@ARGV == 3);
|
|
|
28 my @lines;
|
|
|
29 my $last_chr = "";
|
|
|
30 my $last_pos = "";
|
|
|
31 my $last_alt = "";
|
|
|
32 my %buffered_F;
|
|
|
33 my %buffered_id_rank;
|
|
|
34 my $header = <IN>;
|
|
|
35 print OUT $header;
|
|
|
36
|
|
|
37 my ($chr_column, $pos_column, $alt_column, $cdna_hgvs_column, $aa_hgvs_column, $phase_column, $transcript_column, $transcript_length_column, $exon_dist_column, $ftype_column, $mtype_column, $splicing_score_column, $splicing_effect_column, $sources_column);
|
|
|
38 chomp $header;
|
|
|
39 my @headers = split /\t/, $header;
|
|
|
40 for(my $i = 0; $i <= $#headers; $i++){
|
|
|
41 if($headers[$i] eq "Chr"){
|
|
|
42 $chr_column = $i;
|
|
|
43 }
|
|
|
44 elsif($headers[$i] eq "Feature type"){
|
|
|
45 $ftype_column = $i;
|
|
|
46 }
|
|
|
47 elsif($headers[$i] eq "Variant type"){
|
|
|
48 $mtype_column = $i;
|
|
|
49 }
|
|
|
50 elsif($headers[$i] eq "DNA From" or $headers[$i] eq "DNA Pos"){
|
|
|
51 $pos_column = $i;
|
|
|
52 }
|
|
|
53 elsif($headers[$i] eq "Obs base"){
|
|
|
54 $alt_column = $i;
|
|
|
55 }
|
|
|
56 elsif($headers[$i] eq "Transcript HGVS"){
|
|
|
57 $cdna_hgvs_column = $i;
|
|
|
58 }
|
|
|
59 elsif($headers[$i] eq "Protein HGVS"){
|
|
|
60 $aa_hgvs_column = $i;
|
|
|
61 }
|
|
|
62 elsif($headers[$i] eq "Phase"){
|
|
|
63 $phase_column = $i;
|
|
|
64 }
|
|
|
65 elsif($headers[$i] eq "Selected transcript"){
|
|
|
66 $transcript_column = $i;
|
|
|
67 }
|
|
|
68 elsif($headers[$i] eq "Transcript length"){
|
|
|
69 $transcript_length_column = $i;
|
|
|
70 }
|
|
|
71 elsif($headers[$i] eq "Closest exon junction (AA coding variants)"){
|
|
|
72 $exon_dist_column = $i;
|
|
|
73 }
|
|
|
74 elsif($headers[$i] eq "Splicing alteration potential"){
|
|
|
75 $splicing_score_column = $i;
|
|
|
76 }
|
|
|
77 elsif($headers[$i] eq "Splicing alteration details"){
|
|
|
78 $splicing_effect_column = $i;
|
|
|
79 }
|
|
|
80 elsif($headers[$i] eq "Sources"){
|
|
|
81 $sources_column = $i;
|
|
|
82 }
|
|
|
83 }
|
|
|
84 if(not defined $chr_column){
|
|
|
85 die "Could not find Chr header in $ARGV[0], aborting.\n";
|
|
|
86 }
|
|
|
87 if(not defined $pos_column){
|
|
|
88 die "Could not find 'DNA From' header in $ARGV[0], aborting.\n";
|
|
|
89 }
|
|
|
90 if(not defined $alt_column){
|
|
|
91 die "Could not find 'Obs base' header in $ARGV[0], aborting.\n";
|
|
|
92 }
|
|
|
93 if(not defined $cdna_hgvs_column){
|
|
|
94 die "Could not find 'Transcript HGVS' header in $ARGV[0], aborting.\n";
|
|
|
95 }
|
|
|
96 if(not defined $aa_hgvs_column){
|
|
|
97 die "Could not find 'Protein HGVS' header in $ARGV[0], aborting.\n";
|
|
|
98 }
|
|
|
99 if(not defined $phase_column){
|
|
|
100 die "Could not find Phase header in $ARGV[0], aborting.\n";
|
|
|
101 }
|
|
|
102 if(not defined $transcript_column){
|
|
|
103 die "Could not find 'Selected transcript' header in $ARGV[0], aborting.\n";
|
|
|
104 }
|
|
|
105 if(not defined $transcript_length_column){
|
|
|
106 die "Could not find 'Transcript length' header in $ARGV[0], aborting.\n";
|
|
|
107 }
|
|
|
108 #if(not defined $mtype_column){
|
|
|
109 # die "Could not find 'Variant type' header in $ARGV[0], aborting.\n";
|
|
|
110 #}
|
|
|
111 if(not defined $ftype_column){
|
|
|
112 die "Could not find 'Feature type' header in $ARGV[0], aborting.\n";
|
|
|
113 }
|
|
|
114 if(not defined $exon_dist_column){
|
|
|
115 die "Could not find 'Closest exon junction (AA coding variants)' header in $ARGV[0], aborting.\n";
|
|
|
116 }
|
|
|
117
|
|
|
118 while(<IN>){
|
|
|
119 chomp;
|
|
|
120 my @F = split /\t/, $_, -1;
|
|
|
121 my $chr = $F[$chr_column];
|
|
|
122 my $pos = $F[$pos_column];
|
|
|
123 my $alt = $F[$alt_column];
|
|
|
124 my $cdna_hgvs = $F[$cdna_hgvs_column];
|
|
|
125 my $hgvs = $F[$aa_hgvs_column];
|
|
|
126 my $ftype = $F[$ftype_column];
|
|
|
127 my $mtype;
|
|
|
128 $mtype = $F[$mtype_column] if defined $mtype_column;
|
|
|
129 $hgvs =~ s/\d+//g; # look only at the non-position parts of the AA syntax
|
|
|
130 my $exon_edge_distance = $F[$exon_dist_column];
|
|
|
131 my $phase = $F[$phase_column] || "";
|
|
|
132 # in the case of large indels (i.e. CNVs), their positions may not be the same, but effectively they should be reported as such
|
|
|
133 if($phase =~ /CNV-\S+?:(\S+)/){
|
|
|
134 $pos = $1;
|
|
|
135 }
|
|
|
136 my $preferred_id = $succinct ? ($F[$transcript_column] =~ /^$succinct/o) : 0;
|
|
|
137 my $id_rank = $F[$transcript_column];
|
|
|
138 $id_rank =~ s/^.*?0*(\d+)(\.\d+)?$/$1/; # look at only trailing non-padded number (and no .version suffix)
|
|
|
139
|
|
|
140 my $collapse_key = $succinct ? "$chr:$pos:$alt" : "$chr:$pos:$hgvs:$exon_edge_distance:$phase";
|
|
|
141 # Same variant
|
|
|
142 if($chr eq $last_chr and $pos eq $last_pos and $alt eq $last_alt){
|
|
|
143 # same AA effect
|
|
|
144 if(not exists $buffered_F{$collapse_key}){
|
|
|
145 $buffered_F{$collapse_key} = \@F;
|
|
|
146 $buffered_id_rank{$collapse_key} = $id_rank;
|
|
|
147 }
|
|
|
148 else{
|
|
|
149 $buffered_F{$collapse_key}->[$sources_column] .= "; $F[$sources_column]" if defined $sources_column;
|
|
|
150 $ftype_rank{$ftype} = 0 if not defined $ftype_rank{$ftype};
|
|
|
151 $mtype_rank{$mtype} = 0 if defined $mtype and not exists $mtype_rank{$mtype};
|
|
|
152 my $buf_ftype = $buffered_F{$collapse_key}->[$ftype_column];
|
|
|
153 $ftype_rank{$buf_ftype} = 0 if not defined $ftype_rank{$buf_ftype};
|
|
|
154 my $buf_mtype;
|
|
|
155 $buf_mtype = $buffered_F{$collapse_key}->[$mtype_column] if defined $mtype_column;
|
|
|
156 $mtype_rank{$buf_mtype} = 0 if defined $buf_mtype and not exists $mtype_rank{$buf_mtype};
|
|
|
157 # see if this transcript is "earlier" than the other, based on ID #
|
|
|
158 if($ftype_rank{$ftype} > $ftype_rank{$buf_ftype} or
|
|
|
159 (defined $mtype and $mtype_rank{$mtype} > $mtype_rank{$buf_mtype}) or
|
|
|
160 $preferred_id or
|
|
|
161 ($id_rank =~ /^\d+$/ and $buffered_id_rank{$collapse_key} =~ /^\d+$/ and $id_rank < $buffered_id_rank{$collapse_key})
|
|
|
162 or $F[$transcript_column] lt $buffered_F{$collapse_key}->[$transcript_column]){ # alphabetical as backup
|
|
|
163 # make this the first one reported (and use its HGVS syntax)
|
|
|
164 $buffered_F{$collapse_key}->[$transcript_length_column] = $F[$transcript_length_column]."; ".$buffered_F{$collapse_key}->[$transcript_length_column];
|
|
|
165 $buffered_F{$collapse_key}->[$transcript_column] = $F[$transcript_column]."; ".$buffered_F{$collapse_key}->[$transcript_column];
|
|
|
166 $buffered_F{$collapse_key}->[$cdna_hgvs_column] = $F[$cdna_hgvs_column];
|
|
|
167 $buffered_F{$collapse_key}->[$aa_hgvs_column] = $F[$aa_hgvs_column];
|
|
|
168 $buffered_id_rank{$collapse_key} = $id_rank;
|
|
|
169 }
|
|
|
170 else{
|
|
|
171 # just append it to the list of IDs
|
|
|
172 $buffered_F{$collapse_key}->[$transcript_length_column] .= "; $F[$transcript_length_column]";
|
|
|
173 $buffered_F{$collapse_key}->[$transcript_column] .= "; $F[$transcript_column]";
|
|
|
174 }
|
|
|
175 if($succinct and defined $splicing_score_column and $F[$splicing_score_column] ne "NA" and $F[$splicing_score_column] > $buffered_F{$collapse_key}->[$splicing_score_column]){
|
|
|
176 $buffered_F{$collapse_key}->[$splicing_score_column] = $F[$splicing_score_column];
|
|
|
177 $buffered_F{$collapse_key}->[$splicing_effect_column] = $F[$splicing_effect_column] ." (transcript model ".$F[$transcript_column].")";
|
|
|
178 }
|
|
|
179 }
|
|
|
180 }
|
|
|
181 # Different variant from the last line, dump any buffered data
|
|
|
182 else{
|
|
|
183 for my $collapse_key (keys %buffered_F){
|
|
|
184 if(defined $sources_column){
|
|
|
185 my %seen;
|
|
|
186 $buffered_F{$collapse_key}->[$sources_column] = join("; ", grep {not $seen{$_}++} split(/; /, $buffered_F{$collapse_key}->[$sources_column]));
|
|
|
187 }
|
|
|
188 print OUT join("\t", @{$buffered_F{$collapse_key}}), "\n";
|
|
|
189 }
|
|
|
190 undef %buffered_F;
|
|
|
191 $last_chr = $chr;
|
|
|
192 $last_pos = $pos;
|
|
|
193 $last_alt = $alt;
|
|
|
194 $buffered_F{$collapse_key} = \@F;
|
|
|
195 $buffered_id_rank{$collapse_key} = $id_rank;
|
|
|
196 }
|
|
|
197 }
|
|
|
198 for my $collapse_key (keys %buffered_F){
|
|
|
199 print OUT join("\t", @{$buffered_F{$collapse_key}}), "\n";
|
|
|
200 }
|
|
|
201 close(IN);
|
