perl怎么提取基因组所有基因的启动子序列
perl怎么提取基因组所有基因的启动子序列
这篇“perl怎么提取基因组所有基因的启动子序列”文章的知识点大部分人都不太理解，所以小编给大家总结了以下内容，内容详细，步骤清晰，具有一定的借鉴价值，希望大家阅读完这篇文章能有所收获，下面我们一起来看看这篇“perl怎么提取基因组所有基因的启动子序列”文章吧。
脚本运行命令：
perlgene_promoter.pl-faDonkey_Hic_genome.20180408.fa-gffDonkey_Hic_genome.20180408.gff3-outgene_promoter.fa-n2000
其中 -fa 后跟基因组染色体序列；-gff 后跟基因组gff文件；-n后跟数字，表示要提取基因上游多少bp的序列。
脚本代码：
#!/usr/bin/perl-wusestrict;usewarnings;useGetopt::Long;useData::Dumper;useConfig::General;useCwdqw(abs_pathgetcwd);useFindBinqw($Bin$Script);useFile::Basenameqw(basenamedirname);useBio::SeqIO;useBio::Seq;my$version="1.3";##prepareparameters#########################################################################-------------------------------------------------------------------------------------------##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<<"UsageEnd.";Description:$version:lefseanalysisUsageForcedparameter:-gffgfffile<infile>mustbegiven-outoutdir<outfile>mustbegiven-nnum<int>-fagenomefastafile<infile>mustbegivenOtherparameter:-hHelpdocumentUsageEnd.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"openfile$opts{gff}faild.\n";open(OUT,">$opts{out}")||die"openfile$opts{out}faild.\n";while(<IN>){nextif(/^#/);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);printOUT">$name\n$gene\n";}elsif($line[6]eq"-"){my$gene=substr($fasta{$line[0]},$line[4],$n);$gene=&reverse_complement_IUPAC($gene);printOUT">$name\n$gene\n";}}}close(OUT);close(IN);subreverse_complement_IUPAC{my$dna=shift;#reversetheDNAsequencemy$revcomp=reverse($dna);#complementthereversedDNAsequence$revcomp=~tr/ABCDGHMNRSTUVWXYabcdghmnrstuvwxy/TVGHCDKNYSAABWXRtvghcdknysaabwxr/;return$revcomp;}subreverse_complement{my$dna=shift;#reversetheDNAsequencemy$revcomp=reverse($dna);#complementthereversedDNAsequence$revcomp=~tr/ACGTacgt/TGCAtgca/;return$revcomp;}
以上就是关于“perl怎么提取基因组所有基因的启动子序列”这篇文章的内容，相信大家都有了一定的了解，希望小编分享的内容对大家有帮助，若想了解更多相关的知识内容，请关注恰卡编程网行业资讯频道。
