生物信息 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
";

模块名是Getopt::Long，不能写错；使用时，得用GetOptions函数，函数中间是个哈希结构，将--first传给$first变量，以此类推。

这个模块有点复杂，具体参照：在Perl中使用Getopt::Long模块来接收用户命令行参数

2.读表格式文件，切分，取出目的字段，建立数据结构

open IN, "test.txt" or die $!;
while(<IN>){
    chomp;
    @list = split / /;
    next if $list[0] =~ /char1|char4|char6/;
    $name = $list[2];
    $value = $list[1];
    $hash{$name} = $value;
    foreach (keys %hash){
        print "$_ => $hash{$_}
";
    }
}
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 "
reads_len = $total_len
" if $seq_num == 1;
    }
    print OUT "Total num of reads is $seq_num
";
}
elsif($infile =~ /(fasta|fa)$/){ # easy way, no "OR"
    my $chr_len = 0;
    while(<IN>){
        chomp;
        my $line = $_;
        if ($line =~ /^>(S+)/){
            print OUT "$chr_len
" if $chr_len != 0;
            print OUT "$1	";
            $chr_len = 0;
        }else{
            $len = length($line) if $total_len == 0;
            $chr_len += length($line);
            $total_len += length($line);
        }
    }
    print OUT "$chr_len
";
    print OUT "one line has $len
";
}
print "The total length is $total_len
";
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
}