zoukankan      html  css  js  c++  java
  • 基因组处理

        去除基因组序列中的未定位的scaffold、Contig序列和线粒体序,主要针对NCBI提供refseq基因组序列,组装到染色体级别的物种基本都通用。将所有碱基统一成大写字母,并计算每条染色体长度,每80个字符换行。

    处理脚本如下:

    image

      1 use strict;
      2 open A,"$ARGV[0]";
      3 open B,">$ARGV[1]";
      4 open C,">$ARGV[2]";
      5 my $help=<<USAGE;
      6 Usage: perl $0 genome.fa new.fa chrlen.list
      7 
      8 USAGE
      9 die "$help",unless(@ARGV==3);
     10 
     11 $/=">";
     12 <A>;
     13 my %chrlen;
     14 while(<A>){
     15 	chomp;
     16 	my @line=split /
    +/,$_;
     17 	my $seqName=shift @line;
     18 	my $chr=(split /s+/,((split /,/,$seqName)[0]))[-1];
     19 	next if $chr=~ /scaffold/;
     20 	next if $chr=~ /Contig/;
     21 	next if $chr=~ /mitochondrion/;
     22 	$chr="chr".$chr;
     23 	my $seq=join "",@line;
     24 	$seq=~s/
    //g;
     25 	$seq=uc($seq);
     26 	my $len=length($seq);
     27 	$chrlen{$chr}=$len;
     28 	$seq=~ s/(w{80})/$1
    /g;
     29 	if($len % 80 == 0){
     30 		print B ">$chr
    $seq";
     31 	}
     32 	else{
     33 		print B ">$chr
    $seq
    ";
     34 	}
     35 	print C "$chr	$chrlen{$chr}
    ";
     36 }
    本文版权归作者和博客园共有,欢迎转载,但未经作者同意必须保留此段声明,且在文章页面明显位置给出原文连接,否则保留追究法律责任的权利.
  • 相关阅读:
    [COCI2013]DLAKAVAC
    [TJOI2013]最长上升子序列
    AGC040E Prefix Suffix Addition
    AGC010E Rearranging
    AGC021F Trinity
    AGC002F Leftmost Ball
    JOISC2019D ふたつのアンテナ
    LOJ6210 「美团 CodeM 决赛」tree
    Luogu P3781 [SDOI2017]切树游戏
    Problem. M
  • 原文地址:https://www.cnblogs.com/mmtinfo/p/12165080.html
Copyright © 2011-2022 走看看