搜索
查看: 3249|回复: 1

【bioamin-perl练习1】注释文件.gff和序列文件.fna提取16srRNA

[复制链接]

2

主题

10

帖子

65

积分

注册会员

Rank: 2

积分
65
发表于 2017-7-26 15:45:16 | 显示全部楼层 |阅读模式
本帖最后由 bioamin 于 2017-7-28 20:45 编辑

代码比较烂,欢迎各位批评指正;
1.先了解注释文件.gff的文件格式:.gff文件共有9列,以\t为分隔符,每列含义如下:1.序列的编号2.注释信息的来源3.注释信息类型4基因开始位置5.基因结束位置6.得分7.序列方向8.“phase”9.注释信息描述;本次提取信息使用第4.5.9列;(详细.gff信息请进入传送门http://boyun.sh.cn/bio/?p=1602
2.了解.fna文件格式:.fna文件相当于只有一条信息的基因文件,第一部分为编号等信息,第二部分为序列信息;
>NZ_CP010106.1 Bacillus thuringiensis serovar indiana strain HD521, complete genome
TAGCCACTTTTTTTTGATATTATAGTTGTGTTTTCACTTTGAATAAGTTTTCCACACCTTTATTTTATCCACAGTTTGTG
TATAACATGTGGACAGTTTTAATCACATGTGGGTAAATAGTTGTCCACATTTGTTTTTTTTGTCGAAAACCCTATCTCAT
ATACAAACGACGTTTTTAGGTTTTAAAATACGTTTCGTATAAATATACATTTTATATTTATTGGGTTGTACATTTGTTGC
GCAACCTTATTCTTTTACCATCTTAGTAAAGGAGGGACACCTTTGGAAAATATCTCTGATTTATGGAATAGTGCCTTAAA
AGAATTAGAAAAAAAGGTAAGCAAGCCTAGTTATGAGACATGGTTAAAATCAACAACGGCTCATAACTTGAAGAAAGACG
TATTAACGATTACAGCTCCAAATGAATTTGCTCGTGACTGGCTAGAATCTCATTACTCCGAACTAATTTCAGAAACACTA
TACGATTTAACAGGGGCAAAATTAGCAATTCGCTTTATTATTCCCCAAAGTCAAGCTGAAGAGGATATTGATCTGCCTCC


3.代码思路:首先输入.fna文件将序列信息整合到一个变量中;其次输入.gff文件,以\t为分隔符赋予数组,提取出起始、终止、和注释信息;最后对序列信息进行切割,输出。
open(IN1,"genomic.gff") or die ("can not open gff file");
open(IN2,"genomic.fna") or die ("can not open fna file");
open (out,">out");
$gene;
%hash;
while (<IN2>)#按行输入.fna
{   if(/^>/){next;}
    else
   {chomp($_);$gene.=$_;#将.fna文件里的序列整合为一条,并赋予变量$gene;
    }
}

while (<IN1>)#.gff按行输入
{ chomp($_);
   @arr=split("\t",$_);#以tab为分隔符切割
   $name=$arr[8];$begin=$arr[3];$end=$arr[4];$len=$end-$begin+1;#rna开始为第4列,终止为第5列,计算长度,注释信息第9列;
       if($name=~m/product=16S ribosomal RNA/i)#如果注释信息匹配到16SRNA,输入匹配到31个
         {
            $result=substr($gene,$begin,$len);     #以开始位置,长度对染色体序列进行切割,切割结果放在变量$result;
            $hash{$begin}=$name."\t".$result;   #由于存在重复(起始,终止都一样),所以起始位置为键,“注释信息\t切割的序列结果”为values构建哈希
          }
  }
while (($key,$value)=each %hash)#遍历hash
{ @val=split("\t",$value);  #以\t为分割符,对哈希的内容进行分割
  print out ">"."$val[0]\n";  #将注释信息输出到文件out
  print out "$val[1]\n";     #将切割结果的序列信息输出到文件out
}
close IN1;close IN2;close out;




上一篇:从UCSC上下载人类基因组TSS(转录起始位点)区域bed文件实战
下一篇:为什么reads的每个位置的ATCG碱基理论上比例是一样的?
回复

使用道具 举报

2

主题

10

帖子

65

积分

注册会员

Rank: 2

积分
65
 楼主| 发表于 2017-7-26 15:50:09 | 显示全部楼层
师兄一行确定16S RNA有多少条
less /home/genomic.gff|grep 'product=16S ribosomal RNA'|awk '{print $4,$5}'|sort -u|wc -l
输出13
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-10-14 22:20 , Processed in 0.038317 second(s), 25 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.