Perl的没落

目录: 生物信息 | 标签: Perl | 发表时间: Jan 7, 2020

Perl一直是生信领域工作人员习惯使用的编程语言，在处理序列文本方面有着得天独厚的优势。但是，目前生信领域可使用 Perl 的场景也都可以使用 Python，另外 Python 还在是数据分析机器学习相关方法的重要编程语言，我认为 Perl 会逐渐退出生信分析的历史舞台。由于Perl 5 和 Perl 6 不相互兼容，Perl 6在去年正式改名叫Raku，不少人也会从 Perl 转向其他社区活跃的编程语言。早在2010年左右，国内的 Perl China 还在组织一些定期的交流活动，但是现在已经找不到这个网站。使用Perl的这批人，已经进入了中年，没有时间和精力活跃在社区。Perl China 邮件列表里上一封邮件的发送日期还是2019年5月，目前估计都转战其他联系平台（微信群和slack）。

那么为啥现在还有那么多生信初学者在学习Perl呢？个人观点是由于目前做到PI的人，在他们学习生信的时期，主要用的是 Perl 语言，所以由于这种惯性，我们还会在很多 lab 招聘信息上看到要求应聘的生信博士们要会Perl。常用的由Perl编写的最有名生信软件是Annovar。除此之外，bcftools和samtools里一些小工具是Perl编写的，如果你用不到的话，可能根本就没注意过这些Perl代码。

翻了翻我自己写的 Perl 代码，基本上就做了以下10件事情，以前写东西特别 naive，最后一个代码超级长且都是不必要的重复。最近，我写 Perl 代码的时候越来越少了，除了修改代码 bug 或者写一行命令搞定个小测试外，已经用不到 Perl。

1.计算序列中某些信息的出现频率（GC content）

#!usr/bin/perl -w
$name="Homo_sapiens.GRCh37.66.dna.chromosome";
$input="/data/GENE_transcripLoc.bed";
open FILE,"$input"|| die "Couldn't open $input: $!";

while (defined ($eachline =<FILE>)) {
  chomp $eachline;
  ($chr,$begin,$end,$strand,$geneID,$transID,$other)=();
  ($chr,$begin,$end,$strand,$geneID,$transID,$other)=split /\t/, $eachline;
  $data="/data/data4/yu/homo_sapiens_dna/$name.$chr.fa";
  open DATA,"$data" or die "Couldn't open $data: $!";
  chomp(@DNA = <DATA>);
  $line=join("",@DNA);
        $str=substr($line,$begin,$end); 
        my $num_a    = $str =~ tr/A//;
        my $num_c    = $str =~ tr/C//;
        my $num_g    = $str =~ tr/G//;
        my $num_t    = $str =~ tr/T//;
        #my $num_o    = length($str) - $num_a - $num_c - $num_g - $num_t;
        my $gc       =( $num_c + $num_g ) / length($str)  * 100 ;
        $outpromoter="/data/data4/yu/GCcontent_genebody.bed";
        open PLOG,">>$outpromoter" or die "Couldn't open $outputpromoter: $!";
        print PLOG "$chr\t$begin\t$end\t$strand\t$geneID\t$transID\t$gc\n";
        close PLOG;
  close DATA;
}
close FILE;

2.剪切序列（trim）

#!/usr/bin/perl -w

while (<DATA>) {
    s/^(?!>)...(.*)..../$1/;
    #       ^^^    ^^^^  
    #       x=3    y=4
    print;
}

__DATA__
>HWI-EAS158_40_3_1_46_535
GTGAATGCGTGATACAGGAATGTTCGTTGTGACCAT
>HWI-EAS158_40_3_1_47_579
AAAGTGAATGCGTGATACAGGAATGTTCGTTGTGAC
>HWI-EAS158_40_3_1_46_731
GTGTCATGCGTGATACAGGAATGTTCGTTGTGAAAA
GTGTCATGCGTGATACAGGAATGTTCGTTGTGAAAA

3. 密码子转氨基酸 (抄自《Beginning Perl for Bioinformatics》)

#!/usr/bin/perl
use strict;
use warnings;
#Initialize variables
my $dna='CGACGTCTCGTACGG';
my $protein=''
my $codon;
#Translate each three-base codon into an amino acid,and append to a protein
for(my$i=0;$i<(length($dna)-2);$i+=3){
  $codon=substr($dna,$i,3);
    $protein.=codon2aa($codon);
}
printf"I translated the DNA\n\n$dna\n\n into the protein\n\n$protein\n\n"
exit;
#codon2aa
#A subroutine to translate a DNA 3- character codon to n amino acid
sub codon2aa{
  my($codon)=@_;
    if ($codon=~/GC./i){return'A'} #Alanine
      elsif ($codon=~/TG[TC]/i){return'C'} #Cysteine
        elsif ($codon=~/GA[TC]/i){return'D'} #Aspartic Acid
        elsif ($codon=~/GA[AG]/i){return'E'} #Glutamic Acid
        elsif ($codon=~/TT[Tc]/i){return'F'} #Phenylalanine
        elsif ($codon=~/GC./i){return'G'} #Glycine
        elsif ($codon=~/CA[TC]/i){return'H'} #Histidine
        elsif ($codon=~/AT[TCA]/i){return'I'} #Isoleucine
        elsif ($codon=~/TT[AG]|CT./i){return'L'} #Leucine
        elsif ($codon=~/ATG/i){return'M'} #Methionine
        elsif ($codon=~/AA[TC]/i){return'N'} #Asparagine
        elsif ($codon=~/CC./i){return'P'} #Proline
        elsif ($codon=~/CA[AG]/i){return'Q'} #Glutamine
        elsif ($codon=~/CG.|AG[AG]/i){return'R'} #Arginine
        elsif ($codon=~/TC.|AG[TC]/i){return'S'} #$erine
        elsif ($codon=~/AC./i){return'T'} #Threonine
        elsif ($codon=~/GT./i){return'V'} #Valine
        elsif ($codon=~/TGG/i){return'W'} #Tryptophan
        elsif ($codon=~/TA[TC]/i){return'Y'} #Tyrosine
        elsif ($codon=~/TA[AG]|TGA/i){return'_'} #Stop
    else{
      print STDERR"Bad codon\"$codon"!!\n";
        exit;
    }
}

4.下载序列信息（wget）

超级不建议用这个程序，要安装LWP模块，其实wget加个循环就可以下载了。

#!/usr/bin/perl
#use strict;
#use warnings;

my @down;
while($line=<DATA>){
  chomp $line;
  my $url= "http://www.ncbi.nlm.nih.gov/nuccore/$line?report=gilist&log\$=seqview&format=text";
  use LWP::Simple;
  my $content = get $url;
  print "Got $line Gi ID!\n";
  if($content=~/<pre>(.*)\n<\/pre>/s){
    push @download, $1;
  }
} 

getstore("http://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nucleotide&rettype=fasta&retmode=text&id=".join(",",@download),"seqs.fasta");
print "Download all sequences!\n";
__DATA__
GU132860
GU132861
GU132862
GU132863

5.计算kmer组合

这个代码可以学习push和shift这样的操作。

#!/usr/bin/perl

$k = shift;
@bases = ('A','C','G','T');
@words = @bases;
for $i (1..$k-1){
  undef @newwords;
  foreach $w (@words){
    foreach $b (@bases){
      push (@newwords,$w.$b);
    }
  }
  undef @words;
  @words = @newwords;
}
foreach $w (@words){
  print "$w\n";
}

6.分割FASTA文件

这段代码可以学习$/，$1这样的标志是什么意思。

#!/usr/bin/perl   
usestrict;   
my($file, $dir) = @ARGV; #取得输入文件和输出文件夹   
open FILE, $file  
    or die "cannotopen$file:$!";    #打开文件   
local $/ = ">"; #设置读取结束标记为>,默认为\n   
my $i = 1;  #为输出文件编号   
while(<FILE>){   
    if(/^>$/){   
        next;   
    }   
    s/>//;  #去掉序列结束位置的>号   
    s/(\w+)/>$1/;   #为分割出来的序列第一行添加>号   
    my $name = $1;  #以fasta序列第一行的一个单词作为文件名   
    open OUTPUT, ">", $dir . '/' . $name . '_' . $i . '.fa';   
    print OUTPUT$_;   
    close OUTPUT;   
    $i++;   
}   
close FILE;  

7.将一段长序列变成FASTA格式

#!usr/bin/perl -w
$name="FRT.40A.chromosome";
$data="/share/test";
open DATA,"$data" or die "Couldn't open $data: $!";
chomp(@DNA = <DATA>);
$line=join("",@DNA);
my $length=length($line);
my $re=$length%50;
my $num=($length-$re)/50;
 
for ($i=1;$i<=$num;$i=$i+1){
    $str=substr($line,($i-1)*50+1,50);
    $output="/share/disk5-1/zhangzhhgroup/yulij/result/$name.fa";
    open LOG,">>$output" or die "Couldn't open $output: $!";
    print LOG "$str\n";
    close LOG;
  }

close DATA;

8.查找子序列位置

同KMP算法功能类似。

#usr/bin/perl -w
# initialize the original string
my $input = "Perl index, Perl chop, Perl hex";
my $pos = 0;

my @positions = ();

while(1){
  $pos = index(lc $input, "perl", $pos);
  last if($pos < 0);
  push(@positions, $pos++);
}       

if(scalar (@positions) > 0){
  print "substring found at positions: ", "@positions\n";
}else{
  print "substring not found\n";
}

9.找到与已知信息匹配的行

其实就是个grep命令。

#!usr/bin/perl -w
#$gen_file="/root/testperl.txt";
#open DATA, $gen_file or die "cannot open file";
$id="590168";
while (<DATA>) {
    print $_ if $_ =~ /$id/;
}
#close DATA;

__DATA__
COG1110 604354.TSIB_0530        1       1212
COG1110 604354.TSIB_0530        1       1212
COG1110 604354.TSIB_0530        1       1212
COG1110 598659.NAMH_1071        1       1060 
COG1110 593117.TGAM_1679        1       1216 
COG1110 590168.Tnap_0803        1       1104 
COG1110 579137.Metvu_1264       1       865 
COG1110 573064.Mefer_0049       6       869 
COG1110 563041.HPG27_964        428     531
COG1110 590168.Tnap_0803        1       1104

10. 用Perl提取网页信息

这是用perl写过的最长的一个程序，抓取和讯网站上的董秘信息，可能现在改改代码还能用。如果让我现在写这个程序，肯定直接用Python美丽汤（Beautiful Soup）。

#!usr/bin/perl -w
#use WWW::Mechanize;
#use HTTP::Cookies;

use   LWP;
use   LWP::Simple; 
use LWP::UserAgent;
use Encode qw/encode decode from_to/;


use encoding "utf-8";
binmode(STDIN, ':utf8');
binmode(STDOUT, ':utf8');
binmode(STDERR, ':utf8');


print "股票代码|公司简称|公司全称|法定代表人|董秘|联系电话|传真|注册地址|办公地址|邮编|邮箱\n";
my @k="300001".."300999";
my $k;
foreach $k(@k){
my $url = "http://stockdata.stock.hexun.com/gszl/s$k.shtml";#"http://stockdata.stock.hexun.com/gszl/s000001.shtml";#'http://vip.stock.finance.sina.com.cn/corp/go.php/vCI_CorpInfo/stockid/000001.phtml';

my $content = get $url;
die "Couldn't get $url" unless defined $content;

my $begin="\<\!--平安银行公司简介--\>";
my $end="\<\!--平安银行公司简介--\>";


if($content=~ m/$begin(.*)$end/sgm){
  my $a="公司简称";
  my $b="股票代码";
  my $yyy=$1;
  
  if ($yyy=~m#$a(.*?)$b#sgm){
    my $zzz=$1;
    if ($zzz=~m/<td>(.*?)<\/td>/sgm){
        print "$1|";
      }else {
  print "匹配失败|";
}
  }else {
  print "匹配失败|";
}

}else {
  print "匹配失败|";
}
 
$content = get $url;
die "Couldn't get $url" unless defined $content;

$begin="\<\!--平安银行公司简介--\>";
$end="\<\!--平安银行公司简介--\>";

if($content=~ m/$begin(.*)$end/sg){
  my $a="股票代码";
  my $b="公司全称";
  my $yyy=$1;
  if ($yyy=~m#$a(.*?)$b#sg){
    $zzz=$1;
    if ($zzz=~m/<td>(.*?)<\/td>/sg){
        print "$1|";
      }else {
  print "匹配失败|";
}
  }else {
  print "匹配失败|";
}
}else {
  print "匹配失败|";
}


$content = get $url;
die "Couldn't get $url" unless defined $content;

$begin="\<\!--平安银行公司简介--\>";
$end="\<\!--平安银行公司简介--\>";

if($content=~ m/$begin(.*)$end/sg){
  my $a="公司全称";
  my $b="公司英文名称";
  my $yyy=$1;
  if ($yyy=~m#$a(.*?)$b#sg){
    $zzz=$1;
    if ($zzz=~m/<td>(.*?)<\/td>/sg){
        print "$1|";
      }else {
  print "匹配失败|";
}
  }else {
  print "匹配失败|";
}
}else {
  print "匹配失败|";
}



$content = get $url;
die "Couldn't get $url" unless defined $content;

$begin="\<\!--平安银行公司简介--\>";
$end="\<\!--平安银行公司简介--\>";

if($content=~ m/$begin(.*)$end/sg){
  my $a="法定代表人";
  my $b="独立董事";
  my $yyy=$1;
  if ($yyy=~m#$a(.*?)$b#sg){
    $zzz=$1;
    if ($zzz=~m/<td>(.*?)<\/td>/sg){
        print "$1|";
      }else {
  print "匹配失败|";
}
  }else {
  print "匹配失败|";
}
}else {
  print "匹配失败|";
}




$content = get $url;
die "Couldn't get $url" unless defined $content;

$begin="\<\!--平安银行公司简介--\>";
$end="\<\!--平安银行公司简介--\>";

if($content=~ m/$begin(.*)$end/sg){
  my $a="联系人";
  my $b="邮政编码";
  my $yyy=$1;
  if ($yyy=~m#$a(.*?)$b#sg){
    $zzz=$1;
    if ($zzz=~m/<td>(.*?)<\/td>/sg){
        print "$1|";
      }else {
  print "匹配失败|";
}
  }else {
  print "匹配失败|";
}
}else {
  print "匹配失败|";
}


$content = get $url;
die "Couldn't get $url" unless defined $content;

$begin="\<\!--平安银行公司简介--\>";
$end="\<\!--平安银行公司简介--\>";

if($content=~ m/$begin(.*)$end/sg){
  my $a="董秘";
  my $b="公司传真";
  my $yyy=$1;
  if ($yyy=~m#$a(.*?)$b#sg){
    $zzz=$1;
    if ($zzz=~m/<td>(.*?)<\/td>/sg){
        print "$1|";
      }else {
  print "匹配失败|";
}
  }else {
  print "匹配失败|";
}
}else {
  print "匹配失败|";
}



$content = get $url;
die "Couldn't get $url" unless defined $content;

$begin="\<\!--平安银行公司简介--\>";
$end="\<\!--平安银行公司简介--\>";

if($content=~ m/$begin(.*)$end/sg){
  my $a="公司传真";
  my $b="电子邮箱";
  my $yyy=$1;
  if ($yyy=~m#$a(.*?)$b#sg){
    $zzz=$1;
    if ($zzz=~m/<td>(.*?)<\/td>/sg){
        print "$1|";
      }else {
  print "匹配失败|";
}
  }else {
  print "匹配失败|";
}
}else {
  print "匹配失败|";
}




$content = get $url;
die "Couldn't get $url" unless defined $content;

$begin="\<\!--平安银行公司简介--\>";
$end="\<\!--平安银行公司简介--\>";

if($content=~ m/$begin(.*)$end/sg){
  my $a="注册地址";
  my $b="所得税率";
  my $yyy=$1;
  if ($yyy=~m#$a(.*?)$b#sg){
    $zzz=$1;
    if ($zzz=~m/<td>(.*?)<\/td>/sg){
        print "$1|";
      }else {
  print "匹配失败|";
}
  }else {
  print "匹配失败|";
}
}else {
  print "匹配失败|";
}



$content = get $url;
die "Couldn't get $url" unless defined $content;

$begin="\<\!--平安银行公司简介--\>";
$end="\<\!--平安银行公司简介--\>";

if($content=~ m/$begin(.*)$end/sg){
  my $a="办公地址";
  my $b="主要产品";
  my $yyy=$1;
  if ($yyy=~m#$a(.*?)$b#sg){
    $zzz=$1;
    if ($zzz=~m/<td>(.*?)<\/td>/sg){
        print "$1|";
      }else {
  print "匹配失败|";
}
  }else {
  print "匹配失败|";
}
}else {
  print "匹配失败|";
}


$content = get $url;
die "Couldn't get $url" unless defined $content;

$begin="\<\!--平安银行公司简介--\>";
$end="\<\!--平安银行公司简介--\>";

if($content=~ m/$begin(.*)$end/sg){
  my $a="邮政编码";
  my $b="公司简介";
  my $yyy=$1;
  if ($yyy=~m#$a(.*?)$b#sg){
    $zzz=$1;
    if ($zzz=~m/<td>(.*?)<\/td>/sg){
        print "$1|";
      }else {
  print "匹配失败|";
}
  }else {
  print "匹配失败|";
}
}else {
  print "匹配失败|";
}



$content = get $url;
die "Couldn't get $url" unless defined $content;

$begin="\<\!--平安银行公司简介--\>";
$end="\<\!--平安银行公司简介--\>";

if($content=~ m/$begin(.*)$end/sg){
  my $a="mailto:";
  my $b="公司网址";
  my $yyy=$1;
  if ($yyy=~m#$a(.*?)$b#sg){
    $zzz=$1;
    if ($zzz=~m/target=\"_blank\">(.*?)<\/a>/sg){
        print "$1\n";
      }else {
  print "匹配失败\n";
}
  }else {
  print "匹配失败\n";
}
}else {
  print "匹配失败\n";
}


}

常用下载方式总结

目录: 软件世界 | 标签: aria2, wget, rsync, bypy, scp | 发表时间: Jan 6, 2020

rsync断点续传

网络传输速度很大程度上会影响我的工作。服务器与服务器之间的通连，有时scp不是最好的选择，例如：网络出现问题，两个服务器之间的连接断开后，使用scp再次连接时无法进行断点续传。

使用rsync做服务器之间的文件传输（备份），可以做到断点续传。

rsync -avzP -e 'ssh -p 8888' abc@xxx.xxx.xxx.xxxx:/home/abc/filename ./

wget并行下载

-P 8八线程，-t 0无限重试，-nv 下载时只显示更新和出错信息，不显示指令的详细执行过程

就是输入的地址urls.txt，每一行用一个wget命令进行下载，一次最多同时下载8个文件。

xargs -P 8 -n 1 wget -nv -c -t 0 <urls.txt

wget下载动态连接的文件

wget -c 'http://some.site.com/download?id=234&status=download' -O output_filename

随时上网检索一下，上面这个方法不一定每个动态连接的文件都能下载，例如google网盘的。

aria2分段下载大文件

-s 3 分成3份，-x 3 用3个连接从服务器下载文件

aria2c  -s 3 -x 3 "https://baidu.com/d/datadump.current.zip"

bypy将百度云盘数据下载到服务器

安装信息：(bypy)[https://github.com/houtianze/bypy]

测试发现，大文件上传比慢。文件只能上传到全部文件>apps>bypy这个文件夹。

下载一系列FTP上的文件

这是压箱底的老程序，filetype这里其实不需要使用一个循环，如果只有一个类型的话，直接在下载的循环中用grep筛出这个变量。

#!/bin/bash

for filetype in pdf; 
do 
  for net in `grep "$filetype" file.list2 | sed "s/\<td\>\<a href=\"//g" | sed "s/\"//g" | sed "s/<TD//g"`; 
  do 
    if [ ! -f $net ]
    then
    echo $net
    #wget  $path/$net
    fi
  done
done

Fibonacci数列

目录: 计算方法 | 标签: Fibonacci, 递归, 时间复杂度, 通项公式 | 发表时间: Jan 3, 2020

接上篇内容，翻倒远古时期用MATLAB写的Fibonacci数列第N项的计算方法，直觉上写一个递归函数就可以搞定问题，但是其计算时间复杂度是O(2^n)。

当然，如果需要降低时间复杂度，应该写成循环的形式。我一开始以为自己写的就是递归函数的版本，查着资料，写着博客，猛然发现自己写的就是循环。

%古老的脚本，从第三项开始计算数列的内容
function y=Fibonacci(k)
    a=[1 1];
        for i=3:k
        a(i)=a(i-1)+a(i-2);
        end
    y=a(k);
end

如果还想偷懒，可以写成矩阵计算的样子。

\[\begin{bmatrix}f_{n+1} \\ f_{n}\end{bmatrix}=\begin{bmatrix}1 & 1 \\ 1 & 0\end{bmatrix}\cdot\begin{bmatrix}f_{n} \\ f_{n-1}\end{bmatrix}\]

function y=Fibonacci(n)
  a=[1 0];
  b=[1 1;1 0];
  y=a*b^(n)*[0;1];
end

到这里，其实还没进入我想说的内容。我为了查询Fibonacci数列的算法，搜索了一些资源，基本总结有4种方法（递归，循环，矩阵，通项），最后一种就是用通项公式进行计算。

看到第四种方法时，我懵了，Fibonacci数列还有通项公式？等我再看到通项公式的内容，才发现，这个公式在高中时期就已经学过了。对于公式的内容，自己一直没有记住过。温故知新，这次从头学了一遍通项公式的计算。

具体方法，即从上面的矩阵方法中，求$\begin{bmatrix}1 & 1 \\ 1 & 0\end{bmatrix}$的特征值和特征向量，然后带回到方程中，已知$f_{0}=\begin{bmatrix}1 \\ 0\end{bmatrix}$ 然后可求通项¹。

网上还有一些初级数学的方法可以求解，但是我忘记高中时候是怎么推导的了（很可能是用差分方法），如果有机会，翻翻之前的笔记，尽量找出来。

Fibonacci number: Matrix form ↩

写代码

目录: 聚类不能 | 标签: 代码 | 发表时间: Jan 2, 2020

最近在整理硬盘时发现了一个久远的代码文件夹，里面有我从大学开始写的数学计算程序和研究生时期的一些处理数据的脚本。回顾了一下，发现数学相关的程序很多都记不清算法了，研究生时期的脚本写得又特别的青涩（low）╮(￣▽￣)╭。

为了腾出几M的硬盘空间。我决定把一些常用的方法都挪到博客上，温故知新，重新学习一遍（学了也会忘记的，不经常用的方法是不会记在脑内的，就当个松鼠，收集整理好，存在网络上）。

有些完全没有技巧性的，我就挪到英文博客了。中文博客记录一些可以展开说明的内容。

如无特殊说明，所有程序均按照初期使用的语言编写，其中MATLAB替换成Octave，Python是3.6以上版本，R是最新版本，Perl语言用Perl5，c和c++用电脑里GCC最新版，操作系统是CentOS或者Fedora。

在这里，还要说一下写代码的习惯，最开始的时候，仗着自己用Linux系统，我连代码程序文件的后缀.pl，.sh都不写，这样造成的重大问题是，光看文件名，我自己都不清楚这个程序是用什么语言写的。在自己研究生时期的研究工作硬盘彻底报废的时候，抢救回来的数据，很多都是文本格式的乱码，如果当时写了文件名后缀，那么在数据恢复的时候，软件可以根据后缀来选择恢复的数据类型。

其次，写代码要有注释，下面是一个规范的有注释的MATLAB代码。我的其他MATLAB代码，都没写注释，由于MATLAB都是矩阵运算，读起来就很费力了。

function dijkstra(V,o)
%V为邻接矩阵，o为零点标号，输出向量d是零点到每个点的距离
A=V;%定义可变矩阵
[m,n]=size(A);
d=zeros(1,m);%零点到个点的距离
d(:)=inf;
Q=A(o,:);%零点所在行的矩阵
d(o)=0;%零点到零点的距离为0
p=10;%循环标志
while(min(Q)~=inf&p~=0)%判断条件矩阵为inf或距离没有inf项时停止
    [a id]=min(Q);%求最小值及其标号
    Q(id)=inf;%消去已经求出的值
    A(o,id)=inf;%消去已经求出的值
    o=id;%转换行列
    if d(id)>a%所求距离比已知小的时候代换
        d(id)=a;
    end
    for j=1:m%改写矩阵的距离
        if (Q(j)>d(id)+A(id,j))
            Q(j)=d(id)+A(id,j);
        end
    end
    p=0;%循环标志设0值
    for i=1:m%检查是否有inf项
        if d(i)==inf
            p=p+1;
        end
    end
end
d
end

最后，程序里要写metainfo，例如编写时间，修改时间，修改内容，作者，联系方式等等。尤其在开源软件中，写metainfo会方便记录代码是否有多人进行维护（有点类似于git commit）

2020年5月3日再次提醒自己，写代码一定要有注释，否则，你根本不记得自己三个月前写了什么。。

染色体坐标系

目录: 生物信息 | 标签: 坐标系, 染色体 | 发表时间: May 20, 2019

一开始学生物信息分析的时候，知道BED格式中染色体坐标系是从0开始的。用久了各种工具，突然有一天发现SAM和BAM格式的染色体坐标系起始位置是不一样的。猛然间发现生物信息学的分析工具中各种坐标不统一。读了一些英文资料，觉得有必要在博客里记录一下各种坐标系的位置关系。

以0为起点

sequence	A	C	G	T	A
0-based	0	1	2	3	4	5

以1为起点

sequence	A	C	G	T	A
1-based	1	2	3	4	5	6

除了其实为点的坐标位置不同，终止位置的坐标是否包括最后一个位置，在不同的格式和工具里也是不一样的。

例如上面序列中的GTA的坐标位置：

以0为起始位点，包含最后一个位置: 2-4

以1为起始位点，包含最后一个位置: 3-5

以1为起始位点，不包含最后一个位置: 3-6

下面记录一些工具的染色体坐标系定义位置：

以0为起始位点：SAM, VCF, GFF 和 Wiggle 格式
以1为起始位点：BAM, BCFv2 和 PSL 格式
Ensembl：以1为起始位点，包含最后一个位置。https://asia.ensembl.org/info/docs/api/core/core_tutorial.html
UCSC Genome Browser：内部表示，以0为起始位点，并用以1为起始位点时的终止位点（程序员为了编程方便煞费苦心）。外部显示，以1为起始位点。http://genome.ucsc.edu/FAQ/FAQtracks.html#tracks1
BED 格式：以0为起始位点，不包含最后一个位置。BED格式经常使用，在分析中别弄错起始位置和终点位置非常重要 http://genome.ucsc.edu/FAQ/FAQformat.html#format1
MAF 格式：以1为起始位点，包含最后一个位置。

参考文献

Chromosome coordinate systems: 0-based, 1-based

On Your Mark

Never give up, Dream big