生物信息 perl 脚本实战
索引
1.统计fasta、fa和fastq文件的长度，统计fastq的reads个数，单个reads长度，reads总长度；统计fasta文件中contig的个数，列出名称，单条的长度，以及总长度。
2.1局部组装：创建目录，将比对好的reads按100k为单位，用samtools切，并用awk工具提起reads，分别存放在对应文件夹内
2.2局部组装：用得到的reads_name，去原始的下机reads里提取reads序列
2.3局部组装：使用canu/soapdenovo按窗口组装
基本套路：
1.使用传参数模块
use Getopt::Long; my ($first, $second, $third, $fourth); GetOptions( "first=s" => \$first, "second=s" => \$second, "third=s" => \$third, "fourth=s" => \$fourth, ) or die $!; print "$first-$second-$third-$fourth\n";
模块名是Getopt::Long，不能写错；使用时，得用GetOptions函数，函数中间是个哈希结构，将--first传给$first变量，以此类推。
这个模块有点复杂，具体参照：在Perl中使用Getopt::Long模块来接收用户命令行参数
2.读表格式文件，切分，取出目的字段，建立数据结构
open IN, "test.txt" or die $!; while(<IN>){ chomp; @list = split / /; next if $list[0] =~ /\bchar1\b|char4|char6/; $name = $list[2]; $value = $list[1]; $hash{$name} = $value; foreach (keys %hash){ print "$_ => $hash{$_}\n"; } } close(IN);
3.将原始的下机reads读到哈希结构里，reads名字作为键，其碱基序列作为哈希的值
#Store the original fasta file into a hash named %hash_reads. my %hash_reads; my $name; open IN,"<",$original_fasta or die $!; while(<IN>){ chomp; if(/^>(.*)$/){ $name=$1; $hash_reads{$name}=""; }else{ $hash_reads{$name} .= $_; } } close IN;
1.统计fasta、fa和fastq文件的长度，统计fastq的reads个数，单个reads长度，reads总长度（主要是统计总长度，其他在Linux下很简单就实现了）；统计fasta文件中contig的个数，列出名称，单条的长度，以及总长度。
思路整理：这是个典型的逐行读文件，取字段，计算长度问题。
fastq简单：四行一轮，解决办法很多，可以逐行读取，逐行匹配；也可以一次读两行；输出也少，单reads长度，reads数量，reads长度总和，也没有其他有价值的信息。
fasta略微复杂：没有规律，因为序列被切成短的，只能逐行读，逐行匹配；逐行读有个问题，怎么检测结束？（这是逐行读的一个最大缺陷，你无法操纵最后一次！！！）因此，只能把最后一次放在读循环的外面，每次输出的点就在匹配title那个地方。
代码如下：
#!/usr/bin/perl #Author:zxlee #Function: compute the length of fastq or fasta, fa. #usage: `perl script_name fastq/fasta/fa_file_name`, it will show the total length, also a detail file. use strict; use warnings; my $infile = shift; #give 1st default para to it, you can go on shift to get the 2st para open IN, "<$infile" or die $!; open OUT, ">./result_len.txt" or die $!; our $total_len = 0; our $seq_num = 0; our $len; if($infile =~ /fastq$/){ while(<IN>){ next if not /^@\S+/; my $seq = <IN>; #you cannot use $_ here!!! chomp($seq); $seq_num += 1; $total_len += length($seq); print OUT "\nreads_len = $total_len\n" if $seq_num == 1; } print OUT "Total num of reads is $seq_num\n"; } elsif($infile =~ /(fasta|fa)$/){ # easy way, no "OR" my $chr_len = 0; while(<IN>){ chomp; my $line = $_; if ($line =~ /^>(\S+)/){ print OUT "$chr_len\n" if $chr_len != 0; print OUT "$1\t"; $chr_len = 0; }else{ $len = length($line) if $total_len == 0; $chr_len += length($line); $total_len += length($line); } } print OUT "$chr_len\n"; print OUT "one line has $len\n"; } print "The total length is $total_len\n"; close(IN); close(OUT);
总结：若直接读到数组里，会大大降低操作的复杂度，思维也极容易理解，但是非常耗费计算资源，看数据量决定用哪一种方法。
此脚本主要值得学习的地方：
1.shift的用法，用于逐个取得命令行的参数，让脚本得到封装，不必每次使用时都在里面改文件路径。
2.my、our、local的用法，一般只用my，my管本层及以内的所有层次，my写在最外面就等价于our。
3.$_什么时候会被改变，不是异想天开，while(<IN>)等价于while($_=<IN>)，但是单独<IN>就没有这个效果，这个就是几种局限的用法，不要脑洞大开，随意篡改语法。
3.或的模式匹配，/(fasta|fa)$/)，$就不能写在里面，不知道为什么，记住就好了。
4.算法流程的问题，写程序之前一定要明晰大致的算法，写的时候要是没有算法框架的支持，那你就痛苦死了。
2.1局部组装，创建目录，将比对好的reads按100k为单位，用samtools切，并用awk工具提起reads，分别存放在对应文件夹内
#!/usr/bin/perl #Author: zxlee #Function: ref, origin reads, bwa_sam_filterd_sored.bam, I want split reads per 100k, for local assembly. #Usage: {perl script_name --size=100000 --bam=your_sorted.bam --jobs=200 --outpath=/split_reads_name --sample=DHX_Y255} use strict; use warnings; use Getopt::long; my ($size, $bam, $outpath, $sample_name, $jobs); # get the perl script para GetOptions( "size=s" => \$size, "bam=s" => \$bam, "jobs=s" => \$jobs, "sample=s" => \$sample_name, "outpath=s" => \$outpath, ) or die $1; # my hash, chr => chr_size my %hash = ( 'Chr1' => 43270923, 'Chr2' => 35937250, 'Chr3' => 36413819, 'Chr4' => 35502694, 'Chr5' => 29958434, 'Chr6' => 31248787, 'Chr7' => 29697621, 'Chr8' => 28443022, 'Chr9' => 23012720, 'Chr10' => 23207287, 'Chr11' => 29021106, 'Chr12' => 27531856, ); # iterate each chr, mkdir, create pbs, submit pbs foreach my $chr (keys %hash){ my $split_bam_dir = $outpath."/${chr}_split_bam/"; # dir storage the split bam, abs dir my $split_pbs_dir = $outpath."/${chr}_split_pbs/"; my $split_log_dir = $outpath."/${chr}_split_log/"; my $split_read_names = $outpath."/${chr}_split_read_names/"; system("mkdir $split_bam_dir") if not -e $split_bam_dir; #create dir system("mkdir $split_pbs_dir") if not -e $split_pbs_dir; system("mkdir $split_log_dir") if not -e $split_log_dir; system("mkdir $split_read_names" if not -e $split_read_names; # create split psb for(my $start=1; $start<$hash{$chr}; $start+=$size){ my $end = $start + $size; my $startHK = int($start/$size); # for name convinient my $endHK = int($end/$size); my $bam_name = "$split_bam_dir/${sample_name}_${chr}_${startHK}hk_${endHK}hk.bam"; my $read_name = "$split_read_names_dir/${sample_name}_${chr}_${startHK}hk_${endHK}hk_reads_name.txt"; my $pbs_name = "$split_pbs_dir/${sample_name}_${chr}_${startHK}hk_${endHK}hk.pbs"; open PBS, ">$pbs_name"; print PBS <<SET; #PBS -N split_${startHK}_${endHK} #PBS -j oe #PBS -l nodes=1:ppn=1 #PBS -q low #PBS -o $split_log_dir hostname date cd \$PBS_O_WORKDIR export LD_LIBRARY_PATH=\$LD_LIBRARY_PATH:/public/software/htslib-1.3/lib /public/software/samtools-1.3/bin/samtools view -b $bam $chr:$start-$end -o $bam_name /public/software/samtools-1.3/bin/samtools view $bam_name | awk '{arr[\$1]} END{for(key in arr) {print key} }' > $read_name date SET close PBS; while(1){ my @run_jobs = `qstat -u zxli | grep "split`; if(@run_jobs < $jobs){ lase; }else{ sleep(1); } } system("qsub $pbs_name"); last; # for test } last; # for test }
