Hi Everyone, I have been given a task to compare one gene sequence among 50 strain of E.coli. For this study i have 50 genome scaffold file and one gene sequence file. Now my work will be compare the gene sequence among all the genome and compute a phylogenetic tree of the gene among all sequence. If anybody could point me in the right direction, I would be thankful!.
2 answers
As far as I understand your question, for creating of that tree, you'll need to extract 50 sequences, one from each file representing a strain genome. Once you'll have these 50, you will need to align all them with clastalw and then make a tree. To get 50 sequences you can use blast - in this case i would build a db for each of 50 files and blast the gene against them, then pick a sequence which "seems like a gene i am looking for". This part - picking a sequence that looks like a gene - can be a bit tricky if you have only scaffolds. Another way to get 50 sequences would be to run prodigal on each of scaffolds fasta files getting sets of ORFs and then blast the gene sequence against predicted ORFs.
sounds good. but, as i have mentioned, tblastn will give you a set of hits, which you'll need to somehow process in order to extract gene sequences - i.e eyeball and pick the best hit manually assigning start stop _if scaffolds are of good quality_... but what if there are gaps and mis-assemblies - what will be that best hit? using gene prediction may be better option: you run prodigal which will yield you protein sequences, then blastp those with your gene - pick a best hit putative gene sequence from each of 50 and proceed with the tree construction.
# Predicted gene
for file in ./*.fasta
do
echo $file
./prodigal.v2_60.linux -a out_trans -d out_nuc -i $file
mv out_trans $file_trans
mv out_nuc $file_nuc
done
Hi i am trying to run all 50 genome like that but could you please check why it is not rename the out put file according to input file name.
shouldn't it be "$file"_trans and "$file"_nuc?
There's a pipeline for this. It's called hal. However, in my opinion, you get more reliable trees when you only include conserved genes that are unlikely to transfer horizontally (e.g. ribosomal proteins). Also, the super alignments hal constructs from bacterial genomes are ridiculously long (100k-200k aa), so if you go this route you can only make ML and MP trees (no Bayesian).
that is a cool stuff, thanks!
I have extracted all the translated gene sequence from all 50 genome using Prodigal. Now i have to run a bath mode blast of all the sequence vs my query sequence and extract the max score gene and finally generate the tree. Could you please suggest me something , how to proceed next step. Although i had a look "hal" as you mention i want to try both the option.
Log in to answer this question.
-
How can I confidently determine if my genome sequence Is truly Anaplasma marginale?
written by ashiqullah024 •
-
Chloroplast genome phylogenetic analysis
written by a.bibek52 •
-
Using BioPerl to create Phylogenetic Tree
written by Andrea •
-
How To Find Core Genome Snp
written by HG
-
How To Extract Max Score Blast Hit Among A Huge Data Set?
written by HG
-
How To Extract Common Cds Among Huge Data Set?
written by HG
-
How To Extract Core Genome From Whole Genome Sequence Data?
written by HG
-
How To Calculate Coverage
written by HG
-
How To Find Suitable Reference Genome After Genome Assembly
written by HG
-
Compare Two Protein Sequences Using Local Blast
written by satsurae •
