Переименовать части заголовка fasta в соответствии с .csv

Question

Я хочу изменить части моего заголовка fastta со списком с частями .tsv.

Я не Биоинформатик, а просто Микробиолог с начинающими навыками в bash и python. Спасибо за помощь.

Пример:

Заголовок:

Prevalence_Sequence_ID: 1 | ARO: 3003072 | RES: mphL | Белковый гомолог Модель

с

.tsv

ARO: 3003072 mphL mphL - это кодируемые хромосомой макролидные фосфотрансферазы, которые инактивируют 14- и 15-членные макролиды, такие как эритромицин, кларитромицин, азитромицин.

до

Новый заголовок

Prevalence_Sequence_ID: 1 | mphL mphL - это кодируемые хромосомой макролидные фосфотрансферазы, которые инактивируют 14- и 15-членные макролиды, такие как эритромицин, кларитромицин, азитромицин. | RES: mphL Protein

Возможно, что ARO в заголовке fasta не указан в .tsv, тогда просто проигнорируйте его.

Примеры фаст

>Prevalence_Sequence_ID:1|ARO:3003072|RES:mphL|Protein Homolog Model
MTTLKVKQLANKKGLNILEDS
>gb|ARO:3004145|RES:AxyZ|Achromobacter_insuavis_AXX-A_
MARKTKEESQRTRDRILDAAEHVFLSKG
>Prevalence_Sequence_ID:31298|ARO:3000777|RES:adeF|Protein Homolog Model
MDFSRFFIDRPIFAAVLSILIFI

Пример .tsv

ARO:3003072 mphL    mphL is a chromosomally-encoded macrolide phosphotransferases that inactivate 14- and 15-membered macrolides such as erythromycin, clarithromycin, azithromycin.
ARO:3004145 AxyZ    AxyZ is a transcriptional regulator of the AxyXY-OprZ efflux pump system.
ARO:3000777 adeF    AdeF is the membrane fusion protein of the multidrug efflux complex AdeFGH.

Answer 1

Я сохранил данные вашего примера Fasta и TSV в example.fasta и example.tsv. Вот содержимое входного файла -

$ cat example.fasta
>Prevalence_Sequence_ID:1|ARO:3003072|RES:mphL|Protein Homolog Model
MTTLKVKQLANKKGLNILEDS
>gb|ARO:3004145|RES:AxyZ|Achromobacter_insuavis_AXX-A_
MARKTKEESQRTRDRILDAAEHVFLSKG
>Prevalence_Sequence_ID:31298|ARO:3000777|RES:adeF|Protein Homolog Model
MDFSRFFIDRPIFAAVLSILIFI

$ cat example.tsv
ARO:3003072 mphL    mphL is a chromosomally-encoded macrolide phosphotransferases that inactivate 14- and 15-membered macrolides such as erythromycin, clarithromycin, azithromycin.
ARO:3004145 AxyZ    AxyZ is a transcriptional regulator of the AxyXY-OprZ efflux pump system.
ARO:3000777 adeF    AdeF is the membrane fusion protein of the multidrug efflux complex AdeFGH.

# import biopython, bioython needs to be installed in your environment/machine
from Bio.SeqIO.FastaIO import SimpleFastaParser as sfp

# read in the tsv data into a dict
with open("example.tsv") as tsvdata:
    tsv_data = {line.strip().split("\t")[0]: " ".join(line.strip().split("\t")[1:])
                for line in tsvdata}

# read input fasta file contents and write to a separate file in real time
with open("example_out.fasta", "w") as outfasta:
    with open("example.fasta") as infasta:
        for header, seq in sfp(infasta):
            aro = header.strip().split("|")[1]  # get ARO for header
            header = header.replace(aro, tsv_data.get(aro, aro))  # lookup ARO in dict and replace if found, otherwise ignore it
            outfasta.write(">{0}\n{1}\n".format(header, seq))

Вот содержимое выходного файла -

$ cat example_out.fasta
>Prevalence_Sequence_ID:1|mphL mphL is a chromosomally-encoded macrolide phosphotransferases that inactivate 14- and 15-membered macrolides such as erythromycin, clarithromycin, azithromycin.|RES:mphL|Protein Homolog Model
MTTLKVKQLANKKGLNILEDS
>gb|AxyZ AxyZ is a transcriptional regulator of the AxyXY-OprZ efflux pump system.|RES:AxyZ|Achromobacter_insuavis_AXX-A_
MARKTKEESQRTRDRILDAAEHVFLSKG
>Prevalence_Sequence_ID:31298|adeF AdeF is the membrane fusion protein of the multidrug efflux complex AdeFGH.|RES:adeF|Protein Homolog Model
MDFSRFFIDRPIFAAVLSILIFI

Answer 2

Если вы используете awk , вы можете выполнить следующие шаги:

1. Считать полный файл tsv и сохранить все значения в массив, который индексируетсяпервый столбец.
2. Разбор файла fasta:
   • если вы встретите заголовок (начинается с >), то
     • извлеките ключ из заголовка (первая строка после|)
     • заменить ключ содержимым массива
   • напечатать текущую строку

Эти шаги также могут быть выполнены в python, но вы можете легко сделать это в awk со следующей строкой:

awk '# Read the TSV file only
     (NR==FNR) { key=$1; sub(key"[[:space:]]*",""); a[key]=$0; next }
     # From here process the fasta file
     # - if header, perform action:
     /^>/ { match($0,"[|][^|]+"); key=substr($0,RSTART+1,RLENGTH-1);sub(key,a[key]) }
     # print the current line
     1' index.tsv file.fasta > file.new.fasta

Answer 3

import pandas as pd
from Bio import SeqIO

tsvdata = pd.read_csv('example.tsv', sep='/t', header=None, names=['aro','_', 'description'])

for record in SeqIO.parse("example.fasta", "fasta"): 
    fasta_record = str(record).split('|')
    key = fasta_record[1]
    fasta_record[1]=tsvdata[tsvdata['aro']==key]['description'].values[0] 
    print('|'.join(fasta_record))

Answer 4

Если последовательность не нужно сортировать, мы можем отсортировать fasta со вторым полем, заменить первый пробел внутри .tsv на | в качестве разделителя и отсортировать его по первому полю, а затем объединить, используя правильный формат вывода:

cat <<EOF >fasta
>Prevalence_Sequence_ID:1|ARO:3003072|RES:mphL|Protein Homolog Model
MTTLKVKQLANKKGLNILEDS
>gb|ARO:3004145|RES:AxyZ|Achromobacter_insuavis_AXX-A_
MARKTKEESQRTRDRILDAAEHVFLSKG
>Prevalence_Sequence_ID:31298|ARO:3000777|RES:adeF|Protein Homolog Model
MDFSRFFIDRPIFAAVLSILIFI
EOF
cat <<EOF >tsv
ARO:3003072 mphL    mphL is a chromosomally-encoded macrolide phosphotransferases that inactivate 14- and 15-membered macrolides such as erythromycin, clarithromycin, azithromycin.
ARO:3004145 AxyZ    AxyZ is a transcriptional regulator of the AxyXY-OprZ efflux pump system.
ARO:3000777 adeF    AdeF is the membrane fusion protein of the multidrug efflux complex AdeFGH.
EOF

join -t'|' -12 -21 -o1.1,2.2,1.3 <(
    <fasta sort -t'|' -k2) <(
    <tsv sed 's/ /|/' | sort -t'|' -k1)

Если вам нужно отсортировать вывод в соответствии с fasta, мы можем пронумеровать строки с помощью nl -w1, затем объединить, а затем отсортировать выходные данные по номерам и удалить числа:

join -t'|' -12 -21 -o1.1,2.2,1.3 <(
    <fasta nl -w1 | sort -t'|' -k2) <(
    <tsv sed 's/ /|/' | sort -t'|' -k1) |
sort -t $'\t' -n -k2 | cut -f2-