scripts/edit_all_alignments.pl at master · rhilt13/scripts · GitHub

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
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
#! /usr/bin/perl -w

use Data::Dumper;
use Bio::SeqIO;

open(IN,$ARGV[0]);		# GT_sub_fam.txt
while(<IN>){
	chomp;
	$hash{$_}=1;
}

$new=Bio::SeqIO->new(-file=>$ARGV[1], -format=>"fasta");	# gta_all.cons.mafft.edit.nogaps

while($seq=$new->next_seq){
	$i=$seq->id;
	$s=$seq->seq;
	# print "$i\n";
	@a=split(/\|/,$i);
# =ccc
	# print "$a[1]\n";
	if (scalar(@a)==2){
		next;
	}elsif (scalar(@a)==3){
# =ccc
# 		if (defined $hash{$a[1]}){
# next;
# # =ccc
# 			$filename="sub-sub/".$a[1].".cons";
# 			# print "$filename\n";
# 			$new1=Bio::SeqIO->new(-file=>"$filename", -format=>"fasta");
# 			open(OUT,">$filename.edit");
# # print "aa\n";
# 			while($seq1=$new1->next_seq){
# 				# print $seq1->id,"\n";
# 				$i1=$seq1->id;
# 				$s1=$seq1->seq;
# 				if ($i1 eq $i){
# 					($start, $len, $new_seq)=Parse_sub($s1,$s);
# 					print OUT ">$i1\n$new_seq\n";
# 				}else{
# 					$new_seq=substr $s1, $start, $len;
# 					print OUT ">$i1\n$new_seq\n";
# 					$start_hash{$i1}=$start;
# 					$len_hash{$i1}=$len;
# 					# print "hash\t$i1\t$start\t$len\n";
# 				}
# 			}
# 			close OUT;
# # =cut
# 		}else{
			# next;
			$filename="new_sub/".$a[1].".FASTA.press";
			# print "$filename\n";
			$new1=Bio::SeqIO->new(-file=>"$filename", -format=>"fasta");
# print "aa\n";
			open(OUT,">$filename.edit");
			while($seq1=$new1->next_seq){
				# print $seq1->id,"\n";
				$i1=$seq1->id;
				$s1=$seq1->seq;
				if ($i1 eq $i){
					($start, $len, $new_seq)=Parse_sub($s1,$s);
					print OUT ">$i1\n$new_seq\n";
				}else{
					$new_seq=substr $s1, $start, $len;
					print OUT ">$i1\n$new_seq\n";
				}
			}
			close OUT;
		# }
	}elsif (scalar(@a)==4){
next;
# =ccc
		$filename="sub-sub/".$a[1]."/".$a[2].".FASTA.press";
		$new1=Bio::SeqIO->new(-file=>"$filename", -format=>"fasta");
		open(OUT,">$filename.edit");
		while($seq1=$new1->next_seq){
			$i1=$seq1->id;
			$s1=$seq1->seq;
			if ($i1 eq $i){
				($start, $len, $new_seq)=Parse_sub($s1,$s);
				print OUT ">$i1\n$new_seq\n";
			}else{
				# print "$i1\t$start\t$len\n";
				$new_seq=substr $s1, $start, $len;
				print OUT ">$i1\n$new_seq\n";
			}
		}
		close OUT;
	}
# =cut
}

sub Parse_sub {
	($s1,$s)=@_;
	@gaps=();
	$s1_nogap=$s1;
	$s1_nogap=~s/-//g;
	# print "NOGAP: $s1_nogap\nGAP  : $s1\n";
	$g="-";
	@gap_pos= all_match_positions($g,$s1);
	# print "NEED : $s\n";
	@loc= all_match_positions($s,$s1_nogap);
	foreach $m(@loc){
		$start=@$m[0];
		$end=@$m[1];
		# print "@$m[0]---@$m[1]\n";
	}
	# print "$start\t---$end\n";
	foreach $p(@gap_pos){
		push @gaps,@$p[0];
		# print "@$p[0]---@$p[1]\n";
	}
	foreach $g(@gaps){
		if ($start >=$g){
			$start++;
		}
		if ($end >=$g){
			$end++;
		}
	}
	# print "$start\t$end\n";
	# print "\n";
	$len=$end-$start+1;
	$new_seq=substr $s1, $start, $len;
	return ($start, $len, $new_seq);
}

sub all_match_positions {
    my ($regex, $string) = @_;
    my @ret;
    while ($string =~ /($regex)/g) {
        push @ret,[(pos($string)-length $1),pos($string)-1];
    }
    return @ret
}