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.
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
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.
Terdon
Benutzer2030378