How to extract domain sequence of merging Gene Family by perl
This article is about how perl extracts domain sequences from merged gene families. The editor thinks it is very practical, so share it with you as a reference and follow the editor to have a look.
In gene family analysis, after searching for a domain using hmmsearch, if the gene has more than one domain and all the domain extracts need to be linked together, you can use the following script.
How to use it:
The command is as follows:
Perl domain_hebing.fa.pl hmmsearch.out.txt Arabidopsis_thaliana.TAIR10.31.pep.all.fa domain.fa 0.001
Parameter description:
Hmmsearch.out.txt: the result file of the hmmsearch search.
Arabidopsis_thaliana.TAIR10.31.pep.all.fa: all gene protein sequences.
Domain.fa: outputs the merged domain sequence file.
0.001: sets the e-value value of the filter.
The script code is as follows:
Die "perl $0" unless (@ ARGV = = 4); use Math::BigFloat;use Bio::SeqIO;use Bio::Seq;$in = Bio::SeqIO- > new (- file = > "$ARGV [1]",-format = > 'Fasta'); $out = Bio::SeqIO- > new (- file = > "> $ARGV [2]",-format = >' Fasta'); my% fasta While (my $seq = $in- > next_seq ()) {my ($id, $sequence, $desc) = ($seq- > id, $seq- > seq, $seq- > desc); $fasta {$id} = $seq;} $in- > close (); my% keep = (); open IN, "$ARGV [0]" or die "$!"; while () {chomp;next if / ^ # /; my @ a = split /\ shands next if $a [6] > $ARGV [3] My $subseq = $fasta {$a [0]}-> subseq ($a [17], $a [18]); if (exists $keep {$a [0]}) {$keep {$a [0]}. = $subseq;} else {$keep {$a [0]} = $subseq;}} close (IN); while (my ($key,$value) = each% keep) {my $newseqobj = Bio::Seq- > new (- seq = > $value,-id = > $key,); $out- > write_seq ($newseqobj);} $out- > close (); Thank you for reading! This is the end of the article on "how to extract the domain sequence of the merged gene family from perl". I hope the above content can be of some help to you, so that you can learn more knowledge. if you think the article is good, you can share it out for more people to see!
