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