Lokale BLAST-Kopienzahl pro Treffer

Ich habe eine Reihe lokaler BLAST-Datenbanken mithilfe von makeblastdb aus metagenomischen Daten erstellt und suche nach dem Vorhandensein eines bestimmten Gens. Während ich die normale BLAST-Analyse durchführen kann, indem ich e-Werte, % Identität usw. betrachte, konnte ich nicht herausfinden, wie ich bestimmen kann, wie oft eine bestimmte Sequenz in der Datenbank vorhanden ist.

Also: Wie bestimme ich bei einer BLAST-Datenbank mit metagenomischen Daten, in der ein bestimmtes Gen möglicherweise mehrfach vorhanden ist, die Kopienzahl eines bestimmten BLAST-Treffers?

Danke im Voraus.

Das hängt davon ab, wie Sie einen Treffer definieren. Sprechen wir über homologe Gene verschiedener Arten? Mehrere, leicht unterschiedliche Kopien desselben Gens in einem einzigen Genom? Mehrere identische Kopien eines einzelnen Gens, das sich in der Datenbank wiederholt?
Ich definiere einen Treffer basierend auf den BLAST-Ergebnissen – ich habe seq1 gegen die Datenbank gesprengt und es hat seq2 ... seqn zurückgegeben. Abhängig davon, welche Cutoffs ich für Bit-Scores / Identitäts-Scores usw. verwende, habe ich möglicherweise drei "Hits" - seq2, seq3 und seq4. Ich möchte wissen, wie oft seq2 in den metagenomischen Daten auftauchte.

Antworten (2)

  • Definieren Sie einen "Hit" (basierend auf einem bestimmten Cutoff-Wert, einer Punktzahl usw.)
  • Erhalten Sie eine Ausgabe im Tabellenformat
  • Zählen Sie die Anzahl der Treffer pro Abfrage – wird normalerweise in der Kopfzeile angegeben; wenn Sie nach einigen ausgewählten Treffern suchen möchten (basierend auf einigen Cutoffs, dann können Sie die Datei parsen und herausfinden)

Beispieldatei (Header):

# BLASTN 2.2.27+
# Query: TCONS_00036712 gene=XLOC_017996
# Database: ../nt_db/nt
# Fields: query id, subject id, % identity, alignment length, mismatches, gap opens, q. start, q. end, s. start, s. end, evalue, bit score
# 1014 hits found

Wenn Sie den Header analysieren möchten, um nach der Anzahl der Treffer zu suchen, können Sie diesen Befehl im Terminal ausführen (wenn Sie awk haben – dort standardmäßig in Linux und anderen Unix-basierten Systemen. Für Windows installieren Sie gnuwin32).

awk -v OFS="\t" '/^# Query/{q=$3 FS $4} /^# .* hits found/{h[q]=$2} END{for(i in h){print i,h[i]}}' blastoutput.txt

Um die Anzahl der Treffer für ein bestimmtes Kriterium zu erhalten (z. B. habe ich gerade Bitscore [12. Spalte] als> 400 definiert)

awk -F "\t" '!/^#/ && $12>400{a[$1]++} END{for(i in a){print i,a[i]}}' blastoutput.txt
Großartig, aber das habe ich schon -- ich weiß, wie viele "Treffer", dh eindeutige Sequenzen, es gibt. Was ich wissen möchte, ist, wie oft die Sequenz dieses bestimmten Treffers in meinem zusammengestellten metagenomischen Datensatz vorhanden ist.
@user2030378 Möchten Sie nach redundanten Sequenzen suchen?

BLAST ist nicht das richtige Programm für metagenomische Analysen, um die Anzahl der Reads zu bestimmen, die auf eine bestimmte Region des Genoms abbilden, da es auf einzelne Treffer statt auf eine große Anzahl verschiedener Treffer optimiert ist.

Die schiere Anzahl von Treffern in einer Metagenomik-Datenbank (die eine große Menge redundanter Daten enthält) führt wahrscheinlich zu einer ineffizienten Skalierung und einer hohen Speichernutzung, wenn die Erkennung aller möglichen Treffer erforderlich ist.

Wenn die Anzahl der erforderlichen Treffer niedrig ist, ist es möglich, diese Treffer zurückzugeben, aber wenn Sie >1000 Treffer betrachten, wird wahrscheinlich empfohlen, ein Mapper-Programm wie BWA oder Bowtie zu verwenden . Sie sind ausschließlich für den Zweck konzipiert, Reads (z. B. aus einem RNASeq-Experiment) auf Genomen abzubilden.