注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

Q超越兔子的蜗牛O--逸云沙鸥Linux

飘飘何所似,天地一沙鸥;落霞与孤鹜齐飞,秋水共长天一色~~

 
 
 

日志

 
 

第一个比较完整的perl程序  

2010-02-04 11:06:32|  分类: Perl编程技巧 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |
内容描述:从一系列.snp 文件中 读入各种深度(depth)的个数(count)信息,给出深度区间,以及一定比列(eg 总体0.95 范围内的 深度区间。)。实现了 参数输入,按需输出。
  1
#!/usr/bin/perl
  2 #===============================================================================
  3 #
  4 #         FILE:  snp_depth_dis_qunero.pl
  5 #
  6 #        USAGE:  ./snp_depth_dis_qunero.pl  [options]  <input *.snp file list>
  7 #
  8 #  DESCRIPTION:  
  9 #
 10 #      OPTIONS:  ---
 11 # REQUIREMENTS:  ---
 12 #         BUGS:  ---
 13 #        NOTES:  ---
 14 #       AUTHOR:  QuNengrong (Qunero), QuNengrong@gnomics.cn
 15 #      COMPANY:  BGI
 16 #      VERSION:  1.1
 17 #      CREATED:  2010年02月02日 12时37分29秒
 18 #     REVISION:  ---
 19 #===============================================================================
 20
 21 =head1 Name
 22
 23     snp_depth_dis_qunero.pl -- list loci count of different depths, and give the depth range whose loci count proportion reach a certain percentage
 24
 25 =head1 Usage
 26
 27     perl snp_depth_dis_qunero.pl [options]  <input *.snp file list>
 28     -h,--help   print this help to screen
 29     -o file     to specify a file to save the result ,default : [stat.info]
 30     -p proportion   to specify the intra_proportion   default : [0.95]  
 31     -s      save kind_depth_dis information ,format: ["kind"_depth_dis.txt]
 32     
 33 =head1 Example
 34
 35     perl snp_depth_dis_qunero.pl test.snp
 36     perl snp_depth_dis_qunero.pl -s -o result.txt test/*.snp
 37
 38 =head1 Version
 39     Verion  :   1.1
 40     Created :   2010年02月02日 12时37分29秒
 41     Updated :   Wed Feb  3 23:10:45 CST 2010
 42
 43
 44 =head1 Contact
 45     Author  :   QuNengrong (Qunero)
 46     E-mail  :   Quner612@qq.com,Quner8@gmail.com
 47     Company :   BGI
 48
 49
 50 =cut
 51
 52 use strict;
 53 use warnings;
 54 use Getopt::Long;
 55
 56 my ($need_help , $out_file , $intra_prop , $save_dis);
 57
 58 GetOptions(
 59     "help" => \$need_help,
 60     "o=s"   =>  \$out_file,
 61     "p=f" => \$intra_prop,
 62     "s" =>  \$save_dis,
 63 );
 64
 65 die `pod2text $0` if ($need_help);
 66
 67 #===============================================================================
 68 #                   Global Variable
 69 #===============================================================================
 70 $intra_prop ||= 0.95;                           # set default value 0.95
 71 $intra_prop = 0.95  unless (0 < $intra_prop and $intra_prop < 1);
 72 my @kind_depth_2_count_2Darr = ();              # kind_depth_2_count_2Darr[kind][depth],kind is in ALL and UNI
 73 my @output_kind_info_2Darr = ();                # output_kind_info_2Darr[kind][info]: ($sum ,$minD , $maxD , $min_intraD , $max_intraD );
 74 my ($minD , $maxD , $min_intraD , $max_intraD );
 75 my $kind;                                       # UNI ,ALL
 76 my $ALL = 1;
 77 my $UNI = 0;
 78 my $lo;                                         # var for loop
 79 my $lines_count = 0;
 80 $output_kind_info_2Darr[$UNI]->[0] = 0;
 81 $output_kind_info_2Darr[$ALL]->[0] = 0;
 82
 83
 84 #===============================================================================
 85 #                   Main process
 86 #===============================================================================
 87 if (0 == @ARGV){                                # die if no snp file as input
 88     print "$0 need at least one snp file to process\n";
 89     die "type $0 -h for more information";
 90 }
 91
 92 #step   001:    read files and count depth
 93 foreach my $one_snp_file ( @ARGV ){
 94     &count_depth($one_snp_file);
 95 }
 96
 97 #step   002:    compute minD,maxD,intra_minD,intramaxD
 98 &comput();
 99
100
101 #step   003:    output result to stdout or specified file
102 if(defined $out_file){
103     open STDOUT,">","$out_file";
104 }
105 &output_result();
106
107
108 #step   004:    save depth_dis file if necessary
109 if(defined $save_dis){
110     &save_depth_dis();
111 }
112 close STDOUT;
113 #===============================================================================
114 #               Subroutines
115 #===============================================================================
116
117 #001:   define sub count_depth
118 sub count_depth{
119     my ($one_snp_file) = @_;
120     if (defined $one_snp_file){
121         open SNP,"<",$one_snp_file or die "can't open $one_snp_file for read!";
122     }
123     my ($uni_once , $all_once ,$kind);
124     while(<SNP>){
125         my @split_line= split;
126         $lines_count += 1;
127         ($uni_once ,$all_once) = @split_line[7,8];
128         $output_kind_info_2Darr[$UNI]->[0] += 1;
129         $output_kind_info_2Darr[$ALL]->[0] += 1;
130         if(exists $kind_depth_2_count_2Darr[$UNI][$uni_once]){
131             $kind_depth_2_count_2Darr[$UNI][$uni_once] ++ ;
132             #print "test005:$uni_once\t$kind_depth_2_count_2Darr[$UNI][$uni_once] \n"; # test005
133         }else{
134             $kind_depth_2_count_2Darr[$UNI][$uni_once] = 1;
135             #print "test005:$uni_once\t$kind_depth_2_count_2Darr[$UNI][$uni_once] \n"; # test005
136         }
137         if(exists $kind_depth_2_count_2Darr[$ALL][$all_once]){
138             $kind_depth_2_count_2Darr[$ALL][$all_once] ++ ;
139             #print "test005_all :$all_once\t${kind_depth_2_count_2Darr[$ALL][$all_once]}\n";
140         }else{
141             $kind_depth_2_count_2Darr[$ALL][$all_once] = 1;
142             #print "test005_all :$all_once\t${kind_depth_2_count_2Darr[$ALL][$all_once]}\n";
143         }
144     }
145     close SNP;
146 }
147
148
149 #002:   define sub comput
150 sub comput{
151     #minD,maxD  
152     foreach  $kind ($UNI,$ALL){
153         $maxD = scalar @{$kind_depth_2_count_2Darr[$kind]} - 1; # found maxD **maxD is the offset ,i.e. maxD=total_of_arr - 1
154         #print "test006: $maxD\n";
155         $output_kind_info_2Darr[$kind]->[2] = $maxD;
156         for($lo = 0; $lo <= $maxD; $lo ++){      # found minD
157             if(exists $kind_depth_2_count_2Darr[$kind][$lo]){
158                 $output_kind_info_2Darr[$kind]->[1] = $lo;
159                 last;
160             }
161         }
162     }
163
164     #min_intraD,max_intraD
165     my ($cutoff ,$tmp_count)=(0 , 0);
166     foreach $kind ($UNI ,$ALL){
167         $tmp_count = 0;
168         $cutoff =  ($output_kind_info_2Darr[$kind]->[0] * (1 - $intra_prop) / 2);
169         #print "test007: cutoff = $cutoff\n";    # test007
170         for($lo = $output_kind_info_2Darr[$kind]->[1]; $lo <= $output_kind_info_2Darr[$kind]->[2]; $lo ++){
171             if(exists $kind_depth_2_count_2Darr[$kind][$lo]){
172                 $tmp_count += $kind_depth_2_count_2Darr[$kind][$lo];
173             }
174             if($tmp_count >= $cutoff){
175                 $output_kind_info_2Darr[$kind]->[3] = $lo;
176                 last;
177             }
178         }
179         $tmp_count = 0;
180         for($lo = $output_kind_info_2Darr[$kind]->[2]; $lo > $output_kind_info_2Darr[$kind]->[3]; $lo --){
181             if(exists $kind_depth_2_count_2Darr[$kind][$lo]){
182                 $tmp_count += $kind_depth_2_count_2Darr[$kind][$lo];
183                 #print "test008: lo=$lo ,tem_count = $tmp_count \n";
184             }
185             if($tmp_count >= $cutoff){
186                 $output_kind_info_2Darr[$kind]->[4] = $lo;
187                 last;
188             }
189         }
190     }
191 }
192
193 #003:   define sub output_result
194 sub output_result{
195     #write to file or STDOUT
196     print STDOUT "id(0-->UNI,1-->ALL)\ttotal_count\tmin_depth\tmax_depth\tintra_min_depth\tintra_max_depth\n";
197     foreach $kind ($UNI,$ALL){
198         print STDOUT "$kind\t$lines_count\t";
199         for($lo = 1; $lo <5; $lo++){
200             print STDOUT "$output_kind_info_2Darr[$kind]->[$lo]\t";
201         }
202         print STDOUT "\n";
203     }
204 }
205
206
207 #define sub save_depth_dis
208 sub save_depth_dis{
209     foreach $kind ($UNI , $ALL){
210         open OUT,">","${kind}_depth_dis.txt" or die "Can't open ${kind}_depth_dis.txt to save depth_dis";
211         for($lo = $output_kind_info_2Darr[$kind][1] ; $lo <= $output_kind_info_2Darr[$kind][2] ; $lo ++ ){
212             if(exists $kind_depth_2_count_2Darr[$kind][$lo]){
213                 print OUT "$lo\t$kind_depth_2_count_2Darr[$kind][$lo]\n";
214             }
215         }
216     }
217
218 }
219 
  评论这张
 
阅读(1493)| 评论(0)
推荐 转载

历史上的今天

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2017