@Fatima, great - that should work. Thank you :)
• 0 views
•
link
written
by
jamie.pike
Is there a way to produce/convert nhmmer output to a bed file?
I have recently used nhmmer to search for a sequence in a FASTA file, I'd like to compare the output of the nhmmer search to a bed file I have and think it would be easy to do if the nhmmer output is also in bed format.
Is there a tool to convert the nhmmer output? I have had a look in the nhmmer manual and can't find any specific flags for generating either a bed file output or a table output indicating locations of the hits.
Example of nhmmer output:
Annotation for each hit (and alignments):
>> VMNF01000014.1_Fusarium_oxysporum_f._sp._cubense_strain_TR4_isolate_UK0001_scf_28419_14,_whole_genome_shotgun_sequence
score bias Evalue hmmfrom hmm to alifrom ali to envfrom env to sq len acc
------ ----- --------- ------- ------- --------- --------- --------- --------- --------- ----
! 134.7 1.6 3.9e-39 1 206 [] 3509332 3509126 .. 3509332 3509126 .. 3744637 0.99
Alignment:
score: 134.7 bits
All_mimps_comb_trimmed.aln 1 AcagTgggatGcaAtaAGtttgaatacgtgtgtaagTagg 40
AcagTgggatGcaA+aAGt+t++++a+gtgt taagT
VMNF01000014.1_Fusarium_oxysporum_f._sp._cubense_strain_TR4_isolate_UK0001_scf_28419_14,_whole_genome_shotgun_sequence 3509332 ACAGTGGGATGCAAAAAGTATTCGCAGGTGTATAAGTC-- 3509295
79***********************************8.. PP
All_mimps_comb_trimmed.aln 41 ttcccaagctagcctagtta.tggggtaccttagcttagg 79
+ccc+agct +ctagtt t+g ta+ctt+ cttagg
VMNF01000014.1_Fusarium_oxysporum_f._sp._cubense_strain_TR4_isolate_UK0001_scf_28419_14,_whole_genome_shotgun_sequence 3509294 -GCCCTAGCTCCTCTAGTTCgTAGCTTAGCTTGACTTAGG 3509256
.*************************************** PP
All_mimps_comb_trimmed.aln 80 gctttgctacttagatAtcacactgaaattagtataatat 119
+ +ttg+tact +g+tAtc c ctga+ t agt+taatat
VMNF01000014.1_Fusarium_oxysporum_f._sp._cubense_strain_TR4_isolate_UK0001_scf_28419_14,_whole_genome_shotgun_sequence 3509255 TACTTGATACTAGGTTATCCCCCTGATCTGAGTGTAATAT 3509216
**************************************** PP
All_mimps_comb_trimmed.aln 120 taaggtaTcgagtacctaaaaaa.aactcctagcTagcta 158
aaggta +gag acct+a aa+ ac+cct++cTagct
VMNF01000014.1_Fusarium_oxysporum_f._sp._cubense_strain_TR4_isolate_UK0001_scf_28419_14,_whole_genome_shotgun_sequence 3509215 -AAGGTACTGAGGACCTGACAAGcTACCCCTGACTAGCTT 3509177
.8************************************** PP
All_mimps_comb_trimmed.aln 159 gcgagtgtc...ggtacttcaacacgtattcatacttatt 195
cgag++tc tac+t acac+t++++atactt+tt
VMNF01000014.1_Fusarium_oxysporum_f._sp._cubense_strain_TR4_isolate_UK0001_scf_28419_14,_whole_genome_shotgun_sequence 3509176 CCGAGCATCagaCCTACCTGCACACCTGCGAATACTTTTT 3509137
************99************************** PP
All_mimps_comb_trimmed.aln 196 gCatcccaCtg 206
gCatcccaCtg
VMNF01000014.1_Fusarium_oxysporum_f._sp._cubense_strain_TR4_isolate_UK0001_scf_28419_14,_whole_genome_shotgun_sequence 3509136 GCATCCCACTG 3509126
*********96 PP
Thanks,
Jamie
• 3,173 views
•
link
updated
by
colindaven
•
written
by
jamie.pike
•
written
by
jamie.pike
2 answers
I think if you use the commandline version, and add --noali option, that would be easier.
Then you can use awk to extract the columns that you need:
awk '($1!~"#"){print $1,$2,$3}'
Just change 1,2,3 to columns that you're interested in.
I assume you need to add strand information too.
@Fatima, Strand information would be useful also but I can't work out where it is?
• 0 views
•
link
written
by
jamie.pike
I assume since you're aligning the whole genome, it would be + for all the entries.
awk '($1!~"#"){print "VMNF01000014.1\t"$7"\t"$8"\t"$1"\t.\t+"}'
. for score column
I usually get the genes in the genome, and then use hmmer. I haven't used nhmmer.
@Fatima, thank you. I'm trying to find a transposable element on specific scaffolds so there is variation in strandedness. Perhaps I will try annotating first using hmmer. Thank you for your help.
Ps I have found also that --tblout will provide location information and strandedness.
• 0 views
•
link
written
by
jamie.pike
No, I haven't heard of it! I'll take a look at it. Thank you :)
• 0 views
•
link
written
by
jamie.pike
To make it a bit clearer for nhmmer specifically, I used the following. This is code used in nextflow.
# run nhmmer
nhmmer --cpu $threads --noali --dna --tblout ${prefix}_TIR_out.tsv -o ${prefix}_TIR_out.txt $params.EG_TIRhmm $fasta
sleep 5
# convert to bed
awk '(\$1!~"#"){print \$1,"\t",\$7,"\t",\$8,"\t",".","\t",".","\t",\$12}' ${prefix}_TIR_out.tsv > ${prefix}_nhmmer_TIR_out.bed
• 0 views
•
link
written
by
colindaven
Log in to answer this question.