logo
搜索
我的图书馆

发文章

发文工具

其他工具

留言交流
搜索
分享

该学点编程知识了

 生物_医药_科研 2019-01-11

   大神一句话，菜鸟跑半年。我不是大神，但我可以缩短你走弯路的半年~

   就像歌儿唱的那样，如果你不知道该往哪儿走，就留在这学点生信好不好~

   这里有豆豆和花花的学习历程，从新手到进阶，生信路上有你有我！

豆豆写于19.1.8-1.10  这部分内容学自《Perl在生物信息学中应用》
这次是继上次整理后又学习的全部内容，也是得益于之前零零散散学了一些编程知识，学起来不像刚开始一样费劲，但是也有许多模糊的概念需要巩固

我从不纠结学习哪种语言，因为对于新手，学什么都是困难的。
一切从解决问题出发。关于代码，先看懂，才能创造

常用元素

  • 命令解释

    第一行#起始：#! /usr/bin/perl -w ，告诉计算机这个是perl程序。-w表示遇到错误发出警告

  • 语句

    分号结尾，如赋值语句$DNA='ATCGCTACG'; ，将右边的真实DNA信息存储在perl程序中的DNA变量中

  • 变量

    只由大小写字母、数字、下划线组成，不能数字开头

    $开头的叫标量变量，储存字符串或各种数字，但是一个标量一次只能存储一种类型【不能以数字开头；变量名大小写敏感】

  • 换行

    print $DNA \n\n表示输出DNA后中间空一行再打印其他

  • 单引号与双引号

    单引号属于强引用，其中的任何东西都会输出为字符，例如：'\n'会打印输出换行符，'\n'只会打印\n

  • 文件句柄

    就是处理文件需要用的东西，例如$protein=中的PROTEIN

    大小写没有规定，但一般使用大写表示；

    其中尖括号表示输入操作符，就是利用这个从外部读入数据，一般操作是：

    'NM_021964fragment.pep';
    open( PROTEINFILE, proteinfilenameEE′NMEEfragmentpep′openPROTEINFILE
    # 然后就可以把读入的PROTEINFILE存到数组中

    常用的句柄还有：STDIN（从键盘上读取用户的输入）；STDOUT 把输出内容输出到屏幕（默认）或到其他位置；错误信息定向到STDERR ，例如：print STDERR 'Bad info\n';

  • pop 数组末尾拿掉一个元素；shift开头拿掉一个元素；
    unshift添加一个元素到数组开头；push添加一个元素到数组结尾

  • Perl中有三种主要数据类型：标量变量、数组和散列/哈希（关联数组）

    数组

    @开头，是储存多个标量的变量

    @bases = ('A', 'C', 'G', 'T');
    print @bases; #输出结果没空格
    print '@bases'; #输出结果有空格

    $dna[2] 表示对数组@dna取第3个元素（因为起始索引为0）

    数组可以用=拷贝：@out = @in

    $count = @input 得到的标量变量count是数组input中元素数目

    image.png

    散列/哈希

    散列以%起始，如my %dictionary = (...);。取元素时，散列和数组一样用$，因为查找返回的都是标量值，散列与数组不同之处就是：散列用{}，数组用[]

    获取散列的所有键构成的数组：@keys = keys %my_hash; ；获取散列所有值构成的数组：@values = values %my_hash; 利用键值对可以快速提取信息，例如基因名作为键，基因相关信息作为值【注意：不能把正则表达式或字符集作为散列的键，必须是简单的标量值，比如一个字符串或一个数字】

    对散列进行排序：@sorted_keys = sort keys %my_hash;

    image.png

  • 括号：在复杂的运算中指明运算的顺序，如：(2+3)/3 ，不加括号就先运行3/3

有用的操作

连接

  • 使用.，$DNA3=$DNA1.$DNA2 ，然后print '$DNA3\n';

  • $DNA='$DNA1DNA2'，然后print '$DNA3\n';

  • print $DNA1, $DNA2, '\n';或print '$DNA1$DNA2\n';

  • $protein=join('',@protein)，需要指定一个字符串，会放在数组各个元素之间，当分隔字符串为空字符串时，就将各个元素连起来

替换--s///

image.png

注意到=~，它是绑定操作符，就是将右边的操作应用到左边变量的字符串上；

然后s/T/U/g，s表示替换（substitution），第一个斜线后跟着T，表示字符串中被替换的元素；第二个斜线后跟着U，是要替换成的元素；g表示全局

因此，这句命令的意思就是：把$RNA变量存储的字符串数据中所有T变成U

另外，如果要删除换行符或者空格或制表符，使用$RNA =~ s/\s//g

反转数组--reverse

reverse @bases   ，例如：

@bases = ('A', 'C', 'G', 'T');
@reverse = reverse @bases;
print '@reverse\n\n';
# T G C A

互补--tr

tr 其实也是替换的作用，但是对于求序列的互补比较方便，第一个斜线后是要被替换的内容；第二个斜线后是替换后的内容，且位置对应进行替换

image.png

因为tr函数会返回在字符串中找到特定字符的数目，而且即使替换字符集为空，原始字符集也不会改变，因此还可以作为计数器，如：$base = ($dna =~ tr/ACGTacgt//);就可以统计所有碱基数

获取长度--scalar

scalar @bases ，对于上面的例子，得到的结果就是4

任意位置插入--splice

splice ，例如：

@bases = ('A', 'C', 'G', 'T');
splice ( @bases, 2, 0, 'X' );
print '@bases\n';
# A C X G T

参考：http://blog.sina.com.cn/s/blog_4d3a41f40100m29m.html

拆解--split

@DNA=split('', $DNA)将$DNA拆解成单个字符的@DNA数组【与join相反】

取子集--substr

$base = substr($DNA, $position,1)

变量初始化

声明变量后需要赋一个值，如果不进行初始化，它的值将被定义为'undef'，对于undef，在数字操作中表示0，在字符串操作中表示空字符串

循环

  • foreach循环【对一个数组中的元素进行循环处理时比较方便】

  # 统计碱基数目
  foreach $base (@DNA){
      $count++; #效果等同于$count += 1
  }

当然，如果不指定从数组中读取的标量变量，如foreach (@DNA)，perl会使用特殊变量$_

  • while循环

  while ($DNA =~ /a/ig) {$a++} #统计A碱基的数目(i表示不区分大小写)

  • for循环示例

  # DNA翻译成蛋白序列(codon2aa是一个子程序，记录了密码子与氨基酸的对应)
  for (my $i=0; $i<>length($dna)-2);$i += 3){
      $protein .= condon2aa (substr($dna, $i, 3));
  }
  return $protein;

  • 其他一些常用格式

  while(CONDITION) {BLOCK}
  until(CONDITION) {BLOCK}
  for(INITIALIZATION ; CONDITION ; RE-INITIALIZATION ) {BLOCK}
  foreach VAR (LIST) {BLOCK}
  for VAR (LIST) {BLOCK}
  do {BLOCK} while (CONDITION)
  do {BLOCK} until (CONDITION)

排序

  • 按字母顺序升序：@array = sort @array;

  • 按数字升序：@array = sort {$a <=> $b } @array;

范围操作符

(1../comments/) and next; ，其中(..) 是范围操作符，表示跳过第一行到包含comments的所有行。然后and是一个逻辑操作符，表示跳过后继续进行

逻辑操作符

and可以用&&表示，or 可以用|| 表示，not用!表示

用逻辑操作符来控制流程：比如打开文件

open (FH, $filename) or die 'Cannot open file $filename: $!';

快速调用匹配值

# 例如：'123456' =~ /34/ 可以匹配成功
$&就表示匹配上的内容'34'
$`表示匹配成功前面的内容'12'
$'表示匹配成功后面的内容'56'

正则表达

要想在正则表达式中体现//，可以用/\/\/\n/ ，但总感觉有点奇怪，像这种情况，可以在前面使用m，利用其它的界定符（一般我们使用/，但并不代表只有这一种界定方式），从而避免在正斜线前使用转义符，如m!//\n!  

  • 看一个例子：while( $features =~ /^ {5}\S.*\n(^ {21}\S.*\n)*/gm ) {};

    这是while循环中添加的一个正则表达式匹配：如果features满足一定的要求，就执行while语句下面的内容，那么满足什么要求呢？那就需要看表达式

    不要被看似复杂的结构搞晕，直接看/ /中间包围的那部分：其中\n表示换行符，那我们先看第一个\n前面的内容：
    ^ {5}\S.*\n 表示以5个空格开头^ {5}，然后跟着一个非空白字符\S，再接任意个非换行符字符.*

    然后第二个\n前面的是(^ {21}\S.*\n)* 表示 以21个空格开头，然后加一个非空白字符，再接任意个非换行字符。并且注意到这部分内容被小括号包围，然后跟了星号，表示小括号中的内容可有可无。

    因此，/ /内部整体上说明的是： 匹配第一行以5个空格开头，下面的行中以21个空格开头，但这部分可有可无

  • 点号加量词操作符.{num}表示字符个数

  • (.*?)表示非贪婪匹配或最小化匹配，会尽可能匹配短的字符串（默认*是贪婪的）

  • 关于元字符：\ | ( ) [{ ^ $ * + ? .
    匹配它们需要\转义，使用|进行单一匹配，使用()用于分组
    使用[]设置匹配字符集（如[A-Za-z0-9]表示任意单个字母或数字， [^0-9] 匹配非数字）
    使用.匹配任意一个除换行符外的字符
    使用? 限定量词进行最小化匹配，用在* + {} 后面
    使用* + {MIN,} {MIN,MAX} 表示量词，其中*为0次、1次或多次，+为1次或多次，{}指定几次,可选最少和最多

  • Perl根据小括号从左到右出现的顺序来将匹配到的字符串存储到$1、$2等特殊变量中，例如：

  $alphabet = 'abcdefghijklmnopqrstuvwxyz';
  $alphabet =~ /(((a)b)c)/;
  print 'First pattern = ', $1,'\n';
  print 'Second pattern = ', $2,'\n';
  print 'Third pattern = ', $3,'\n';

  #First pattern = abc
  #Second pattern = ab
  #Third pattern = a

另外小括号将变量括起来还可以捕获变量，例如：

  #如果不加括号，变量只会存储数组中的元素个数
  @array = ('one', 'two');
  $a = @array;
  print $a;
  # 结果为2

  #加上括号，就变成了只有一个元素的列表，结果就是数组在列表中按行求值，并且括号中有几个元素，就能捕获到数组中几个变量
  @array = ('one', 'two');
  ($a,$b) = @array;
  print '$a $b';
  #结果为 one two

模式修饰

常用的模式修饰符：用于全局匹配的/g，不分区大小写的/i，将整个字符串当成单独一行的/s（忽略换行符），将整个字符串作为含有换行符的行/m （/m可以让元字符^和$ 匹配嵌在内部换行符之前和之后的位置）

搭配元字符： ^锚定开头，$锚定结尾，.匹配除了换行符以外的任意字符

例如：要对AAC\nGTT 进行匹配

use warnings;
'AAC\nGTT' =~ /^.*$/; #使用^.*$可以匹配除了换行符外的任意字符
print $&, '\n';
# 结果报错

# 原因就在于：在中间出现了一个换行符，默认的模式修饰不能识别，因此可以转换成m修饰，会从字符串开头一直匹配到第一个换行符(即便出现在内部)
'AAC\nGTT' =~ /^.*$/m;
print $&, '\n';
#结果输出AAC，$&就是匹配到的内容

# 如果用/s，就可以匹配包括换行符在内的所有
'AAC\nGTT' =~ /^.*$/s;
print $&, '\n';
# AAC
# GTT

输出

  • printf：示例：

  my $first = '3.14159265';
  my $second = 76;
  my $third = 'Hello world!';

  printf STDOUT 'A float: %6.4f An integer: %-5d and a string: %s\n', 
  $first, $second, $third;
  #结果生成
  A float: 3.1416 An integer: 76      and a string: Hello world!

基本的指令就是prinft % 标识 转换指令符 ，转换指令符就是说明输出的内容格式，包括：用于浮点数的f、整数d、字符串s。另外在%与转换指令符之间还有标识，对应的是最小字段宽度、精度、长度。

例如：%6.4f 表示输出一个浮点数，最小宽度是6个字符，小数部分是四位，于是可以看到，在结果的3.1416前空了1位，就是为了满足6个字符长度的要求；

%-5d指定输出整数，字段宽度为5，- 表示数字居左

子程序

子程序可以重复使用，可以单独进行测试，子程序本身可以调用其他子程序，总而言之，子程序就是只做一件事并把它做好，然后我们要做的就是：组合它们来完成整个任务

编写子程序示例

子程序一般包括：sub、子程序名字、大括号代码块

子程序从@_ 数组中收集参数值，然后赋给子程序的变量（子程序的变量名可以和主程序重名，并不冲突）。如果只在子程序中写my ($dna)而不加= @_，那么子程序就会使用主程序的同名变量或者报错找不到新变量

sub addbase {
    my ($dna) = @_; #perl的特殊数组是“_”。my用来声明一个新变量，只在包裹的代码块中有效
    $dna .= 'ATGC';
    return $dna
}

切记：子程序中定义任何变量，都要使用my ，可以在程序开始添加use strict; 提醒我们使用my定义变量

例如：计算DNA中G的数目

#!/usr/bin/perl -w
use strict; #强制声明变量

# 主程序
# 定义变量USAGE作为程序帮助文档，其中包含所有命令行参数[这里包括了$0和DNA]
# $0也是特殊变量，表示这个perl脚本名；DNA指的是需要的DNA序列
my ($USAGE) = 'USAGE：$0 DNA\n\n'
# 数组变量@ARGV是内置变量，其中包含输入的参数；如果没有输入参数就打印程序帮助文档USAGE
unless (@ARGV){
    print $USAGE;
    exit;
}
my ($dna) =$ARGV[0]; #提取数组第一个元素，注意要换成$提取(perl的坐标从0开始)
my ($sum_Gs) = countG($dna);
print '\nThe DNA $dna has $sum_Gs G\'s in it!\n\n';
exit;

#子程序
sub countG{
    my ($dna) = @_;
    my ($count) = 0;
    $count = ($dna =~ tr/Gg//);
    return $count;
}

传递数据给子程序

  • 值传递：（参数比较简单）主程序参数复制后传给子程序，子程序中这部分值改变并不会影响主程序相应参数的值，例如：

  #!/usr/bin/perl -w
  use strict;

  my $i=2;
  subvalue($i);
  print 'In the main program, the result is $i\n';
  exit;

  sub subvalue{
      my ($i)=@_;
      $i += 100;
      print 'In the sub program, the result is $i\n';
  }

  #In the sub program, the result is 102
  #In the main program, the result is 2

  • 引用传递：（参数比较复杂）可以向子程序中传递标量、数组等各种形式的参数，子程序中所做的任何事都改变主程序中的对应参数。重点是：在引用的变量名前加一个反斜线

  #!/usr/bin/perl -w
  use strict;

  #主程序
  my @i = ('1','2','3');
  my @j = ('a','b','c');
  print 'Mian-1 program: i = '.'@i\n';
  print 'Mian-1 program: j = '.'@j\n';

  subvalue(\@i,\@j);
  print 'Main-2 program: i = '.'@i\n';
  print 'Main-2 program: j = '.'@j\n';
  exit;

  #子程序
  sub subvalue{
      my($i,$j)=@_;
      print 'Sub program: i = '.'@i\n';
      print 'Sub program: j = '.'@j\n';
      push(@$i,'4');
      shift(@$j);
  }

  #Mian-1 program: i = 1 2 3
  #Mian-1 program: j = a b c
  #Sub program: i = 1 2 3
  #Sub program: j = a b c
  #Main-2 program: i = 1 2 3 4
  #Main-2 program: j = b c

关于指针：只有一个$的变量都是标量变量的指针，作用仅仅是当一个中间桥梁的作用；如果要直接表示指向的那个标量变量，需要用两个$，同理，如果要表示数组变量就用@$，表示散列变量%$

重复使用的子程序—模块/库

子程序确实很方便，但是一些比较好的子程序需要应用到许多脚本中，一旦需要修改就需要修改所有脚本中的子程序，这是非常麻烦的。因此，如果一开始将子程序放到一个方便使用的文件中，如果改动就只改一次，然后可以应用到所有的脚本就好了。因此，模块（库）就诞生了。

例如，我们有一个库叫beginperl.pm（后缀名是perl module的意思），里面存放了所有的子程序。【需要注意的是，.pm文件的最后一行要是1;】。

使用库的时候，直接使用use beginperl; 【注意：这条语句放在靠近use strict;的地方；库的名字不带后缀；库如果和脚本不在一个目录下，需要绝对路径】

例如：File::Find模块帮助定位文件和目录，找到一个好的模块可以减少编程时间

应用

一：提取Genbank序列和注释信息

先对Genbank的信息分布有个大体了解，比如序列信息在最后，从ORIGIN开始到最后一行的//，注释信息是ORIGIN之前的信息

基本思路：加载文件=》遍历数组（文件每一行）=》先匹配// 得到最后一行行号=》再匹配ORIGIN得到分界线行号【之所以先匹配//是因为如果先匹配分界线的话，一旦匹配上就不会再继续进行了，也就得不到最后一行的行号了】=》序列就是ORIGIN行向下，注释就是ORIGIN行向上

#前提得到Genbank文件@Genbankfile
#sub open_file
#    given the filename, return the filehandle
sub open_file {
    my ($filename) = @_;
    my $fh;
    unless (open($fh, $filename)){
        print 'cannot open\n';
        exit;
    }
    return $fh;
}
# 这里只是列一个基本思路
my $linenumber = 0;
my $end;
my $origin;
foreach my $line (@Genbankfile){
    if ($line =~ /^//\n/) {
        $end = $linenumber;
        last;
    }elsif ($line =~ /^ORIGIN/){
        $origin = $linenumber;
    }
    $linenumber++;
}
@annotation = @Genbankfile[0..($origin-1)];
@seq=@Genbankfile[($origin+1)..($end-1)];

关于读取文件：可以在命令行给出要打开的文件名（例如perl a.txt test.pl ），并在perl脚本中第一行使用<>就可以表示读取文件，例如while (<>) {操作}

二、Blast

非编程相关内容

不论是网页还是本地运行的blast，都是基于字符串匹配程序，利用字符串匹配来查找相似，来作为同源的一个指标。相似度（similarity）可以用percent identity来衡量，就是说查询的序列和数据库中的某个序列对应区域完全匹配的碱基数目，它会找到等价密码子或者不改变蛋白功能并具有相似属性的氨基酸残基之间的匹配，表示一个保守的程度。说到序列同源（homology） ，它是指序列在进化上相关，要么同源，要么不同源，没有同源度这种说法

Blast的运行主要基于：原始值S、打分算法、数据库属性。原始值（Raw score S）是对相似性和匹配大小的度量。另外在Blast的结果中还能看到：按照E值排序的hits数。一个匹配的E值/期望值（E value/Expect value）表示一个随机生成的具有同样大小和数据库构成的字符串匹配几率（其中允许空位存在）。E越接近0，表示越不可能随机发生。因此，E越小，匹配结果越好 。一般E小于1是一个比较可靠的结果

资料：

http://barc.wi./education/bioinfo2003/perl_bioinfo.html

https://github.com/Yixf-Education/course_Perl/tree/master/book

想说点啥？

    本站是提供个人知识管理的网络存储空间，所有内容均由用户发布，不代表本站观点。请注意甄别内容中的联系方式、诱导购买等信息，谨防诈骗。如发现有害或侵权内容，请点击一键举报。
    转藏 分享 献花（0）

    来自： 生物_医药_科研 > 《待分类》

    举报/认领

    0条评论

    发表

    请遵守用户 评论公约

    类似文章 更多