tBLASTn can`t find proteins annotated in same genome. Is that possible?
Seems, tBLASTn can't find some of the proteins that were annotated using PGAAP pipeline. Here is the code and the results:
# download Nagasaki genome file
wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA_000600995.1_ASM60099v1/GCA_000600995.1_ASM60099v1_genomic.fna.gz
mv GCA_000600995.1_ASM60099v1_genomic.fna Nagasaki.fna
# download Nagasaki gbk file
wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA_000600995.1_ASM60099v1/GCA_000600995.1_ASM60099v1_genomic.gbff.gz
mv GCA_000600995.1_ASM60099v1_genomic.gbff Nagasaki.gbk
# extract selected proteins from gbk file
python /Users/bernardo/Documents/BioLinux/A0_scripts/parse_gbk.py Nagasaki.gbk > Nagasaki.faa
# count proteins
grep '>' Nagasaki.faa | wc -l
> 2260
# BLASTp proteins against proteome
makeblastdb -in Nagasaki.faa -out H_parasuis_strains_gb_ALL.fna_databaseBLAST -dbtype prot -parse_seqids # blast database
blastp -db H_parasuis_strains_gb_ALL.fna_databaseBLAST -query 'Nagasaki.faa' -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" -max_target_seqs 50
awk '$3 > 99 { print $0 }' HPNK_selected_vs_H_parasuis_strainss.tblastn | awk '$13 > 99 { print $0 }' | wc -l # count sequences with more than 60% identities and 70% query coverage per HSP
> 2545
awk '$3 > 99 { print $0 }' HPNK_selected_vs_H_parasuis_strainss.tblastn | awk '$13 > 99 { print $0 }' > HPNK_selected_vs_H_parasuis_strainss.tblastn.tab
echo 'see the number of queries that had HIT'
cat HPNK_selected_vs_H_parasuis_strainss.tblastn.tab | awk '{print $1}' | sort | uniq -c | sort | wc -l # see the number of queries that had HIT
> 2260
# tBLASTn proteins against genome
# BLASTp proteins against proteome
makeblastdb -in Nagasaki.fna -out H_parasuis_strains_gb_ALL.fna_databaseBLAST -dbtype nucl -parse_seqids # blast database
tblastn -db H_parasuis_strains_gb_ALL.fna_databaseBLAST -query 'Nagasaki.faa' -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" -max_target_seqs 50
awk '$3 > 99 { print $0 }' HPNK_selected_vs_H_parasuis_strainss.tblastn | awk '$13 > 99 { print $0 }' | wc -l # count sequences with more than 60% identities and 70% query coverage per HSP
> 2509
awk '$3 > 99 { print $0 }' HPNK_selected_vs_H_parasuis_strainss.tblastn | awk '$13 > 99 { print $0 }' > HPNK_selected_vs_H_parasuis_strainss.tblastn.tab
echo 'see the number of queries that had HIT'
cat HPNK_selected_vs_H_parasuis_strainss.tblastn.tab | awk '{print $1}' | sort | uniq -c | sort | wc -l # see the number of queries that had HIT
> 2108
Thanks, Bernardo
• 324 views
•
link
0 answers
No answers yet.
Log in to answer this question.
Seems there is a problem is BLAST output parsing step, revising it now.