Subject: Program conseqal.f - conversion of sequence alignments

From: [email-masked] (Clifford Felder)
Date: Wed, 21 Oct 98 12:02:58 EET

Dear Andrej Sali and Marilyn Safran,
  Enclosed is a copy of my program conseqal.f, a FORTRAN program to inter-convert
between some common sequence alignment formats. Please post a copy on your site
or ftp server, and feel to redistribute in my name. Use at your own risk, no
guarantees.

C conseqal.f - program to convert formats of sequence alignment files between
C SwissModel, pir (Sali-Modeller), Biosym; clustal and pileup (the latter only
C to read in). Limits: 25 sequences, 1500 residues per sequence. Runs 
C interactively. WARNING: only writes the sequences, user must edit the header
C information manually afterwards.
C Note: when reading  PIR format you should remove the final * of each sequence
C and when writing it, you should add it back.
C Clifford Felder [email-masked] August 1998.
C Dept. of Structural Biology, Weizmann Institute, Rehovot Israel
C This program removed blanks common to all sequences, great for consolidating
C sequence aignment subsets.
C Program now also calculates sequence homology. Does NOT convert case.
      CHARACTER*1 SEQ(1500,25)
      CHARACTER*12 LAB(25)
      CHARACTER*75 INFIL,OUFIL
      LOGICAL*4 GAP
      WRITE(*,*) 'Enter name of input alignment file:'
      READ(5,100) INFIL
      WRITE(*,*) 'Enter name of output alignment file:'
      READ(5,100) OUFIL
100   FORMAT(A75)
      OPEN(UNIT=2,FILE=INFIL,FORM='FORMATTED',STATUS='OLD')
      WRITE(*,101)
101   FORMAT('Enter format code of input sequence alignment:'/
     1'   1 - Old Swiss-Model'/'   2 - pir (Sali-Modeller)'/
     2'   3 - Biosym (Homology)'/'   4 - clustal'/
     3'   5 - Pileup (GCG)'/'   6 - Pileup (GCG) (wide lines)'/
     4'   7 - Single long editable lines'/'   8 - New Swiss-Model'/
     5'   0 - exit')
      READ(5,103) INFOR
103   FORMAT(I1)
      IF(INFOR.LE.0.OR.INFOR.GT.8) GOTO 99
      WRITE(*,102)
102   FORMAT('Enter format code of output sequence alignment:'/
     1'   1 - Swiss-Model'/'   2 - pir (Sali-Modeller)'/
     2'   3 - Biosym (Homology)'/'   4 - calculate homology'/
     3'   5 - Single long editable lines'/'   6 - New Swiss-Model'/
     4'   0 - Exit')
      READ(5,103) IOFOR
      IF(IOFOR.LE.0.OR.IOFOR.GT.6) GOTO 99
      WRITE(*,*) 'How many sequences are there in the alignment?'
      READ(5,104) NSEQ
104   FORMAT(I2)
      IF(NSEQ.LT.2) NSEQ=2
      IF(NSEQ.GT.25) NSEQ=25
      WRITE(*,*) 'What is the maximum sequence length as read in?'
      READ(5,105) LSEQ
105   FORMAT(I4)
      IF(LSEQ.LT.1) LSEQ=1
      IF(LSEQ.GT.1500) LSEQ=1500
C
C read in the sequence alignment
      IF(INFOR.EQ.1) CALL INSWS(SEQ,LAB,NSEQ,LSEQ)
      IF(INFOR.EQ.2) CALL INPIR(SEQ,LAB,NSEQ,LSEQ)
      IF(INFOR.EQ.3) CALL IBIOS(SEQ,LAB,NSEQ,LSEQ)
      IF(INFOR.EQ.4) CALL INCLU(SEQ,LAB,NSEQ,LSEQ)
      IF(INFOR.EQ.5) CALL INPIL(SEQ,LAB,NSEQ,LSEQ)
      IF(INFOR.EQ.6) CALL INPILW(SEQ,LAB,NSEQ,LSEQ)
      IF(INFOR.EQ.8) CALL INSWN(SEQ,LAB,NSEQ,LSEQ)
      IF(INFOR.EQ.7) THEN
         DO 400 I=1,NSEQ
400      READ(2,401,END=402) (SEQ(J,I),J=1,LSEQ)
401      FORMAT(1500(A1))
         ENDIF
402   CONTINUE
C
      CLOSE(2)
      WRITE(6,106) LSEQ
106   FORMAT(' Computed actual length of sequence=',I4)
C remove common blanks from sequences
      I = 0
20    I = I + 1
c         write(11,*)i,lseq,seq(i,1)
          IF(I.GT.LSEQ) GOTO 24
          GAP=.TRUE.
          DO 21 J=1,NSEQ
             IF(SEQ(I,J).NE.' '.AND.SEQ(I,J).NE.'-') GAP=.FALSE.
21           CONTINUE
          IF(.NOT.GAP) GOTO 20
C                  corection for deleting last residue if blank in all seqs
               IF(I.EQ.LSEQ) GOTO 23
C
               I1=I+1
               DO 22 J=1,NSEQ
               DO 22 K=I1,LSEQ
22                 SEQ(K-1,J)=SEQ(K,J)
C                     adjustment to sequence length
23                 LSEQ=LSEQ-1
                   I = I - 1  
      GOTO 20
24    CONTINUE
      WRITE(6,106) LSEQ
C
C write out the reformatted alignment
      OPEN(UNIT=3,FILE=OUFIL,FORM='FORMATTED')
C
      IF(IOFOR.EQ.1) CALL OUSWS(SEQ,LAB,NSEQ,LSEQ)
      IF(IOFOR.EQ.2) CALL OUPIR(SEQ,LAB,NSEQ,LSEQ)
      IF(IOFOR.EQ.3) CALL OBIOS(SEQ,LAB,NSEQ,LSEQ)      
      IF(IOFOR.EQ.4) CALL SEQHOM(SEQ,LAB,NSEQ,LSEQ)      
      IF(IOFOR.EQ.6) CALL OUSWN(SEQ,LAB,NSEQ,LSEQ)
      IF(IOFOR.EQ.5) THEN
         DO 403 I=1,NSEQ
403      WRITE(3,401) (SEQ(J,I),J=1,LSEQ)
         ENDIF
99    STOP
      END
      SUBROUTINE INSWS(SEQ,LAB,NSEQ,LSEQ)
C read in Old SwissModel sequence alignment
      CHARACTER*1 SEQ(1500,25),QQ
      CHARACTER*12 LAB(25)
C header info
      DO 1 I=1,NSEQ
1     READ(2,100) LAB(I)
100   FORMAT(23X,A12)
C read in each section of sequence, keeping *'s intact
      J1 = 0
C
4     I1 = J1 + 1
      READ(2,101,END=7) QQ
101   FORMAT(A1)
      J1 = I1 + 79
      IF(J1.GT.LSEQ) J1=LSEQ
         DO 5 K=1,NSEQ
5        READ(2,102,END=7) (SEQ(J,K),J=I1,J1)
102      FORMAT(80A1) 
C        write(*,*)i1,j1  
         IF(J1.LT.LSEQ) GOTO 4
C
7     LSEQ=J1
      RETURN
      END
      SUBROUTINE INSWN(SEQ,LAB,NSEQ,LSEQ)
C read in New SwissModel sequence alignment
      CHARACTER*1 SEQ(1500,25),QQ
      CHARACTER*12 LAB(25)
C this program handles only the actual sequence lines, you must remove the 
C homology-star and secondary structure lines before feeding into this program.
C
C read in each section of sequence, keeping *'s intact
      J1 = 0
4     I1 = J1 + 1
      J1 = I1 + 49
      IF(J1.GT.LSEQ) J1=LSEQ
         DO 5 K=1,NSEQ
5        READ(2,102,END=7) LAB(K),(SEQ(J,K),J=I1,J1)
102      FORMAT(10X,A12,3X,5(1X,10A1)) 
C        write(*,*)i1,j1  
      READ(2,101,END=7) QQ
101   FORMAT(A1)
         IF(J1.LT.LSEQ) GOTO 4
C
7     LSEQ=J1
      RETURN
      END
      SUBROUTINE INPIR(SEQ,LAB,NSEQ,LSEQ)
C read sequence in pir or Sali Modeller format
      CHARACTER*1 SEQ(1500,25),QQ
      CHARACTER*12 LAB(25)
C reading in first sequence to get maximum extent
      READ(2,100) LAB(1)
100   FORMAT(4X,A12/)
      J1=0
C
4     I1=J1 + 1
      J1=I1+74
      IF(J1.GT.LSEQ) J1=LSEQ
        READ(2,102,END=6) (SEQ(J,1),J=I1,J1)
102     FORMAT(75A1)
        IF(SEQ(I1,1).EQ.' '.OR.SEQ(I1,1).EQ.'') THEN
               J1=I1-1
               GOTO 5
               ENDIF
           IF(SEQ(I1,1).EQ.'>') THEN
               BACKSPACE 2
               J1=I1-1
               GOTO 5
               ENDIF
C        write(*,*)seq(i1,1)
C        write(*,*)i1,j1  
         IF(J1.LT.LSEQ) GOTO 4
C
6        READ(2,103)QQ
103      FORMAT(A1)
           write(*,*)qq
           IF(QQ.EQ.'>') THEN
               BACKSPACE 2
               GOTO 5
               ENDIF
           IF(QQ.NE.' '.AND.QQ.NE.'') GOTO 6
5        LSEQ=J1
C  reading in the other sequences, assumed to be the same length
      DO 10 I=2,NSEQ
         READ(2,100) LAB(I)
         READ(2,102,END=11) (SEQ(J,I),J=1,LSEQ)
10       READ(2,103,END=11) QQ   
11    RETURN
      END
      SUBROUTINE IBIOS(SEQ,LAB,NSEQ,LSEQ)
C read alignment in Biosym homology format
      CHARACTER*1 SEQ(1500,25),QQ
      CHARACTER*12 LAB(25)
      CHARACTER*5 LB1
      J1=0
C
      DO 1 I=1,NSEQ
1     LAB(I)='            '
4     I1 = J1 + 1
      J1 = I1 + 49
      IF(J1.GT.LSEQ) J1=LSEQ
         DO 5 K=1,NSEQ
            READ(2,102,END=7) LB1,(SEQ(J,K),J=I1,J1)
102         FORMAT(A5,1X,50A1)
            IF(I1.EQ.1) LAB(K)(1:5)=LB1
5           CONTINUE
C        write(*,*)i1,j1  
         READ(2,101,END=7) QQ
101      FORMAT(A1)
         IF(J1.LT.LSEQ) GOTO 4
7        LSEQ=J1
C  convert dots to dashes
      DO 50 I=1,NSEQ
      DO 50 J=1,LSEQ
         IF(SEQ(J,I).EQ.'.') SEQ(J,I)='-'
50    CONTINUE
      RETURN
      END
      SUBROUTINE INCLU(SEQ,LAB,NSEQ,LSEQ)
C read in alignment in clustal format
      CHARACTER*1 SEQ(1500,25)
      CHARACTER*12 LAB(25),LB1
      J1=0
         READ(2,101) QQ
         READ(2,101) QQ
C
4     I1 = J1 + 1
      J1 = I1 + 59
      IF(J1.GT.LSEQ) J1=LSEQ
         READ(2,101,END=7) QQ
101      FORMAT(A1)
         DO 5 K=1,NSEQ
            READ(2,102,END=7) LB1,(SEQ(J,K),J=I1,J1)
102         FORMAT(A12,4X,60A1)
            IF(I1.EQ.1) LAB(K)=LB1
5           CONTINUE
         READ(2,101,END=7) QQ
         IF(J1.LT.LSEQ) GOTO 4
7        LSEQ=J1
      RETURN
      END
      SUBROUTINE INPIL(SEQ,LAB,NSEQ,LSEQ)
C read in alignment in GCG pileup format
      CHARACTER*1 SEQ(1500,25),QQ
      CHARACTER*12 LAB(25),LB1
      CHARACTER*2 QX
C skip over header
1     CONTINUE
      READ(2,100) QX
100   FORMAT(A2)
      IF(QX.NE.'//') GOTO 1
C
      READ(2,102) QQ
      READ(2,102) QQ
      J1=0
C
4     I1 = J1 + 1
      J1 = I1 + 49
      IF(J1.GT.LSEQ) J1=LSEQ
         DO 5 K=1,NSEQ
            READ(2,102,END=7) LB1,(SEQ(J,K),J=I1,J1)
102         FORMAT(A12,5(1X,10A1))
            IF(I1.EQ.1) LAB(K)=LB1
5           CONTINUE
         READ(2,102,END=7) QQ
         READ(2,102,END=7) QQ
101      FORMAT(A1)
         IF(J1.LT.LSEQ) GOTO 4
7        LSEQ=J1
C  convert dots to dashes
      DO 50 I=1,NSEQ
      DO 50 J=1,LSEQ
         IF(SEQ(J,I).EQ.'.') SEQ(J,I)='-'
50    CONTINUE
      RETURN
      END
      SUBROUTINE INPILW(SEQ,LAB,NSEQ,LSEQ)
C read in alignment in Wide-line GCG pileup format
      CHARACTER*1 SEQ(1500,25),QQ
      CHARACTER*12 LAB(25),LB1
      CHARACTER*2 QX
      READ(2,102) QQ
      J1=0
C
4     I1 = J1 + 1
      J1 = I1 + 99
      IF(J1.GT.LSEQ) J1=LSEQ
         DO 5 K=1,NSEQ
            READ(2,102,END=7) LB1,(SEQ(J,K),J=I1,J1)
102         FORMAT(A12,1X,10(1X,10A1))
            IF(I1.EQ.1) LAB(K)=LB1
5           CONTINUE
         READ(2,102,END=7) QQ
         READ(2,102,END=7) QQ
101      FORMAT(A1)
         IF(J1.LT.LSEQ) GOTO 4
7        LSEQ=J1
C  convert dots to dashes
      DO 50 I=1,NSEQ
      DO 50 J=1,LSEQ
         IF(SEQ(J,I).EQ.'.') SEQ(J,I)='-'
50    CONTINUE
      RETURN
      END
      SUBROUTINE OUSWS(SEQ,LAB,NSEQ,LSEQ)
C  write out an Old SwissModel sequence alignment
      CHARACTER*1 SEQ(1500,25),QQ
      CHARACTER*12 LAB(25)
      QQ=''
C header info
      WRITE(3,100) LAB(1)
100   FORMAT('New  structure       : ',A12)
      IF(NSEQ.EQ.2) GOTO 2
      I1=NSEQ-1
      DO 1 I=2,I1
1     WRITE(3,101) LAB(I)
101   FORMAT('3D-Matched structure : ',A12)
2     CONTINUE
      WRITE(3,102) LAB(NSEQ),NSEQ
102   FORMAT('Reference structure  : ',A12/I5)
C
C write out each section of sequence
      J1 = 0
C
4     I1 = J1 + 1
      J1 = I1 + 79
      IF(J1.GT.LSEQ) J1=LSEQ
         DO 5 K=1,NSEQ
5        WRITE(3,105) (SEQ(J,K),J=I1,J1)
105      FORMAT(80A1)
      WRITE(3,106)
106   FORMAT(1X)
      IF(J1.LT.LSEQ) GOTO 4   
      RETURN
      END
      SUBROUTINE OUSWN(SEQ,LAB,NSEQ,LSEQ)
C  write out a New SwissModel sequence alignment (no homology or 2ndard struc)
      CHARACTER*1 SEQ(1500,25)
      CHARACTER*12 LAB(25)
C write out each section of sequence
      J1 = 0
C
4     I1 = J1 + 1
      J1 = I1 + 49
      IF(J1.GT.LSEQ) J1=LSEQ
         DO 5 K=1,NSEQ
5        WRITE(3,105) LAB(K)(1:6),I1,(SEQ(J,K),J=I1,J1)
105      FORMAT('SEQALI',4X,A6,1X,I6,2X,5(1X,10A1))
      WRITE(3,106)
106   FORMAT('SEQALI')
      IF(J1.LT.LSEQ) GOTO 4   
      RETURN
      END
      SUBROUTINE OUPIR(SEQ,LAB,NSEQ,LSEQ)
C
      CHARACTER*1 SEQ(1500,25),QQ
      CHARACTER*12 LAB(25)
      QQ=''
      DO 1 I=1,NSEQ
      WRITE(3,100) LAB(I),LAB(I),LSEQ
100   FORMAT('>P1;',A12/'structureX:',A12,':   4 : : ',I4,' : :XXXX:::')
      WRITE(3,101) (SEQ(J,I),J=1,LSEQ)
1     WRITE(3,102)
102   FORMAT(1X)
101   FORMAT(75A1)
      RETURN
      END
      SUBROUTINE OBIOS(SEQ,LAB,NSEQ,LSEQ)
C  output alignment in Biosym homology format
      CHARACTER*1 SEQ(1500,25),QQ
      CHARACTER*12 LAB(25)
      CHARACTER*5 LB1
      QQ=''
C  convert dots to dashes
      DO 50 I=1,NSEQ
      DO 50 J=1,LSEQ
         IF(SEQ(J,I).EQ.'-') SEQ(J,I)='.'
50    CONTINUE
      J1 = 0
4     I1 = J1 + 1
      J1 = I1 + 49
      IF(J1.GT.LSEQ) J1=LSEQ
         DO 5 K=1,NSEQ
         LB1=LAB(K)(1:5)
5        WRITE(3,105) LB1,(SEQ(J,K),J=I1,J1)
105      FORMAT(A5,1X,50A1)
      WRITE(3,106)
106   FORMAT(1X)
      IF(J1.LT.LSEQ) GOTO 4   
      RETURN
      END
      SUBROUTINE SEQHOM(SEQ,LAB,NSEQ,LSEQ)
C calculate pairwise sequence homology and identity
C the following sets of residues are considered homologous:
C ( A G ), ( T S ), ( D E ), ( K R ), ( F Y ), ( N Q ), ( I V L M )
C Of course, identical residues are also considered homologous.
      CHARACTER*1 SEQ(1500,25),C1,C2
      CHARACTER*12 LAB(25)
      LOGICAL HOM
      IF(NSEQ.LT.2) RETURN
      IIDG=0
      IHOMG=0
      ISEQG=0
      NSEQG=0
      I1 = NSEQ-1
      DO 1 I=1,I1
         J1=I+1
         DO 2 J=J1,NSEQ
            ISEQ=0
            IID=0
            IHOM=0
            DO 3 K=1,LSEQ
               C1=SEQ(K,I)
               C2=SEQ(K,J)
C eliminate from the count all places where both sequences are non-alphabetic
       IF((C1.LT.'A'.OR.C1.GT.'Z').AND.(C2.LT.'A'.OR.C2.GT.'Z')) GOTO 3
               ISEQ=ISEQ+1
               IF(C1.EQ.C2) IID=IID+1
               HOM=.FALSE.
               IF(C1.EQ.'A'.AND.C2.EQ.'G') HOM=.TRUE.
               IF(C1.EQ.'G'.AND.C2.EQ.'A') HOM=.TRUE.
               IF(C1.EQ.'D'.AND.C2.EQ.'E') HOM=.TRUE.
               IF(C1.EQ.'E'.AND.C2.EQ.'D') HOM=.TRUE.
               IF(C1.EQ.'T'.AND.C2.EQ.'S') HOM=.TRUE.
               IF(C1.EQ.'S'.AND.C2.EQ.'T') HOM=.TRUE.
               IF(C1.EQ.'K'.AND.C2.EQ.'R') HOM=.TRUE.
               IF(C1.EQ.'R'.AND.C2.EQ.'K') HOM=.TRUE.
               IF(C1.EQ.'F'.AND.C2.EQ.'Y') HOM=.TRUE.
               IF(C1.EQ.'Y'.AND.C2.EQ.'F') HOM=.TRUE.
               IF(C1.EQ.'N'.AND.C2.EQ.'Q') HOM=.TRUE.
               IF(C1.EQ.'Q'.AND.C2.EQ.'N') HOM=.TRUE.
               IF(C1.EQ.'I'.AND.C2.EQ.'V') HOM=.TRUE.
               IF(C1.EQ.'I'.AND.C2.EQ.'L') HOM=.TRUE.
               IF(C1.EQ.'I'.AND.C2.EQ.'M') HOM=.TRUE.
               IF(C1.EQ.'V'.AND.C2.EQ.'I') HOM=.TRUE.
               IF(C1.EQ.'L'.AND.C2.EQ.'I') HOM=.TRUE.
               IF(C1.EQ.'M'.AND.C2.EQ.'I') HOM=.TRUE.
               IF(C1.EQ.'V'.AND.C2.EQ.'L') HOM=.TRUE.
               IF(C1.EQ.'V'.AND.C2.EQ.'M') HOM=.TRUE.
               IF(C1.EQ.'L'.AND.C2.EQ.'M') HOM=.TRUE.
               IF(C1.EQ.'L'.AND.C2.EQ.'V') HOM=.TRUE.
               IF(C1.EQ.'M'.AND.C2.EQ.'V') HOM=.TRUE.
               IF(C1.EQ.'M'.AND.C2.EQ.'L') HOM=.TRUE.
               IF(HOM) IHOM=IHOM+1
3           CONTINUE
               IHOM=IHOM+IID
               PIID=IID
               PHOM=IHOM
               PSEQ=ISEQ
               PIID = 100.*PIID/PSEQ
               PHOM = 100.*PHOM/PSEQ
               WRITE(6,100) I,LAB(I),J,LAB(J)
               WRITE(6,101) IID,PIID
               WRITE(6,102) IHOM,PHOM,ISEQ
               WRITE(3,100) I,LAB(I),J,LAB(J)
100   FORMAT(' Sequences ',I2,1X,A12,'  and ',I2,1X,A12)
               WRITE(3,101) IID,PIID
101   FORMAT(' have ',I5,' residues or ',F6.1,
     1' percent identical, and')
               WRITE(3,102) IHOM,PHOM,ISEQ
102   FORMAT(6X,I5,' residues or ',F6.1,' percent homologous out of ',
     1I6,' residues total'/)
               IIDG=IIDG+IID
               IHOMG=IHOMG+IHOM
               ISEQG=ISEQG+ISEQ
               NSEQG=NSEQG+1
2        CONTINUE
1     CONTINUE
      JIID=IIDG/NSEQG
      JHOM=IHOMG/NSEQG
      JSEQ=ISEQG/NSEQG
      PIID=IIDG
      PHOM=IHOMG
      PSEQ=ISEQG
      PNSEQ=NSEQG
      PIID=100.*PIID/PSEQ
      PHOM=100.*PHOM/PSEQ
      WRITE(6,103)
      WRITE(6,104) JIID,PIID
      WRITE(6,105) JHOM,PHOM,JSEQ
      WRITE(3,103)
103   FORMAT(' Global sequence homology for all the sequences:')
      WRITE(3,104) JIID,PIID
104   FORMAT(I5,' residues/sequence or ',F6.1,
     1' percent are identical, and')
      WRITE(3,105) JHOM,PHOM,JSEQ
105   FORMAT(I5,' residues/sequence or ',F6.1,
     1' percent are homologous'/
     2 5X,' out of ',I5,' total residues/sequence') 
      RETURN
      END

Sincerely, Clifford Felder <[email-masked]>