perl怎么提取基因组所有基因的启动子序列

这篇"perl怎么提取基因组所有基因的启动子序列"文章的知识点大部分人都不太理解，所以小编给大家总结了以下内容，内容详细，步骤清晰，具有一定的借鉴价值，希望大家阅读完这篇文章能有所收获，下面我们一起来看看这篇"perl怎么提取基因组所有基因的启动子序列"文章吧。

脚本运行命令：

perl gene_promoter.pl  -fa  Donkey_Hic_genome.20180408.fa  -gff  Donkey_Hic_genome.20180408.gff3  -out  gene_promoter.fa -n 2000

其中 -fa 后跟基因组染色体序列；-gff 后跟基因组gff文件；-n后跟数字，表示要提取基因上游多少bp的序列。

脚本代码：

#!/usr/bin/perl -wuse strict;use warnings;use Getopt::Long;use Data::Dumper;use Config::General;use Cwd qw(abs_path getcwd);use FindBin qw($Bin $Script);use File::Basename qw(basename dirname);use Bio::SeqIO;use Bio::Seq;my $version = "1.3";## prepare parameters ######################################################################### -------------------------------------------------------------------------------------------## GetOptionsmy %opts;GetOptions(\%opts,  "gff=s","fa=s", "out=s", "n=s","h");if(!defined($opts{out}) || !defined($opts{gff}) || || !defined($opts{fa}) ||defined($opts{h})){print <<"Usage End.";Description:$version:lefse analysisUsageForced parameter:-gff         gff file                            must be given-out          outdir                                must be given-n        num            -fa        genome fasta file        must be givenOther parameter:-h           Help documentUsage End.exit;}$opts{n} ||= 2000;my $n = $opts{n};my $in  = Bio::SeqIO->new(-file => "$opts{fa}" , -format => 'Fasta');my %fasta;while ( my $seq = $in->next_seq() ) {my($id,$sequence)=($seq->id,$seq->seq);$fasta{$id}=$sequence;}open(IN,"$opts{gff}") ||die "open file $opts{gff} faild.\n";open(OUT,">$opts{out}") ||die "open file $opts{out} faild.\n";while(){next if(/^#/);my @line = split ("\t",$_);if($line[2] eq "gene"){$line[8] =~ /ID=([^;]*);/;my $name = $1;if($line[6] eq "+"){my $gene = substr( $fasta{$line[0]},$line[3]-$n-1, $n);print OUT ">$name\n$gene\n";}elsif($line[6] eq "-"){my $gene = substr( $fasta{$line[0]},$line[4], $n);$gene = &reverse_complement_IUPAC($gene);print OUT ">$name\n$gene\n";}}}close(OUT);close(IN);sub reverse_complement_IUPAC {        my $dna = shift;        # reverse the DNA sequence        my $revcomp = reverse($dna);        # complement the reversed DNA sequence        $revcomp =~ tr/ABCDGHMNRSTUVWXYabcdghmnrstuvwxy/TVGHCDKNYSAABWXRtvghcdknysaabwxr/;        return $revcomp;}sub reverse_complement {        my $dna = shift;        # reverse the DNA sequence        my $revcomp = reverse($dna);        # complement the reversed DNA sequence        $revcomp =~ tr/ACGTacgt/TGCAtgca/;        return $revcomp;}

以上就是关于"perl怎么提取基因组所有基因的启动子序列"这篇文章的内容，相信大家都有了一定的了解，希望小编分享的内容对大家有帮助，若想了解更多相关的知识内容，请关注行业资讯频道。