Finden Sie Proteinhomologe mit BLASTp

Ich versuche mit BLASTp Homologe einer Reihe von Proteinen zu finden. Ich arbeite mit benutzerdefinierten Datenbanken.

Ich verwende einen Wert von 0,00001 als Schwellenwert.

Ich möchte Abfragen filtern, die Treffer mit >90 % Identitäten haben. Da die BLASTp-Ausgabe auf HSPs basiert, kann ich nicht nach %identities/query filtern, sondern nur nach HSP.

Ich würde gerne wissen, wie das geht und ob ich eine vernünftige Strategie verfolge.

Hier ist ein Ausrichtungsbeispiel: qcovs=100 aber qcovhsp niedriger.

qseqid           sseqid         pident length mismatch gapopen qstart   qend  sstart    send   evalue bitscore qcovs qcovhsp
HPNK_01698      HAPS_0519       81.88   596     75      5       630     1225    615     1177    0.0      889    100     49
HPNK_01698      HAPS_0519       49.17   301     115     8       84      366     201     481     2e-56    214    100     23
HPNK_01698      HAPS_0519       53.64   261     61      6       436     684     616     828     6e-49    191    100     20
HPNK_01698      HAPS_0519       46.61   251     62      3       332     510     584     834     6e-46    181    100     15
HPNK_01698      HAPS_0519       53.27   214     79      4       1       194     1       213     1e-45    180    100     16
HPNK_01698      HAPS_0519       55.96   218     60      8       550     764     643     827     1e-40    164    100     18
HPNK_01698      HAPS_0519       51.56   225     61      7       516     731     642     827     1e-38    157    100     18
HPNK_01698      HAPS_0519       49.57   230     77      6       484     713     643     833     1e-37    154    100     19
HPNK_01698      HAPS_0519       57.89   76      26      1       364     433     760     835     1e-13   76.3    100     6

Code verwendet

Datenbank machen

makeblastdb -in $Hparasuisfastadatabase -out H_parasuis_strains_gb_ALL.fna_databaseBLAST -dbtype prot -parse_seqids 

Führen Sie BLAST aus

blastp -db H_parasuis_strains_gb_ALL.fna_databaseBLAST -query 'out_2.fasta' -out HPNK_selected_vs_H_parasuis_strainss.tblastn -evalue 0.00001 -outfmt "6 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore qcovs qcovhsp" -max_target_seqs 50

Danke, Bernhard

Wenn ich mich nicht irre, ist qcovhsp eine perc_identity für EINEN HSP. Ich müsste manuell für jeden Treffer als Durchschnitt für alle HSPs perc_identity für einen Treffer berechnen.
Ich mag hit_perc_identity wegen seiner reduktiven Informationsannahme nicht einmal. Der Durchschnitt wird normalerweise von Extremwerten beeinflusst.
Welche Art vergleichst du? Wie weit sind sie? Wenn sie nahe genug sind, sollten "echte" Homologe in der Lage sein, ein einziges HSP zu bilden, das die meisten Zielsequenzen umfasst.
Sie können die durchschnittliche HSP-Abdeckung für die gesamte Sequenz (aus allen Alignments) berechnen. Sie können auch die durchschnittliche Punktzahl pro Rückstand berechnen. Die andere Möglichkeit besteht darin, eine End-to-End-Ausrichtung mit Werkzeugen wie Stretcher durchzuführen .
Schauen Sie sich meine neue Frage zu BLAST an: biology.stackexchange.com/questions/23958/…

Antworten (1)

Zuallererst, wenn Sie 90%Identität wollen, können Sie diesen Treffer verwerfen. Keiner der HSP überschreitet diese Schwelle. Da Sie mit Proteinen arbeiten, gibt es darüber hinaus keine Spleißprobleme, und Sie sollten in der Lage sein, ein einziges HSP zu erhalten, das die meisten Abfrage- und Subjektsequenzen umfasst. Vorausgesetzt natürlich, Sie haben ein echtes Homolog.

In Ihrer Ausgabe sehe ich viele kleine, sich überschneidende HSPs, von denen die meisten eine geringe Identität haben. Ich kann nicht sicher sein, ohne die Sequenz zu sehen, aber es ist eine sichere Wette, dass das, was Sie dort haben, Bereiche mit geringer Komplexität/sich wiederholende Bereiche sind, und das ist der Grund, warum Sie so viele separate HSP haben. Das nur halbwegs anständige beginnt an Position 630 der Abfragesequenz und ist nur 595 Reste lang, weniger als die Hälfte Ihres Abfrageproteins. Entweder haben Sie eine sehr divergierende N-terminale Region oder Ihr HSP ist nur eine konservierte Domäne. Auch hier müsste ich die tatsächliche Sequenzausrichtung sehen, um sicher zu sein, aber das sieht nicht wie ein echtes Homolog aus (vorausgesetzt, Ihre Arten sind sich ziemlich ähnlich, was sie sein müssen, wenn Sie eine Identitätsschwelle von 90% verwenden).

Also, immer vorausgesetzt, dass Ihre Spezies nahe genug ist, um anständige Homologe zu erwarten, würde ich einfach die kürzeren HSPs ignorieren und mich mit denen befassen, die mehr als sagen wir 80 % der Länge meiner Abfrage bei> = 90 % Identität darstellen. Kürzere Treffer werden meistens konservierte Domänen oder Regionen mit sich wiederholender/niedriger Komplexität sein. Die von Ihnen gewählten Schwellenwerte hängen von der Art ab, die Sie untersuchen.

Wenn Ihre Arten nicht so nah sind, verwenden Sie BLASTP überhaupt nicht. Stattdessen können Sie etwas wie verwenden hmmer. Sammeln Sie für jedes Ihrer Suchproteine ​​eine Reihe von Homologen aus verschiedenen Arten, erstellen Sie daraus eine Matrix und verwenden Sie diese Matrix, um Ihre Datenbank zu durchsuchen. Sie könnten auch Selenoprofiles verwenden, das einen ähnlichen Ansatz verwendet.

Sie sind Stämme derselben Bakterienart. Ich stimme den von Ihnen vorgeschlagenen Schwellenwerten zu. Vielen Dank