1#!/usr/bin/perl -w
2
3#
4# Copyright 2011, Ben Langmead <langmea@cs.jhu.edu>
5#
6# This file is part of Bowtie 2.
7#
8# Bowtie 2 is free software: you can redistribute it and/or modify
9# it under the terms of the GNU General Public License as published by
10# the Free Software Foundation, either version 3 of the License, or
11# (at your option) any later version.
12#
13# Bowtie 2 is distributed in the hope that it will be useful,
14# but WITHOUT ANY WARRANTY; without even the implied warranty of
15# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16# GNU General Public License for more details.
17#
18# You should have received a copy of the GNU General Public License
19# along with Bowtie 2.  If not, see <http://www.gnu.org/licenses/>.
20#
21
22package DNA;
23use strict;
24use warnings;
25use Carp;
26
27my %revcompMap = (
28	"A" => "T", "a" => "t",
29	"T" => "A", "t" => "a",
30	"C" => "G", "c" => "g",
31	"G" => "C", "g" => "c",
32	"R" => "Y", "r" => "y",
33	"Y" => "R", "y" => "r",
34	"M" => "K", "m" => "k",
35	"K" => "M", "k" => "m",
36	"S" => "S", "s" => "s",
37	"W" => "W", "w" => "w",
38	"B" => "V", "b" => "v",
39	"V" => "B", "v" => "b",
40	"H" => "D", "h" => "d",
41	"D" => "H", "d" => "h",
42	"N" => "N", "n" => "n"
43);
44
45my %compat = (
46	"A" => "A",
47	"T" => "T",
48	"C" => "C",
49	"G" => "G",
50	"R" => "AG",
51	"Y" => "CT",
52	"M" => "AC",
53	"K" => "GT",
54	"S" => "CG",
55	"W" => "AT",
56	"B" => "CGT",
57	"V" => "ACG",
58	"H" => "ACT",
59	"D" => "AGT",
60	"N" => "N"
61);
62
63my %incompat = (
64	"A" => "CGT",
65	"T" => "ACG",
66	"C" => "AGT",
67	"G" => "ACT",
68	"R" => "CT",
69	"Y" => "AG",
70	"M" => "GT",
71	"K" => "AC",
72	"S" => "AT",
73	"W" => "CG",
74	"B" => "A",
75	"V" => "T",
76	"H" => "G",
77	"D" => "C",
78	"N" => "N"
79);
80
81my %unambigSet = (
82	"A" => 1, "a" => 1,
83	"C" => 1, "c" => 1,
84	"G" => 1, "g" => 1,
85	"T" => 1, "t" => 1
86);
87
88##
89# Return the complement, incl. if it's IUPAC.
90#
91sub comp($) {
92	my $ret = $revcompMap{$_[0]} || die "Can't reverse-complement '$_[0]'";
93	return $ret;
94}
95
96##
97# Return the complement, incl. if it's IUPAC.
98#
99sub revcomp {
100	my ($ret) = @_;
101	$ret = reverse $ret;
102	for(my $i = 0; $i < length($ret); $i++) {
103		substr($ret, $i, 1) = comp(substr($ret, $i, 1));
104	}
105	return $ret;
106}
107
108##
109# Return true iff it's unambiguous.
110#
111sub unambig($) {
112	return $unambigSet{$_[0]};
113}
114
115##
116# Manipulate DNA in an integer-indexed fashion.
117#
118sub plus($$) {
119	my ($c, $amt) = @_;
120	my %ctoi = ("A" => 0, "C" => 1, "G" => 2, "T" => 3);
121	my %itoc = (0 => "A", 1 => "C", 2 => "G", 3 => "T");
122	$c = uc $c;
123	defined($ctoi{$c}) || die;
124	return $itoc{($ctoi{$c}+$amt) % 4};
125}
126
1271;
128