search

veg / hyphy-analyses

HyPhy standalone analyses
MIT License
37 stars 17 forks source link

Question pre-msa.bf sequence order #36

Open kullrich opened 1 year ago

kullrich commented 1 year ago

Hi,

I have encountered that the sequence order for the processed sequences differ from the original input fasta file.

Is it possible to retain the order directly with pre-msa.bf?

Thank you in anticipation

hyphy --version
HYPHY 2.5.31(MP) for Linux on x86_64

grep ">" example.fas 
>P26
>P18_T-tropic
>P22
>P10_M-tropic
>P13_M-tropic
>C06_T-tropic
>P21
>P24
>C27
>C17
>C25
>C13_M-tropic
>C04
>C26_frameshift
>C10
>C23
>C11
>C24
>C20
>C08
>C07
>C03
>C09
>C18
>C01_M-tropic
>C16
>C21
>C21-copy

hyphy pre-msa.bf --input example.fas

grep ">" example.fas_nuc.fas 
>C26_frameshift
>C20
>C21
>C21_copy
>P26
>P18_T_tropic
>P22
>P10_M_tropic
>P13_M_tropic
>C06_T_tropic
>P21
>P24
>C27
>C17
>C25
>C13_M_tropic
>C04
>C10
>C23
>C11
>C24
>C08
>C07
>C03
>C09
>C18
>C01_M_tropic
>C16

grep ">" example.fas_protein.fas 
>C26_frameshift
>C20
>C21
>C21_copy
>P26
>P18_T_tropic
>P22
>P10_M_tropic
>P13_M_tropic
>C06_T_tropic
>P21
>P24
>C27
>C17
>C25
>C13_M_tropic
>C04
>C10
>C23
>C11
>C24
>C08
>C07
>C03
>C09
>C18
>C01_M_tropic
>C16

E.g. muscle per default re-orders the sequences and in the new version muscle5 the -stable option has been removed.

muscle -version
MUSCLE v3.8.31 by Robert C. Edgar

muscle -in example.fas_protein.fas -out example.fas_protein.msa

grep ">" example.fas_protein.msa
>P24
>C27
>C25
>C08
>C26_frameshift
>C17
>C04
>C13_M_tropic
>C10
>C23
>C11
>C24
>C03
>C20
>C07
>C01_M_tropic
>C16
>C21
>C21_copy
>C09
>C18
>C06_T_tropic
>P21
>P26
>P18_T_tropic
>P22
>P10_M_tropic
>P13_M_tropic

and

muscle5.1.linux_intel64 -version
muscle 5.1.linux64 [12f0e2]
Built Jan 13 2022 23:17:13

muscle5.1.linux_intel64 -align example.fas_protein.fas -output example.fas_protein.msa

grep ">" example.fas_protein.msa
>C21
>C04
>C17
>C21_copy
>P21
>C27
>P26
>C16
>C25
>C20
>C06_T_tropic
>P10_M_tropic
>P22
>C24
>C01_M_tropic
>C23
>P18_T_tropic
>C13_M_tropic
>C07
>P13_M_tropic
>C26_frameshift
>C11
>C08
>C03
>C09
>P24
>C10
>C18

mafft keeps the sequence order

MAFFT v7.505 (2022/Apr/10)

mafft example.fas_protein.fas > example.fas_protein.msa

>C26_frameshift
>C20
>C21
>C21_copy
>P26
>P18_T_tropic
>P22
>P10_M_tropic
>P13_M_tropic
>C06_T_tropic
>P21
>P24
>C27
>C17
>C25
>C13_M_tropic
>C04
>C10
>C23
>C11
>C24
>C08
>C07
>C03
>C09
>C18
>C01_M_tropic
>C16

Thank you in anticipation

Best regards

Kristian

spond commented 1 year ago

Dear @kullrich,

There is no option in pre-msa.bf to retain sequence ordering. For HyPhy, sequence ordering is immaterial. However, if you need to maintain such ordering for downstream analyses, perhaps the easiest thing to do is to re-order the output file using a post-processing script. I attach one for your reference that is written in python and assumes that BioPython is installed and the files are in FASTA format.

$python3.9 reorder.py -r original.fas -i scrambled.fas > unscrambled.fas

Best, Sergei

reorder.py.zip

kullrich commented 1 year ago

Thx for the quick reply and the downstream script.

I was just wondering, why e.g. pal2nal still fails using the constructed codon alignment and the protein alignment and complain about inconsistency between pep and nuc sequences:

hyphy pre-msa.bf --input example.fas
muscle -in example.fas_protein.fas -out example.fas_protein.unordered.msa
python reorder.py -r example.fas_protein.fas -i example.fas_protein.unordered.msa > example.fas_protein.msa
hyphy post-msa.bf --protein-msa example.fas_protein.msa --nucleotide-sequences example.fas_nuc.fas --output example-all.msa --compress No
# re-check with pal2nal
sed 's/?/-/g' example-all.msa > example-all.noquestion.msa
pal2nal example.fas_protein.msa example-all.noquestion.msa

Of course might be related how pal2nal handles gaps or codons that are only represented partially, like 'A--' or something alike.

Best regards

Kristian

#---  ERROR: inconsistency between the following pep and nuc seqs  ---#
>C26_frameshift
MRVKGIKKNCQGLWKWGMMLLGILMICSATEKLWVTVYYGVPVWKEATTTLFCASDAKAY
ETEVHNVWATHACVPTDPNPQEIVLENVTENFNMWKNNMVEQMQEDIISLWDQSLKPCVK
LTPLCVTLNCTKMMNVTNTNSSATTNTSSSENPMEEMKNCSFNITTHLRDQVXEXNMQLY
NLDLVPISDKNDSKYMLASCNTSVITQACPKVSFEPIPIHYCAPAGFAILKCNNKTFNGK
GPCTNVSTVQCTHGIKPVVSTQLLLNGSLAEKEIVIRSENLTNNAKTIIVQLNESVIINC
TRPNNNTRKSIHIQPGRAFYATGEIIGNIRQAYCILNETKWNNTLKQIVDKLREEFKNKT
IKFNPSSGGDPEIVMHTFNCGGEFFYCNTTKLFNSTWSVDGTWNGTKENSTIILPCKIKQ
IINMWQEVGKAMYAPPIKGQINCSSYITGLLLTRDGGYESRNGTEIFRPGGGDMRDNWRS
ELYKYKVVKIEPIGIAPTKAKRRVVQREKRAVGIGAVFLGFLGAAGSTMGAASITLTVQA
RQLLSGIVQQQNNLLRAIEAQQHLLQLTVWGIKQLQARLLAVERHLKDQQLLGIWGCSGK
LICTTAVPWNTSWSNKSLNQIWNNMTWMEWEREIDSYTGLIYSLIEESQNQQDKNEQELL
ALDHWASLWNWFSITNWLWYIKIFIMIVGGLIGLRIVFTVLSIVNRVRQGYSPLSFQTLL
PTQRGPDRPGGIEEEGGERDRDRSGQSVNGFLAIIWVDLRSLCLFLYRHLRDFLLIVTRT
VELLGLRGWEALKYLWNLLQYWIQELKNSAVSLLNAIAIAVAEGTDRVIEALQRACRAIL
HIPRRIRQGIERAVL
>C26_frameshift
ATGAGAGTGAAGGGGATCAAGAAGAATTGTCAGGGCTTGTGGAAATGGGGCATGATGCTC
CTTGGGATATTGATGATCTGTAGTGCTACAGAAAAATTGTGGGTCACAGTCTATTATGGG
GTACCTGTGTGGAAAGAAGCAACCACCACTCTATTCTGTGCATCAGATGCTAAAGCATAT
GAGACAGAGGTACATAATGTTTGGGCCACACATGCCTGTGTACCCACAGACCCCAACCCA
CAAGAAATAGTCTTGGAAAATGTGACAGAAAATTTTAATATGTGGAAAAATAACATGGTA
GAACAGATGCAAGAGGACATAATCAGTTTATGGGATCAAAGTCTAAAGCCCTGTGTAAAA
TTAACCCCACTCTGTGTCACTTTAAATTGCACTAAGATGATGAATGTTACTAACACCAAT
AGTAGCGCTACTACTAACACTAGTAGTAGCGAGAACCCGATGGAGGAAATGAAGAACTGC
TCTTTCAATATCACCACACACCTAAGAGATCAGGTGAAGAAAGAATATGCAACTTTATAA
CCTTGATTTAGTACCAATAAGTGATAAGAATGACAGCAAATATATGTTGGCAAGTTGTAA
CACCTCAGTTATTACACAAGCCTGTCCAAAGGTATCCTTTGAACCGATTCCCATACATTA
TTGTGCCCCGGCTGGTTTTGCGATTCTAAAGTGTAATAATAAAACGTTTAATGGAAAAGG
ACCATGTACAAATGTCAGCACAGTACAATGTACACATGGAATTAAGCCAGTAGTATCGAC
TCAACTGCTGTTAAATGGCAGTCTAGCAGAAAAGGAGATAGTAATTAGATCTGAAAATCT
CACAAACAATGCTAAAACTATAATAGTACAACTAAATGAATCTGTAATAATTAATTGTAC
AAGACCCAACAATAATACAAGAAAAAGTATACATATACAACCAGGGAGAGCATTTTATGC
AACAGGAGAAATAATAGGAAATATAAGACAAGCATATTGTATCCTCAATGAAACAAAATG
GAACAACACTTTAAAACAGATAGTTGATAAATTAAGAGAAGAGTTTAAGAATAAAACAAT
AAAATTTAATCCATCCTCAGGAGGAGACCCAGAAATTGTAATGCACACTTTTAATTGTGG
AGGGGAATTTTTCTACTGTAATACAACAAAACTGTTTAATAGTACTTGGAGTGTTGATGG
CACTTGGAATGGTACTAAAGAAAATAGCACCATCATACTCCCATGCAAAATAAAACAAAT
TATAAACATGTGGCAGGAAGTAGGAAAAGCAATGTATGCCCCTCCCATTAAAGGACAAAT
TAACTGTTCATCATATATTACGGGGCTGCTATTAACAAGAGATGGTGGTTATGAGAGTAG
GAACGGGACAGAGATTTTCAGACCTGGAGGAGGAGATATGAGGGACAATTGGAGAAGTGA
ATTATACAAATATAAAGTAGTAAAAATTGAACCAATAGGAATAGCACCCACCAAGGCAAA
GAGAAGAGTGGTGCAGAGAGAAAAAAGAGCAGTGGGAATAGGAGCTGTGTTCCTTGGGTT
CTTGGGAGCAGCAGGAAGCACTATGGGCGCAGCATCAATAACGCTGACGGTACAGGCCAG
ACAATTATTGTCTGGTATAGTGCAGCAGCAGAACAATCTGCTGAGGGCTATTGAGGCGCA
ACAGCATCTGTTGCAACTCACAGTCTGGGGCATCAAGCAGCTCCAGGCAAGACTCCTGGC
TGTGGAAAGACACCTAAAGGATCAACAGCTCCTGGGGATTTGGGGTTGCTCTGGAAAACT
CATTTGCACCACTGCTGTGCCTTGGAATACTAGTTGGAGTAATAAATCTCTGAATCAAAT
TTGGAATAACATGACCTGGATGGAGTGGGAAAGAGAGATTGACAGTTACACAGGCTTAAT
ATATTCCTTAATTGAAGAATCGCAGAACCAACAAGATAAAAATGAACAAGAATTATTAGC
ATTGGATCATTGGGCAAGTTTGTGGAATTGGTTTAGCATAACAAACTGGCTGTGGTATAT
AAAAATCTTTATAATGATAGTAGGAGGTTTAATAGGTTTAAGAATAGTTTTTACTGTGCT
TTCTATAGTGAATAGAGTTAGGCAGGGATACTCACCATTATCATTTCAGACCCTCCTCCC
AACCCAGAGGGGACCCGACAGGCCCGGAGGAATAGAAGAAGAAGGTGGAGAGAGAGACAG
AGACAGATCCGGACAATCAGTGAACGGATTCTTAGCAATTATCTGGGTGGACCTGCGGAG
CCTGTGCCTCTTCCTCTACCGCCACTTGAGAGACTTTCTCTTGATTGTAACGAGGACTGT
GGAACTCCTGGGACTCAGGGGGTGGGAAGCCCTCAAATACTTGTGGAACCTTCTACAGTA
TTGGATTCAGGAACTAAAGAATAGTGCTGTTAGCTTGCTCAATGCCATAGCCATAGCAGT
AGCTGAGGGAACAGATAGGGTCATAGAAGCATTACAAAGAGCTTGTAGAGCTATCCTCCA
CATACCTAGAAGAATCAGACAGGGCATCGAAAGGGCTGTGCTA

        ###-----   Check the following TBLASTN output:           -----###
        ###-----       your pep -> 'Query'                       -----###
        ###-----       your nuc -> 'Sbjct'                       -----###

Query= C26_frameshift
         (855 letters)

>C26_frameshift
          Length = 2563

 Score = 1385 bits (3584), Expect = 0.0
 Identities = 680/680 (100%), Positives = 680/680 (100%)
 Frame = +2

Query: 176  NMQLYNLDLVPISDKNDSKYMLASCNTSVITQACPKVSFEPIPIHYCAPAGFAILKCNNK 235
            NMQLYNLDLVPISDKNDSKYMLASCNTSVITQACPKVSFEPIPIHYCAPAGFAILKCNNK
Sbjct: 524  NMQLYNLDLVPISDKNDSKYMLASCNTSVITQACPKVSFEPIPIHYCAPAGFAILKCNNK 703

Query: 236  TFNGKGPCTNVSTVQCTHGIKPVVSTQLLLNGSLAEKEIVIRSENLTNNAKTIIVQLNES 295
            TFNGKGPCTNVSTVQCTHGIKPVVSTQLLLNGSLAEKEIVIRSENLTNNAKTIIVQLNES
Sbjct: 704  TFNGKGPCTNVSTVQCTHGIKPVVSTQLLLNGSLAEKEIVIRSENLTNNAKTIIVQLNES 883

Query: 296  VIINCTRPNNNTRKSIHIQPGRAFYATGEIIGNIRQAYCILNETKWNNTLKQIVDKLREE 355
            VIINCTRPNNNTRKSIHIQPGRAFYATGEIIGNIRQAYCILNETKWNNTLKQIVDKLREE
Sbjct: 884  VIINCTRPNNNTRKSIHIQPGRAFYATGEIIGNIRQAYCILNETKWNNTLKQIVDKLREE 1063

Query: 356  FKNKTIKFNPSSGGDPEIVMHTFNCGGEFFYCNTTKLFNSTWSVDGTWNGTKENSTIILP 415
            FKNKTIKFNPSSGGDPEIVMHTFNCGGEFFYCNTTKLFNSTWSVDGTWNGTKENSTIILP
Sbjct: 1064 FKNKTIKFNPSSGGDPEIVMHTFNCGGEFFYCNTTKLFNSTWSVDGTWNGTKENSTIILP 1243

Query: 416  CKIKQIINMWQEVGKAMYAPPIKGQINCSSYITGLLLTRDGGYESRNGTEIFRPGGGDMR 475
            CKIKQIINMWQEVGKAMYAPPIKGQINCSSYITGLLLTRDGGYESRNGTEIFRPGGGDMR
Sbjct: 1244 CKIKQIINMWQEVGKAMYAPPIKGQINCSSYITGLLLTRDGGYESRNGTEIFRPGGGDMR 1423

Query: 476  DNWRSELYKYKVVKIEPIGIAPTKAKRRVVQREKRAVGIGAVFLGFLGAAGSTMGAASIT 535
            DNWRSELYKYKVVKIEPIGIAPTKAKRRVVQREKRAVGIGAVFLGFLGAAGSTMGAASIT
Sbjct: 1424 DNWRSELYKYKVVKIEPIGIAPTKAKRRVVQREKRAVGIGAVFLGFLGAAGSTMGAASIT 1603

Query: 536  LTVQARQLLSGIVQQQNNLLRAIEAQQHLLQLTVWGIKQLQARLLAVERHLKDQQLLGIW 595
            LTVQARQLLSGIVQQQNNLLRAIEAQQHLLQLTVWGIKQLQARLLAVERHLKDQQLLGIW
Sbjct: 1604 LTVQARQLLSGIVQQQNNLLRAIEAQQHLLQLTVWGIKQLQARLLAVERHLKDQQLLGIW 1783

Query: 596  GCSGKLICTTAVPWNTSWSNKSLNQIWNNMTWMEWEREIDSYTGLIYSLIEESQNQQDKN 655
            GCSGKLICTTAVPWNTSWSNKSLNQIWNNMTWMEWEREIDSYTGLIYSLIEESQNQQDKN
Sbjct: 1784 GCSGKLICTTAVPWNTSWSNKSLNQIWNNMTWMEWEREIDSYTGLIYSLIEESQNQQDKN 1963

Query: 656  EQELLALDHWASLWNWFSITNWLWYIKIFIMIVGGLIGLRIVFTVLSIVNRVRQGYSPLS 715
            EQELLALDHWASLWNWFSITNWLWYIKIFIMIVGGLIGLRIVFTVLSIVNRVRQGYSPLS
Sbjct: 1964 EQELLALDHWASLWNWFSITNWLWYIKIFIMIVGGLIGLRIVFTVLSIVNRVRQGYSPLS 2143

Query: 716  FQTLLPTQRGPDRPGGIEEEGGERDRDRSGQSVNGFLAIIWVDLRSLCLFLYRHLRDFLL 775
            FQTLLPTQRGPDRPGGIEEEGGERDRDRSGQSVNGFLAIIWVDLRSLCLFLYRHLRDFLL
Sbjct: 2144 FQTLLPTQRGPDRPGGIEEEGGERDRDRSGQSVNGFLAIIWVDLRSLCLFLYRHLRDFLL 2323

Query: 776  IVTRTVELLGLRGWEALKYLWNLLQYWIQELKNSAVSLLNAIAIAVAEGTDRVIEALQRA 835
            IVTRTVELLGLRGWEALKYLWNLLQYWIQELKNSAVSLLNAIAIAVAEGTDRVIEALQRA
Sbjct: 2324 IVTRTVELLGLRGWEALKYLWNLLQYWIQELKNSAVSLLNAIAIAVAEGTDRVIEALQRA 2503

Query: 836  CRAILHIPRRIRQGIERAVL 855
            CRAILHIPRRIRQGIERAVL
Sbjct: 2504 CRAILHIPRRIRQGIERAVL 2563

 Score =  366 bits (939), Expect = e-105
 Identities = 172/172 (100%), Positives = 172/172 (100%)
 Frame = +1

Query: 1   MRVKGIKKNCQGLWKWGMMLLGILMICSATEKLWVTVYYGVPVWKEATTTLFCASDAKAY 60
           MRVKGIKKNCQGLWKWGMMLLGILMICSATEKLWVTVYYGVPVWKEATTTLFCASDAKAY
Sbjct: 1   MRVKGIKKNCQGLWKWGMMLLGILMICSATEKLWVTVYYGVPVWKEATTTLFCASDAKAY 180

Query: 61  ETEVHNVWATHACVPTDPNPQEIVLENVTENFNMWKNNMVEQMQEDIISLWDQSLKPCVK 120
           ETEVHNVWATHACVPTDPNPQEIVLENVTENFNMWKNNMVEQMQEDIISLWDQSLKPCVK
Sbjct: 181 ETEVHNVWATHACVPTDPNPQEIVLENVTENFNMWKNNMVEQMQEDIISLWDQSLKPCVK 360

Query: 121 LTPLCVTLNCTKMMNVTNTNSSATTNTSSSENPMEEMKNCSFNITTHLRDQV 172
           LTPLCVTLNCTKMMNVTNTNSSATTNTSSSENPMEEMKNCSFNITTHLRDQV
Sbjct: 361 LTPLCVTLNCTKMMNVTNTNSSATTNTSSSENPMEEMKNCSFNITTHLRDQV 516

 Score = 20.0 bits (40), Expect = 0.69
 Identities = 9/44 (20%), Positives = 23/44 (52%), Gaps = 2/44 (4%)
 Frame = -2

Query: 238  NGKGPCTNVSTVQCTHG--IKPVVSTQLLLNGSLAEKEIVIRSE 279
            NG+ PC  + T++ T    +KP+    +++   +   + V+ ++
Sbjct: 2139 NGEYPCLTLFTIESTVKTILKPIKPPTIIIKIFIYHSQFVMLNQ 2008

 Score = 17.7 bits (34), Expect = 3.4
 Identities = 9/21 (42%), Positives = 12/21 (57%)
 Frame = -2

Query: 835 ACRAILHIPRRIRQGIERAVL 855
           +C   LH+   +   IERAVL
Sbjct: 534 SCIFFLHLIS*VCGDIERAVL 472

 Score = 17.3 bits (33), Expect = 4.5
 Identities = 6/10 (60%), Positives = 7/10 (70%)
 Frame = -1

Query: 385  FYCNTTKLFN 394
            FYC   KLF+
Sbjct: 1084 FYCFILKLFS 1055

 Score = 17.3 bits (33), Expect = 4.5
 Identities = 7/20 (35%), Positives = 9/20 (45%)
 Frame = -3

Query: 243 CTNVSTVQCTHGIKPVVSTQ 262
           C N+ T  C H    +V  Q
Sbjct: 602 CYNLPTYICCHSYHLLVLNQ 543

 Score = 16.9 bits (32), Expect = 5.9
 Identities = 9/30 (30%), Positives = 15/30 (50%)
 Frame = -2

Query: 242 PCTNVSTVQCTHGIKPVVSTQLLLNGSLAE 271
           PC + + +   HG  P     LLL+  +A+
Sbjct: 759 PCVHCTVLTFVHG--PFPLNVLLLHFRIAK 676

 Score = 16.5 bits (31), Expect = 7.6
 Identities = 5/8 (62%), Positives = 7/8 (87%)
 Frame = -1

Query: 234  NKTFNGKG 241
            N +FNG+G
Sbjct: 1321 NLSFNGRG 1298

Lambda     K      H
   0.320    0.134    0.420 

Gapped
Lambda     K      H
   0.267   0.0410    0.140 

Matrix: BLOSUM62
Gap Penalties: Existence: 11, Extension: 1
Number of Sequences: 1
Number of Hits to DB: 4259
Number of extensions: 88
Number of successful extensions: 13
Number of sequences better than 10.0: 1
Number of HSP's gapped: 9
Number of HSP's successfully gapped: 9
Length of query: 855
Length of database: 854
Length adjustment: 42
Effective length of query: 813
Effective length of database: 812
Effective search space:   660156
Effective search space used:   660156
Neighboring words threshold: 13
Window for multiple hits: 40
X1: 16 ( 7.4 bits)
X2: 38 (14.6 bits)
X3: 64 (24.7 bits)
S1: 30 (16.7 bits)
S2: 30 (16.2 bits)