forked from rorycraig337/mouse_mm10_conserved_elements
-
Notifications
You must be signed in to change notification settings - Fork 0
/
neutralise_species.pl
55 lines (45 loc) · 1.16 KB
/
neutralise_species.pl
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
#!/usr/bin/perl
use warnings;
use strict;
use Getopt::Long qw(GetOptions);
#script replaces all bases for given species with Ns, useful for neutralising sequences prior to phastCons conserved element identification
#also cleans up headers to contain only species names (required for phastCons)
#usage: perl neutralise_species.pl --mfa in.mfa --species species1,species2 --out out.mfa
my $mfa;
my @species;
my $out;
GetOptions(
'mfa=s' => \$mfa,
'species=s' => \@species,
'out=s' => \$out,
) or die "missing input\n";
@species = split(/,/,join(',',@species));
open (IN, "$mfa") or die;
open (OUT, ">$out") or die;
my $mask_flag = 0;
while (my $line = <IN>) {
chomp $line;
if ($line =~ /^>/) {
$mask_flag = 0;
my @header = split /\//, $line;
print OUT "$header[0]\n"; #print only header before "/"
my $match = substr $header[0], 1; #remove ">"
foreach my $target (@species) {
if ($match eq "$target") {
$mask_flag = 1; #following sequence needs to be masked
}
}
}
else {
if ($mask_flag == 1) {
$line =~ tr/A-Za-z/N/; #change all text characters to Ns
print OUT "$line\n";
}
else {
print OUT "$line\n";
}
}
}
close IN;
close OUT;
exit;