API tools faq
paste
Advertisement
Guest User

Untitled

a guest
Jan 17th, 2017
121
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Java 8.45 KB | None | 0 0
raw download clone embed print report
  1. import java.util.*;
  2. import java.io.*;
  3. import java.lang.*;
  4.  
  5. public class GlobalAlign {
  6.  
  7.     // Gap penalties defined as negative values because addition is easier to work with in the algorithm
  8.     static int d = -10;
  9.     static int e = -5;
  10.     // BLOSUM50 matrix
  11.     static int[][] blosum50 = {
  12.         { 5, -2, -1, -2, -1, -1, -1,  0, -2, -1, -2, -1, -1, -3, -1,  1,  0, -3, -2,  0, -2, -2, -1, -1, -5},
  13.         {-2,  7, -1, -2, -4,  1,  0, -3,  0, -4, -3,  3, -2, -3, -3, -1, -1, -3, -1, -3, -1, -3,  0, -1, -5},
  14.         {-1, -1,  7,  2, -2,  0,  0,  0,  1, -3, -4,  0, -2, -4, -2,  1,  0, -4, -2, -3,  5, -4,  0, -1, -5},
  15.         {-2, -2,  2,  8, -4,  0,  2, -1, -1, -4, -4, -1, -4, -5, -1,  0, -1, -5, -3, -4,  6, -4,  1, -1, -5},
  16.         {-1, -4, -2, -4, 13, -3, -3, -3, -3, -2, -2, -3, -2, -2, -4, -1, -1, -5, -3, -1, -3, -2, -3, -1, -5},
  17.         {-1,  1,  0,  0, -3,  7,  2, -2,  1, -3, -2,  2,  0, -4, -1,  0, -1, -1, -1, -3,  0, -3,  4, -1, -5},
  18.         {-1,  0,  0,  2, -3,  2,  6, -3,  0, -4, -3,  1, -2, -3, -1, -1, -1, -3, -2, -3,  1, -3,  5, -1, -5},
  19.         { 0, -3,  0, -1, -3, -2, -3,  8, -2, -4, -4, -2, -3, -4, -2,  0, -2, -3, -3, -4, -1, -4, -2, -1, -5},
  20.         {-2,  0,  1, -1, -3,  1,  0, -2, 10, -4, -3,  0, -1, -1, -2, -1, -2, -3,  2, -4,  0, -3,  0, -1, -5},
  21.         {-1, -4, -3, -4, -2, -3, -4, -4, -4,  5,  2, -3,  2,  0, -3, -3, -1, -3, -1,  4, -4,  4, -3, -1, -5},
  22.         {-2, -3, -4, -4, -2, -2, -3, -4, -3,  2,  5, -3,  3,  1, -4, -3, -1, -2, -1,  1, -4,  4, -3, -1, -5},
  23.         {-1,  3,  0, -1, -3,  2,  1, -2,  0, -3, -3,  6, -2, -4, -1,  0, -1, -3, -2, -3,  0, -3,  1, -1, -5},
  24.         {-1, -2, -2, -4, -2,  0, -2, -3, -1,  2,  3, -2,  7,  0, -3, -2, -1, -1,  0,  1, -3,  2, -1, -1, -5},
  25.         {-3, -3, -4, -5, -2, -4, -3, -4, -1,  0,  1, -4,  0,  8, -4, -3, -2,  1,  4, -1, -4,  1, -4, -1, -5},
  26.         {-1, -3, -2, -1, -4, -1, -1, -2, -2, -3, -4, -1, -3, -4, 10, -1, -1, -4, -3, -3, -2, -3, -1, -1, -5},
  27.         { 1, -1,  1,  0, -1,  0, -1,  0, -1, -3, -3,  0, -2, -3, -1,  5,  2, -4, -2, -2,  0, -3,  0, -1, -5},
  28.         { 0, -1,  0, -1, -1, -1, -1, -2, -2, -1, -1, -1, -1, -2, -1,  2,  5, -3, -2,  0,  0, -1, -1, -1, -5},
  29.         {-3, -3, -4, -5, -5, -1, -3, -3, -3, -3, -2, -3, -1,  1, -4, -4, -3, 15,  2, -3, -5, -2, -2, -1, -5},
  30.         {-2, -1, -2, -3, -3, -1, -2, -3,  2, -1, -1, -2,  0,  4, -3, -2, -2,  2,  8, -1, -3, -1, -2, -1, -5},
  31.         { 0, -3, -3, -4, -1, -3, -3, -4, -4,  4,  1, -3,  1, -1, -3, -2,  0, -3, -1,  5, -3,  2, -3, -1, -5},
  32.         {-2, -1,  5,  6, -3,  0,  1, -1,  0, -4, -4,  0, -3, -4, -2,  0,  0, -5, -3, -3,  6, -4,  1, -1, -5},
  33.         {-2, -3, -4, -4, -2, -3, -3, -4, -3,  4,  4, -3,  2,  1, -3, -3, -1, -2, -1,  2, -4,  4, -3, -1, -5},
  34.         {-1,  0,  0,  1, -3,  4,  5, -2,  0, -3, -3,  1, -1, -4, -1,  0, -1, -2, -2, -3,  1, -3,  5, -1, -5},
  35.         {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -5},
  36.         {-5, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5,  1}
  37.     };
  38.     // Negative infinity
  39.     static int NEG_INF = -99999999;
  40.  
  41.     public static void main(String[] args) {
  42.         // System.out.println(blosum50.length);
  43.         // System.out.println(blosum50[0].length);
  44.         // System.out.println(getScore('F', 'M'));
  45.         if (args.length != 2) {
  46.             System.out.println("Please run the program with the reference sequence file as the first argument and the fasta file as the second argument.");
  47.             System.exit(0);
  48.         }
  49.  
  50.         // Initialize variables
  51.         String reference = "";
  52.         String temp = "";
  53.         ArrayList<FastaSeq> fastas = new ArrayList<FastaSeq>();
  54.  
  55.         // Load input files
  56.         File referenceFile = new File(args[0]);
  57.         File fastaFile = new File(args[1]);
  58.  
  59.         // Read the reference file
  60.         try {
  61.             Scanner sc = new Scanner(referenceFile);
  62.             reference = sc.nextLine();
  63.         } catch (Exception e) {
  64.             System.out.println("Reference file not found.");
  65.             System.exit(0);
  66.         }
  67.  
  68.         // Parse all of the fasta sequences into the Fasta array
  69.         try {
  70.             Scanner sc = new Scanner(fastaFile);
  71.             temp = sc.nextLine();
  72.             temp = temp.substring(1);
  73.             int counter = 0;
  74.             FastaSeq fasta = new FastaSeq(counter, temp);
  75.             while (sc.hasNextLine()) {
  76.                 temp = sc.nextLine();
  77.                 if (temp.startsWith(">")) {
  78.                     fastas.add(fasta);
  79.                     counter++;
  80.                     temp = temp.substring(1);
  81.                     fasta = fasta = new FastaSeq(counter, temp);
  82.                 } else {
  83.                     fasta.sequence += temp;
  84.                 }
  85.             }
  86.             counter++;
  87.             fastas.add(fasta);         
  88.         } catch (Exception e) {
  89.             System.out.println("Fasta file not found.");
  90.             System.exit(0);
  91.         }
  92.  
  93.         // Begin aligning each fasta sequence to the reference
  94.         for (int i = 0; i < fastas.size(); i++) {
  95.             FastaSeq f = fastas.get(i);
  96.  
  97.             // Declare DP matrices
  98.             int[][] m = new int[reference.length() + 1][f.sequence.length() + 1];
  99.             int[][] gx = new int[reference.length() + 1][f.sequence.length() + 1];
  100.             int[][] gy = new int[reference.length() + 1][f.sequence.length() + 1];
  101.  
  102.             // Initialize the matrices
  103.             for (int k = 0; k <= reference.length(); k++) {
  104.                 m[k][0] = NEG_INF;
  105.                 gy[k][0] = NEG_INF;
  106.             }
  107.             for (int k = 0; k <= f.sequence.length(); k++) {
  108.                 m[0][k] = NEG_INF;
  109.                 gx[0][k] = NEG_INF;
  110.             }
  111.             for (int k = 1; k <= reference.length(); k++) {
  112.                 gx[k][0] = d + ((k-1)*e);
  113.             }
  114.             for (int k = 1; k <= f.sequence.length(); k++) {
  115.                 gy[0][k] = d + ((k-1)*e);
  116.             }
  117.             m[0][0] = 0;
  118.  
  119.             // Begin filling the matrices using affine gap penalty
  120.             for (int k = 1; k <= reference.length(); k++) {
  121.                 for (int j = 1; j <= f.sequence.length(); j++) {
  122.                     m[k][j] = getScore(reference.charAt(k-1), f.sequence.charAt(j-1)) + Math.max(m[k - 1][j - 1], Math.max(gx[k - 1][j - 1], gy[k - 1][j - 1]));
  123.                     gx[k][j] = Math.max(m[k - 1][j] + d, Math.max(gx[k - 1][j] + e, gy[k - 1][j] + d));
  124.                     gy[k][j] = Math.max(m[k][j - 1] + d, Math.max(gx[k][j - 1] + d, gy[k][j - 1] + e));
  125.                 }
  126.             }
  127.  
  128.             // Set score to the maximum value from the 3 DP matrices
  129.             f.score = Math.max(m[reference.length()][f.sequence.length()], Math.max(gx[reference.length()][f.sequence.length()], gy[reference.length()][f.sequence.length()]));
  130.         }
  131.  
  132.         // Sort the fastas according to their scores
  133.         // FastaSeq implements Comparable, which enables sorting via Java Collections
  134.         Collections.sort(fastas);
  135.  
  136.         // Print the result
  137.         for (int i = 0; i < 3; i++) {
  138.             FastaSeq f = fastas.get(i);
  139.             System.out.print("Index=" + f.index);
  140.             System.out.print(" Name=" + f.name);
  141.             System.out.println(" Score=" + f.score);
  142.         }
  143.        
  144.     }
  145.  
  146.     // Given 2 amino acids, return the BLOSUM50 score of the amino acids
  147.     public static int getScore(char a, char b) {
  148.         int first = 0;
  149.         int second = 0;
  150.         switch(a) {
  151.             case 'A': first = 0; break;
  152.             case 'R': first = 1; break;
  153.             case 'N': first = 2; break;
  154.             case 'D': first = 3; break;
  155.             case 'C': first = 4; break;
  156.             case 'Q': first = 5; break;
  157.             case 'E': first = 6; break;
  158.             case 'G': first = 7; break;
  159.             case 'H': first = 8; break;
  160.             case 'I': first = 9; break;
  161.             case 'L': first = 10; break;
  162.             case 'K': first = 11; break;
  163.             case 'M': first = 12; break;
  164.             case 'F': first = 13; break;
  165.             case 'P': first = 14; break;
  166.             case 'S': first = 15; break;
  167.             case 'T': first = 16; break;
  168.             case 'W': first = 17; break;
  169.             case 'Y': first = 18; break;
  170.             case 'V': first = 19; break;
  171.             case 'B': first = 20; break;
  172.             case 'J': first = 21; break;
  173.             case 'Z': first = 22; break;
  174.             case 'X': first = 23; break;
  175.         }
  176.         switch(b) {
  177.             case 'A': second = 0; break;
  178.             case 'R': second = 1; break;
  179.             case 'N': second = 2; break;
  180.             case 'D': second = 3; break;
  181.             case 'C': second = 4; break;
  182.             case 'Q': second = 5; break;
  183.             case 'E': second = 6; break;
  184.             case 'G': second = 7; break;
  185.             case 'H': second = 8; break;
  186.             case 'I': second = 9; break;
  187.             case 'L': second = 10; break;
  188.             case 'K': second = 11; break;
  189.             case 'M': second = 12; break;
  190.             case 'F': second = 13; break;
  191.             case 'P': second = 14; break;
  192.             case 'S': second = 15; break;
  193.             case 'T': second = 16; break;
  194.             case 'W': second = 17; break;
  195.             case 'Y': second = 18; break;
  196.             case 'V': second = 19; break;
  197.             case 'B': second = 20; break;
  198.             case 'J': second = 21; break;
  199.             case 'Z': second = 22; break;
  200.             case 'X': second = 23; break;
  201.         }
  202.         return blosum50[first][second];
  203.     }
  204.  
  205.     // Object representation of a fasta sequence
  206.     public static class FastaSeq implements Comparable<FastaSeq> {
  207.         public int index;
  208.         public String name;
  209.         public String sequence;
  210.         public int score;
  211.  
  212.         public FastaSeq(int i, String n) {
  213.             index = i;
  214.             name = n;
  215.             sequence = "";
  216.         }
  217.  
  218.         public int compareTo(FastaSeq other) {
  219.             return (other.score - score);
  220.         }
  221.     }
  222. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!