Hello everyone,
I am trying to determine the abundance of a gene in a set of metagenomic short reads using the hmm of this gene, then use this method to compare the abundance of the gene across several samples after normalizing for data size. What I would do is to first translate the genes into ORFs, then run hmm against the reads. Since the amino acid sequence of this HMM is longer (~220 amino acids) than the read length (150 bases, which translates to ~50 amino acids), I am wondering if hmm can accurately detect the partial sequences and still report hits?
Thank you!
Rui
1 answer
If you use HMMer, the answer is yes - unless you change the default glocal scoring. However, there is going to be less signal in 50 residues compared to the whole protein, so you may miss some divergent sequences. Also, some reads will have less than 49-50 residues if they are at either end of the protein, which may also leave them out.
Why not assemble instead and determine the abundance properly?
Thanks for the comments! I was just curious about how would the results differ between reads and assembly based analysis. I also assembled the reads with MEGAHIT, after retaining the contigs longer than 500 bases, I remapped the reads back onto the contigs with bowtie2, but the alignment rate is only about 30% (probably because my samples are derived from soil). So I just thought only using the contigs would result in a large chunk of the data being wasted. However, the problem is even the reads cannot be guaranteed to capture the whole picture of the functional diversities in the environment. I think using assemblies would be a better idea in my case.
If you still care about the reads not included in the final assembly you might want have a look at the raw reads analysis pipeline v5.0 from MGnify. It does not involve any assembly step and it is specifically designed for the functional annotation of raw sequencing data
Thank you for the suggestion, I will take a look into it! I've tried similar tools before such as HUMANN3 that is based on reads data for functional annotation, but the results were not satisfactory, ~90% of the reads cannot be classified. While looking for a solution, I found this paper interesting as it benchmarked various analysis approaches. It seems assemblies are preferred over reads when the alignment rates are high.
Log in to answer this question.
-
Identify unique amino acids in a seuence of interest present in a MSA block
written by sushmitapaul007 •
-
How to Detect Amino Acid Changes from a GATK VCF File and Reference Genome Using Code?
written by DareDevil
-
Translating a Fasta File
written by Hb •
-
Compare protein functionality between species
written by Howard •
-
How can I retrieve nucleotide sequences from amino acid sequences?
written by Rui •
-
How to read force field parameters of amino acids such as vanderwaal radii and partial charge?
written by physicsvigna •
-
Nucleotide to amino acid conversion of multiple sequences
written by lokraj2003
-
Determinig the number of hidden states in a Hidden Markov Model
written by Tamador •
-
Effect of indel mutation on protein
written by sahar.ava •
-
Theta I And Theta Ii Functional Divergence In Proteins
written by Rui •
