How to match a list of strings to another file and print the preceding line for multiple matches?

user_10329477 0
Viewed 83

I have a file of peptide sequences in peptides.txt that I would like to match to my protein database human_proteins.fasta. I would like to match the peptide list to the protein database and fetch the protein ID which is in the preceding line. Some of the peptides have multiple matches to the protein database.

Ultimately, I would like to produce a table/dataframe like this:

Peptide No. of matches Sequence ID Protein Sequence
AAAAA 2 ENST0001 AAAAABCFMED
AAAAA 2 ENST0002 AAAAAXXX

The first few lines of my hypothetical protein database human_proteins.fasta look like this:

>ENST0001
AAAAABCFMED
>ENST0002
AAAAAXXX
>ENST0003
MGRVSGLVPSR

peptides.txt looks like this:

AAAAA
LSSPATLNSR
HETLTSLNLEK
GGGGNFGPGPGSNFR
VSEQGLIEILK
DFLAGGIAAAISK

I am using the following command in bash

while read line; do printf $line grep -B 1 $line ../databases/human_proteins.fasta < peptides.txt

and I am able to get output like this:

>ENST0001
AAAAABCFMED
>ENST0002
AAAAAXXX

However, I am having trouble processing the output into a table. Is there a nice solution in unix/bash that can solve this?

bash
1 Answers

How about doing it with python:

import pandas as pd
import re

seq = {}
peptide = {}
matches = {}
result = []

with open("human_proteins.fasta") as f:
    while True:
        id = f.readline().rstrip().lstrip(">")
        if not id: break
        protein = f.readline().rstrip()
        seq[id] = protein
        matches[id] = 0

with open("peptides.txt") as f:
    for i in f:
        i = i.rstrip()
        for id in seq:
            protein = seq[id]
            if re.match(i, protein):
                peptide[id] = i
                matches[id] += 1

for id in seq:
    if matches[id]:
        result.append([peptide[id], matches[id], id, seq[id]])

df = pd.DataFrame(result)
df.columns = ["Peptide", "No. of matches", "Sequence ID", "Protein Sequence"]

print(df)

Output for provided files:

  Peptide  No. of matches Sequence ID Protein Sequence
0   AAAAA               1    ENST0001      AAAAABCFMED
1   AAAAA               1    ENST0002         AAAAAXXX

Please note I have assumed the column of No. of matches means the count of the same Protein Sequence. If it should be the count of the same Peptide (then the values above will be 2), please let me know.

[EDIT]
If the column of No. of matches referres to the count of peptide found in the fasta file, here is an alternative:

import pandas as pd
import re

seq = {}
peptide = {}
matches = {}
result = []

with open("human_proteins.fasta") as f:
    while True:
        id = f.readline().rstrip().lstrip(">")
        if not id: break
        protein = f.readline().rstrip()
        seq[id] = protein

with open("peptides.txt") as f:
    for i in f:
        i = i.rstrip()
        matches[i] = 0
        for id in seq:
            protein = seq[id]
            if re.match(i, protein):
                peptide[id] = i
                matches[i] += 1

for id in peptide:
    result.append([peptide[id], matches[peptide[id]], id, seq[id]])

df = pd.DataFrame(result)
df.columns = ["Peptide", "No. of matches", "Sequence ID", "Protein Sequence"]

print(df)
user_85723800
Related