Answers.Body: <p>Came across this looking for a solution myself last night but now I have created my own I will share it.</p>
  <p>Rather then using <code>Bio</code>python (which is a Python package I use a lot), it is more computationally efficient to parse the <code>.fasta</code> files from the protein data bank (<a href="https://www.rcsb.org/" rel="nofollow noreferrer">RCSB</a> in this case) which as far as I can see, contain this information. For example, this is <a href="https://www.rcsb.org/structure/2P16" rel="nofollow noreferrer">2P16</a>'s <code>.fasta</code> file (with 2 unique chains):</p>
  <pre><code>&gt;2P16_1|Chain A|Coagulation factor X (EC 3.4.21.6) (Stuart factor) (Stuart-Prower factor)|Homo sapiens (9606)
  IVGGQECKDGECPWQALLINEENEGFCGGTILSEFYILTAAHCLYQAKRFKVRVGDRNTEQEEGGEAVHEVEVVIKHNRFTKETYDFDIAVLRLKTPITFRMNVAPACLPERDWAESTLMTQKTGIVSGFGRTHEKGRQSTRLKMLEVPYVDRNSCKLSSSFIITQNMFCAGYDTKQEDACQGDSGGPHVTRFKDTYFVTGIVSWGEGCARKGKYGIYTKVTAFLKWIDRSMKT
  &gt;2P16_2|Chain B[auth L]|Coagulation factor X (EC 3.4.21.6) (Stuart factor) (Stuart-Prower factor)|Homo sapiens (9606)
  KLCSLDNGDCDQFCHEEQNSVVCSCARGYTLADNGKACIPTGPYPCGKQTLE
  </code></pre>
  <p>and <a href="https://www.rcsb.org/structure/6won" rel="nofollow noreferrer">6WON</a>'s (with 8 of the same chains):</p>
  <pre><code>&gt;6WON_1|Chains A, B, C, D, E, F, G, H|YohF|Salmonella typhimurium (strain LT2 / SGSC1412 / ATCC 700720) (99287)
  QGMTKVAIVTASDSGIGKACALLLAQNGFDIGITWHSDERGAQETAKKAAQFGVRAETIHLDLSQLPEGAQAIEHLIQRLGRVDVLVNNAGAMTKSAFIDMPFTQWRQIFTVDVDGAFLCAQIAARHMIKQGEGGRIINITSVHEHTPLPQASAYTAAKHALGGLTKSMALELIEHHILVNAVAPGAIATPMNDMDDSDIEPGSEPSIPIARPGSTHEIASLVAWLCSEGASYTTGQSLIVDGGFMLANPQFNAK
  </code></pre>
  <p>So far this always seems to be the case and parsing the <code>.fasta</code> is a lot easier then using BLAST or sequence similarity such as you are or were.</p>
  <p>For the implementation, whereas the <code>.pdb</code>, <code>.cif</code> files etc use the same download URL (<code>https://files.rcsb.org/download/</code>), the <code>.fasta</code> files use this link: <code>https://www.rcsb.org/fasta/entry/</code>. Parsing the downloaded <code>.fasta</code> files and only reading the lines that begin with <code>&gt;</code>, stripping out unnecessary chains (in my specific case I just want the first unique chains in the PDB so <code>A, B</code> being the same and <code>C, D</code> being the same would return me a list with [<code>A, C</code>], you can modify this script to your own ends) and running the Python code gives this output:</p>
  <pre><code>First unique chains in 1d2z (2 unique chains): ['B', 'A']
  First unique chains in 2p16 (2 unique chains): ['A', 'B']
  First unique chains in 6won (1 unique chain):  ['A']
  </code></pre>
  <p>There was an issue with some entries having <code>Chain</code> or <code>Chains</code> so both had to be considered. Additionally, there are commas that can be removed from all but the last chain ID but for my purposes I do not need to strip it, you can put something in <code>_fasta_helper()</code> to this end before <code>return</code>. This is the Python code and hopefully it helps someone out in the future:</p>
  <pre class="lang-py prettyprint-override"><code>import re
  from typing import Union
  import urllib.request
  from urllib.error import (HTTPError,
                            URLError,
                            ContentTooShortError as ContentError)
  
  def rcsb_downloader(pdb: str, ext: str = None, down_path: str = None) -&gt; str:
      &quot;&quot;&quot;Downloads the PDB from the protein databank RCSB.
  
      Parameters
      ----------
      pdb : str
          Protein databank accession 4 letter code i.e. '6won'
      ext : str, default = None
          Extension of the file to download i.e. '.fasta', '.cif' etc. Note: if
          None argument then ext = '.pdb'
  
      Returns
      -------
      str | bool
          Returns str path if download successful, False if failure
      &quot;&quot;&quot;
  
      if ext is None:
          ext = '.pdb'
      if down_path is None:
          download_path = ('/your/path/to/dir/' + pdb + ext)
      else:
          download_path = down_path + pdb + ext
  
      rcsb_path = 'https://files.rcsb.org/download/' + pdb + ext
  
      if ext == '.fasta': # fasta has a different web search then others
          rcsb_path = 'https://www.rcsb.org/fasta/entry/' + pdb
  
      try:
          urllib.request.urlretrieve(rcsb_path, download_path)
          return download_path
      # Try again in case there was a temporary issue with the first attempt
      except (ContentError, HTTPError, URLError):
          try:
              urllib.request.urlretrieve(rcsb_path, download_path)
              return download_path
          except (ContentError, HTTPError, URLError):
              print(f&quot;PDB failed to download {pdb}&quot;)
              return False
  
      def _fasta_helper(self, line: str, single_unique_chain: bool = True) -&gt; Union[str, list]:
          &quot;&quot;&quot;Separate the first chain ID from the .fasta comment line.&quot;&quot;&quot;
  
          # [1] some of the chains had [auth C] etc after so this is removed
          # [2] was originally a list comprehension that filtered based on the word
          # chain, however some of the fasta entries had &quot;Chain&quot; or &quot;Chains&quot; in their
          # description which broke this. The first arrow is removed, split at |, the
          # second entry (1st in Python) entry contains the chain information and the
          # first entry removed (which so far has always been chain, anticipate their
          # being edge cases with older PDB entries). This appears to hold up.
  
          clean_chain = [
              x.replace(&quot;,&quot;, &quot;&quot;)
              for x in list(
                  re.sub(r&quot;\[.*?\]&quot;, &quot;&quot;, line).lstrip(&quot;&gt;&quot;).split(&quot;|&quot;)[1].split(&quot; &quot;)[1:]  # [1]  # [2]
              )
          ]
  
          if not single_unique_chain:
              return clean_chain
          return clean_chain[0]
  
  def fasta_parse(path_: str) -&gt; list:
      &quot;&quot;&quot;Parse .fasta files and analyse lines beginning with '&gt;'.&quot;&quot;&quot;
  
      with open(path_, 'r', encoding='utf-8') as fasta:
          total_unique_chains = 0  # Use as a flag to test if more than one unique chain
          unique_chains       = []
          for line in fasta:
              if line[0] == '&gt;' and total_unique_chains == 0:
                  total_unique_chains += 1
                  unique_chains.append(_fasta_helper(line))
                  continue
              elif line[0] == '&gt;' and total_unique_chains != 0:
                  unique_chains.append(_fasta_helper(line))
      return unique_chains
  
  pdb_diff_chain  = '2p16' # Chain ...
  pdb_diff_chains = '1d2z' # Chains ...
  pdb_same_chains = '6won' # Chains ...
  
  download_p_diff  = rcsb_downloader(pdb_diff_chain, '.fasta')
  download_p_diffs = rcsb_downloader(pdb_diff_chains, '.fasta')
  download_p_same  = rcsb_downloader(pdb_same_chains, '.fasta')
  
  unique_s  = fasta_parse(download_p_same)
  unique_d  = fasta_parse(download_p_diffs)
  unique_ds = fasta_parse(download_p_diff)
  
  print(f'First unique chains in {pdb_diff_chains} (2 unique chains): {unique_ds}\n'
        f'First unique chains in {pdb_diff_chain} (2 unique chains): {unique_d}\n'
        f'First unique chains in {pdb_same_chains} (1 unique chain):  {unique_s}')
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-06-12T15:53:44.630Z
Answers.Text: This is a great way, but unfortunately I have the files as PDBs, so I'd have to re-download everything as FASTA files instead. However, I really like your solution, thank you for it! It's a great way to do it and definitely better than processing sequences from PDBs (due to chain breakings and its processing by BioPython.Polypeptide module), so I'll keep this in mind when I need to do it again.
Answers.UserDisplayName: kormi
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-06-07T13:34:33.653Z
Answers.LastActivityDate: 2023-07-15T15:47:23.660Z
Answers.LastEditDate: 2023-07-15T15:47:23.660Z
Answers.UserDisplayName: piJT
Body: <p>I'm trying to find unique chains or proteins within their PDB files. Many proteins have multiple chains, but very often they are identical (say, the PDB file consists of 4 homodimers, for example, so there is one unique chain represented 8 times total). I am processing certain structural features of the chains, so processing two identical chains is redundant for my pruposes.</p>
  <p>I am using Biopython and the way I'm currently identifying unique chains is looking through each PDB file, checking the sequences of the chains and comparing them to one another using alignments (specifically <code>pairwise2.align.globalxx(seq1, seq2)</code>. I consider a chain non-unique if it shares &gt;90% sequence similarity to another chain. (Most chains have 100+ residues.)</p>
  <p>My question is, is this threshold appropriate? Should I go lower, maybe 80%, or the opposite, higher? For now, it seems to work fine, but I'll be processing a lot of PDB files and want to be on the safer side. If I accidentally process identical chains, that's not much of an issue, but it adds unnecessary time to the processing. Is there maybe a <strong>source/article</strong> to help with what sequence similarity I should be using? I've found multiple articles, but they usually talk about homology (i.e., common ancestry of proteins), but I'm only doing this within the protein itself and its subunits (chains), so this is of no use for me.</p>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2022-06-09T14:18:10.350Z
Comments.Text: I commented yesterday w/ a link to SIFTS —sorry about the terseness: I meant to say why not use SIFTS table of PDB chain to Uniprot if you are just after finding which chains are the same. Alternatively, the PDB API data will tell you in a more convoluted way what is present. Neither will tell you about gaps in sequence and the latter will be weird for C-Alpha traces (as the precise AA is unknown)
Comments.UserDisplayName: Matteo Ferla
Tags: <phylogenetics><proteins><python><sequence>
Title: Identification of unique chains in proteins by sequence similarity
ThreadText: # Identification of unique chains in proteins by sequence similarity
  
  (Asked by: username_0 on 2022-06-08)
  
  I'm trying to find unique chains or proteins within their PDB files. Many proteins have multiple chains, but very often they are identical (say, the PDB file consists of 4 homodimers, for example, so there is one unique chain represented 8 times total). I am processing certain structural features of the chains, so processing two identical chains is redundant for my pruposes.
  
  I am using Biopython and the way I'm currently identifying unique chains is looking through each PDB file, checking the sequences of the chains and comparing them to one another using alignments (specifically `pairwise2.align.globalxx(seq1, seq2)`. I consider a chain non-unique if it shares >90% sequence similarity to another chain. (Most chains have 100+ residues.)
  
  My question is, is this threshold appropriate? Should I go lower, maybe 80%, or the opposite, higher? For now, it seems to work fine, but I'll be processing a lot of PDB files and want to be on the safer side. If I accidentally process identical chains, that's not much of an issue, but it adds unnecessary time to the processing. Is there maybe a **source/article** to help with what sequence similarity I should be using? I've found multiple articles, but they usually talk about homology (i.e., common ancestry of proteins), but I'm only doing this within the protein itself and its subunits (chains), so this is of no use for me.
  
  **Tags:** phylogenetics, proteins, python, sequence
  
  > **Comments:**
  > **username_1**: I commented yesterday w/ a link to SIFTS —sorry about the terseness: I meant to say why not use SIFTS table of PDB chain to Uniprot if you are just after finding which chains are the same. Alternatively, the PDB API data will tell you in a more convoluted way what is present. Neither will tell you about gaps in sequence and the latter will be weird for C-Alpha traces (as the precise AA is unknown)  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_2 on 2023-06-07. Score: 3 | Accepted answer
  
  Came across this looking for a solution myself last night but now I have created my own I will share it.
  
  Rather then using `Bio`python (which is a Python package I use a lot), it is more computationally efficient to parse the `.fasta` files from the protein data bank ([RCSB](https://www.rcsb.org/) in this case) which as far as I can see, contain this information. For example, this is [2P16](https://www.rcsb.org/structure/2P16)'s `.fasta` file (with 2 unique chains):
  
  ```
  >2P16_1|Chain A|Coagulation factor X (EC 3.4.21.6) (Stuart factor) (Stuart-Prower factor)|Homo sapiens (9606)
  IVGGQECKDGECPWQALLINEENEGFCGGTILSEFYILTAAHCLYQAKRFKVRVGDRNTEQEEGGEAVHEVEVVIKHNRFTKETYDFDIAVLRLKTPITFRMNVAPACLPERDWAESTLMTQKTGIVSGFGRTHEKGRQSTRLKMLEVPYVDRNSCKLSSSFIITQNMFCAGYDTKQEDACQGDSGGPHVTRFKDTYFVTGIVSWGEGCARKGKYGIYTKVTAFLKWIDRSMKT
  >2P16_2|Chain B[auth L]|Coagulation factor X (EC 3.4.21.6) (Stuart factor) (Stuart-Prower factor)|Homo sapiens (9606)
  KLCSLDNGDCDQFCHEEQNSVVCSCARGYTLADNGKACIPTGPYPCGKQTLE
  
  ```
  
  and [6WON](https://www.rcsb.org/structure/6won)'s (with 8 of the same chains):
  
  ```
  >6WON_1|Chains A, B, C, D, E, F, G, H|YohF|Salmonella typhimurium (strain LT2 / SGSC1412 / ATCC 700720) (99287)
  QGMTKVAIVTASDSGIGKACALLLAQNGFDIGITWHSDERGAQETAKKAAQFGVRAETIHLDLSQLPEGAQAIEHLIQRLGRVDVLVNNAGAMTKSAFIDMPFTQWRQIFTVDVDGAFLCAQIAARHMIKQGEGGRIINITSVHEHTPLPQASAYTAAKHALGGLTKSMALELIEHHILVNAVAPGAIATPMNDMDDSDIEPGSEPSIPIARPGSTHEIASLVAWLCSEGASYTTGQSLIVDGGFMLANPQFNAK
  
  ```
  
  So far this always seems to be the case and parsing the `.fasta` is a lot easier then using BLAST or sequence similarity such as you are or were.
  
  For the implementation, whereas the `.pdb`, `.cif` files etc use the same download URL (`https://files.rcsb.org/download/`), the `.fasta` files use this link: `https://www.rcsb.org/fasta/entry/`. Parsing the downloaded `.fasta` files and only reading the lines that begin with `>`, stripping out unnecessary chains (in my specific case I just want the first unique chains in the PDB so `A, B` being the same and `C, D` being the same would return me a list with [`A, C`], you can modify this script to your own ends) and running the Python code gives this output:
  
  ```
  First unique chains in 1d2z (2 unique chains): ['B', 'A']
  First unique chains in 2p16 (2 unique chains): ['A', 'B']
  First unique chains in 6won (1 unique chain):  ['A']
  
  ```
  
  There was an issue with some entries having `Chain` or `Chains` so both had to be considered. Additionally, there are commas that can be removed from all but the last chain ID but for my purposes I do not need to strip it, you can put something in `_fasta_helper()` to this end before `return`. This is the Python code and hopefully it helps someone out in the future:
  
  ```py
  import re
  from typing import Union
  import urllib.request
  from urllib.error import (HTTPError,
                            URLError,
                            ContentTooShortError as ContentError)
  
  def rcsb_downloader(pdb: str, ext: str = None, down_path: str = None) -> str:
      """Downloads the PDB from the protein databank RCSB.
  
      Parameters
      ----------
      pdb : str
          Protein databank accession 4 letter code i.e. '6won'
      ext : str, default = None
          Extension of the file to download i.e. '.fasta', '.cif' etc. Note: if
          None argument then ext = '.pdb'
  
      Returns
      -------
      str | bool
          Returns str path if download successful, False if failure
      """
  
      if ext is None:
          ext = '.pdb'
      if down_path is None:
          download_path = ('/your/path/to/dir/' + pdb + ext)
      else:
          download_path = down_path + pdb + ext
  
      rcsb_path = 'https://files.rcsb.org/download/' + pdb + ext
  
      if ext == '.fasta': # fasta has a different web search then others
          rcsb_path = 'https://www.rcsb.org/fasta/entry/' + pdb
  
      try:
          urllib.request.urlretrieve(rcsb_path, download_path)
          return download_path
      # Try again in case there was a temporary issue with the first attempt
      except (ContentError, HTTPError, URLError):
          try:
              urllib.request.urlretrieve(rcsb_path, download_path)
              return download_path
          except (ContentError, HTTPError, URLError):
              print(f"PDB failed to download {pdb}")
              return False
  
      def _fasta_helper(self, line: str, single_unique_chain: bool = True) -> Union[str, list]:
          """Separate the first chain ID from the .fasta comment line."""
  
          # [1] some of the chains had [auth C] etc after so this is removed
          # [2] was originally a list comprehension that filtered based on the word
          # chain, however some of the fasta entries had "Chain" or "Chains" in their
          # description which broke this. The first arrow is removed, split at |, the
          # second entry (1st in Python) entry contains the chain information and the
          # first entry removed (which so far has always been chain, anticipate their
          # being edge cases with older PDB entries). This appears to hold up.
  
          clean_chain = [
              x.replace(",", "")
              for x in list(
                  re.sub(r"\[.*?\]", "", line).lstrip(">").split("|")[1].split(" ")[1:]  # [1]  # [2]
              )
          ]
  
          if not single_unique_chain:
              return clean_chain
          return clean_chain[0]
  
  def fasta_parse(path_: str) -> list:
      """Parse .fasta files and analyse lines beginning with '>'."""
  
      with open(path_, 'r', encoding='utf-8') as fasta:
          total_unique_chains = 0  # Use as a flag to test if more than one unique chain
          unique_chains       = []
          for line in fasta:
              if line[0] == '>' and total_unique_chains == 0:
                  total_unique_chains += 1
                  unique_chains.append(_fasta_helper(line))
                  continue
              elif line[0] == '>' and total_unique_chains != 0:
                  unique_chains.append(_fasta_helper(line))
      return unique_chains
  
  pdb_diff_chain  = '2p16' # Chain ...
  pdb_diff_chains = '1d2z' # Chains ...
  pdb_same_chains = '6won' # Chains ...
  
  download_p_diff  = rcsb_downloader(pdb_diff_chain, '.fasta')
  download_p_diffs = rcsb_downloader(pdb_diff_chains, '.fasta')
  download_p_same  = rcsb_downloader(pdb_same_chains, '.fasta')
  
  unique_s  = fasta_parse(download_p_same)
  unique_d  = fasta_parse(download_p_diffs)
  unique_ds = fasta_parse(download_p_diff)
  
  print(f'First unique chains in {pdb_diff_chains} (2 unique chains): {unique_ds}\n'
        f'First unique chains in {pdb_diff_chain} (2 unique chains): {unique_d}\n'
        f'First unique chains in {pdb_same_chains} (1 unique chain):  {unique_s}')
  
  ```
  
  > **Comments:**
  > **username_0**: This is a great way, but unfortunately I have the files as PDBs, so I'd have to re-download everything as FASTA files instead. However, I really like your solution, thank you for it! It's a great way to do it and definitely better than processing sequences from PDBs (due to chain breakings and its processing by BioPython.Polypeptide module), so I'll keep this in mind when I need to do it again.  
  
  ---
  
  