-
Notifications
You must be signed in to change notification settings - Fork 16
/
concatenate_alignments.pl
151 lines (112 loc) · 5.01 KB
/
concatenate_alignments.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
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
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
#!/usr/bin/env perl
#####################################
## Robin van der Lee ##
## robinvanderlee AT gmail DOT com ##
############################################################################################################
## Genome-scale detection of positive selection in 9 primates predicts human-virus evolutionary conflicts ##
## Robin van der Lee, Laurens Wiel, Teunis J.P. van Dam, Martijn A. Huynen ##
############################################################################################################
use warnings;
use strict;
use Bio::Seq;
use Bio::SeqIO;
use Bio::AlignIO;
use Bio::Align::Utilities;
require('functions.pl');
### output directory for sequence files
my $OUTPUT_ROOT = "sequences";
my $OUTPUT_CDS = "cds";
my $OUTPUT_MASKED_ALIGNMENTS = "$OUTPUT_ROOT/prank-codon-masked";
###
my $BASE_DIRECTORY = "$OUTPUT_MASKED_ALIGNMENTS";
opendir(BASEDIR, $BASE_DIRECTORY) or die $!;
###########################
####### MAIN SCRIPT #######
###########################
print STDERR "Concatenating the alignments in $BASE_DIRECTORY into one alignment...\n\n";
my $count = 0;
my @array_of_alignments = ();
### loop over all sorted alignments results file in the base directory
while (my $basedir_file = readdir(BASEDIR)) {
next unless(-f "$BASE_DIRECTORY/$basedir_file"); # only loop over files
# only loop over relevant files, i.e. the # GUIDANCE & TCS MASKED ALIGNMENT, SORTED BY SPECIES
if($basedir_file =~ /^(ENSG\d+)__cds\.prank\-codon\-guidance\-tcs\-masked\-species\-sorted\.aln\.fa$/){
### PRINT COUNTER
$count++;
print STDERR $count . " MSAs read\r" if (($count % 100) == 0);
my $current_ensembl_id = $1;
### GET ALIGNMENT
my $current_alignment = undef;
eval { $current_alignment = readFastaFile("$OUTPUT_MASKED_ALIGNMENTS" . "/" . $basedir_file) }; # GUIDANCE & TCS MASKED ALIGNMENT, SORTED BY SPECIES
die if ($@);
# print $basedir_file . "\t" . $current_ensembl_id . "\t" . $current_alignment->num_sequences . "\n";
### push to alignment array
push(@array_of_alignments, $current_alignment);
}
}
print STDERR $count . " MSAs read\n";
print scalar(@array_of_alignments) . " alignments in array\n";
closedir BASEDIR;
### concatenate alignments into one alignment
# Title : cat
# Usage : $aln123 = cat($aln1, $aln2, $aln3)
# Function : Concatenates alignment objects. Sequences are identified by id.
# An error will be thrown if the sequence ids are not unique in the
# first alignment. If any ids are not present or not unique in any
# of the additional alignments then those sequences are omitted from
# the concatenated alignment, and a warning is issued. An error will
# be thrown if any of the alignments are not flush, since
# concatenating such alignments is unlikely to make biological
# sense.
# Returns : A new Bio::SimpleAlign object
# Args : A list of Bio::SimpleAlign objects
#my $concatenated_alignment = Bio::Align::Utilities->cat($array_of_alignments[1], $array_of_alignments[10]); ### ERROR Can't locate object method "id" via package "Bio::Align::Utilities" at /mnt/datas/rvdlee/lib/lib/perl5/Bio/Align/Utilities.pm line 372.
my %sequences = ();
foreach my $alignment ( @array_of_alignments ) {
foreach my $seq ($alignment->each_seq) {
my $identifier = $seq->display_id;
my $species = "";
$species = $1 if($identifier =~ /^(\d+__\S+?)__/);
$sequences{$species} .= $seq->seq; # append the sequence
}
}
### write sequences
my $outfile = "$OUTPUT_ROOT/concatenated_alignment__9primates__cds.prank-codon-guidance-tcs-masked.aln.fa";
my $seqout = Bio::SeqIO->new
(-file => ">$outfile",
-format => 'fasta');
while( my ($seqname,$seqstr) = each %sequences ) {
my $seq = Bio::Seq->new(-id => $seqname, -seq => $seqstr);
$seqout->write_seq($seq);
}
print "File $outfile has been generated\n";
### FINISH
print STDERR "\nFINISHED\n";
#########################
####### FUNCTIONS #######
#########################
###########################################################################
sub readFastaFile{
my $file = shift @_;
#stop function if file doesn't exists
eval{ checkIfExists($file) };
die if($@);
#read alignment file
my $AlignIO_in = Bio::AlignIO->new(-file => $file ,
-format => 'fasta');
my $alignment = $AlignIO_in->next_aln(); # Returns : a Bio::Align::AlignI compliant object
$alignment->set_displayname_flat(); # Function : Makes all the sequences be displayed as just their name, not name/start-end
#throw error if file contains another alignment
my $alignment2 = $AlignIO_in->next_aln();
if(defined $alignment2){
die "$file contains multiple alignments";
}
return $alignment;
}
###########################################################################
sub checkIfExists{
my $file = shift @_;
if(! -f $file){
die "File $file does not exist";
}
}