在命令行获取标准输入序列的反互序列,pep序列和长度信息

最近对序列文件处理的比较多,时常要看一些核酸序列的反向互补序列,长度,可能的翻译序列。以前我常常使用seqBuider 来查看,如果能在命令行直接查看,想必是极好的。

这是一个perl脚本,不过我把它的执行路径写入环境变量后,就可以当linux命令直接使用了,很方便的。

这个脚本有四个参数。【-i -r -p -l 】

其中

-i 是必要的参数,用来接收标准输入;

-r 是获得一段序列的反向互补序列(50个字符一行的格式输出);

-p 是提供一段序列的ORF框架序列,即三种可能的pep翻译(50个字符一行的格式输出);

-l 获取一段序列的长度。

如果【-r-p-l】都是缺省状态的话,默认三种结果都输出。

在linux配置文件 ~/.bashrc 文件可以写入:

alias tfa=' perl /yourpath/transfa.pl'

这样以后在linux命令行执行 tfa 命令,出现:

Usage: tfa <STDIN>[-i-r-p-l]

这样的使用提示。

整个代码如下:

#! /usr/bin/perl -w
use strict;
use Getopt::Long;
my ($i,$r,$p,$l);
GetOptions(
	"i!"=>\$i,
	"r!"=>\$r,
	"p!"=>\$p,
	"l!"=>\$l,
);
my $usage = "\nUsage: tfa <STDIN>[-i-r-p-l]\n";
die "$usage\n" unless $i;
print "Please input the nucleotide sequence,and end by ctrl+D.\n\n";
unless($r || $p || $l){
	($r,$p,$l)=(1,1,1);
}
my $fa;
do{local $/;chomp($fa=<STDIN>)};
$fa =~ s/\s+//g;
die "$usage\n" unless $fa;

if($r){
	my $faout = reverse_complement($fa);
	$faout = out_fasta($faout,50);
	print "\n###rc###\n$faout\n";
}
if($p){
	my @fa_arr = cds2pep($fa);
	print "\n###protein###\n";
	$fa_arr[0] = out_fasta($fa_arr[0],50);
	print "ORF1:\n$fa_arr[0]\n";
	$fa_arr[1] = out_fasta($fa_arr[1],50);
        print "ORF2:\n$fa_arr[1]\n";
	$fa_arr[2] = out_fasta($fa_arr[2],50);
        print "ORF3:\n$fa_arr[2]\n";

}
if($l){
	my $len = length $fa;
	print "\n###Length###\n$len\n";
}
#####################
sub out_fasta{
        my ($seq,$num) = @_;
        my $len = length $seq;
        $seq =~ s/([A-Za-z]{$num})/$1\n/g;
        chop($seq) unless $len % $num;
        return $seq;
}
#####################
sub reverse_complement{
        my ($seq)=shift;
        $seq=reverse$seq;
        $seq=~tr/AaGgCcTt/TtCcGgAa/;
        return $seq;
}
#####################
sub cds2pep{
        my $seq=shift;
	##phase0
	my $str0 = $seq;
	$str0 = trans($str0);
	##phase1
	my $str1 = substr($seq,1);
	$str1 = trans($str1);
	##phase0
        my $str2 = substr($seq,2);
        $str2 = trans($str2);
	return ($str0,$str1,$str2);
}
#####################
sub trans{
	my $seq = shift;
	my $p = code();
	my $out;
	for(my $i=0;$i<length$seq;$i+=3){
                my $codon=uc(substr($seq,$i,3));
                last if (length$codon <3);
                $out.= exists $p->{"standard"}{$codon} ? $p->{"standard"}{$codon} : "X";
        }
        return $out;
}
#####################
sub code{
        my $p={
                "standard" =>
                        {
                                'GCA' => 'A', 'GCC' => 'A', 'GCG' => 'A', 'GCT' => 'A',                               # Alanine
                                'TGC' => 'C', 'TGT' => 'C',                                                           # Cysteine
                                'GAC' => 'D', 'GAT' => 'D',                                                           # Aspartic Aci
                                'GAA' => 'E', 'GAG' => 'E',                                                           # Glutamic Aci
                                'TTC' => 'F', 'TTT' => 'F',                                                           # Phenylalanin
                                'GGA' => 'G', 'GGC' => 'G', 'GGG' => 'G', 'GGT' => 'G',                               # Glycine
                                'CAC' => 'H', 'CAT' => 'H',                                                           # Histidine
                                'ATA' => 'I', 'ATC' => 'I', 'ATT' => 'I',                                             # Isoleucine
                                'AAA' => 'K', 'AAG' => 'K',                                                           # Lysine
                                'CTA' => 'L', 'CTC' => 'L', 'CTG' => 'L', 'CTT' => 'L', 'TTA' => 'L', 'TTG' => 'L',   # Leucine
                                'ATG' => 'M',                                                                         # Methionine
                                'AAC' => 'N', 'AAT' => 'N',                                                           # Asparagine
                                'CCA' => 'P', 'CCC' => 'P', 'CCG' => 'P', 'CCT' => 'P',                               # Proline
                                'CAA' => 'Q', 'CAG' => 'Q',                                                           # Glutamine
                                'CGA' => 'R', 'CGC' => 'R', 'CGG' => 'R', 'CGT' => 'R', 'AGA' => 'R', 'AGG' => 'R',   # Arginine
                                'TCA' => 'S', 'TCC' => 'S', 'TCG' => 'S', 'TCT' => 'S', 'AGC' => 'S', 'AGT' => 'S',   # Serine
                                'ACA' => 'T', 'ACC' => 'T', 'ACG' => 'T', 'ACT' => 'T',                               # Threonine
                                'GTA' => 'V', 'GTC' => 'V', 'GTG' => 'V', 'GTT' => 'V',                               # Valine
                                'TGG' => 'W',                                                                         # Tryptophan
                                'TAC' => 'Y', 'TAT' => 'Y',                                                           # Tyrosine
                                'TAA' => 'U', 'TAG' => 'U', 'TGA' => 'U'                                              # Stop
                         }
                ## more translate table could be added here in future
                ## more translate table could be added here in future
                ## more translate table could be added here in future
        };
        return $p;
}

__END__

时间: 2024-12-23 12:47:18

在命令行获取标准输入序列的反互序列,pep序列和长度信息的相关文章

Linux学习之命令行获取公网IP方法详解

本文和大家分享的主要是主要是linux命令行获取公网IP相关内容,一起来看看吧,希望对大家学习linux有所帮助. curl ipinfo.io $ curl ipinfo.io { "ip": "36.10.25.4", "city": "Hangzhou", "region": "Zhejiang", "country": "CN", &quo

用Linux命令行获取本机外网IP地址

用Linux命令行获取本机外网IP地址 $ curl ifconfig.me$ curl icanhazip.com$ curl ident.me$ curl ipecho.net/plain$ curl whatismyip.akamai.com$ curl tnx.nl/ip$ curl myip.dnsomatic.com$ curl ip.appspot.com$ curl -s checkip.dyndns.org | sed 's/.*IP Address: \([0-9\.]*\)

直接从命令行获取MySQL查询语句结果

如果你需要直接从命令行获取MySQL查询语句结果,那么你可以使用-B和-N这两个参数来达到目的. 例:获取MySQL用户数. [[email protected] ~]# mysql -BN -uroot -predhat mysql -e 'select count(*) from user' 6 [[email protected] ~]# -B参数:去掉边框 -N参数:只显示结果

命令行获取mobileprovision文件的UUID

有时候我们需要获取Provisioning Profile的UUID,一般情况下我们是双击mobileprovision文件安装,在老版本Xcode的Organizer窗口里找到对应的UUID,新版本Xcode里不知道为什么这个路径下就找不到uuid了,幸运地我找到了另外一个软件:iPhone配置实用工具,这个软件里也能查看. 只是每次都要借助于GUI工具来查看还是太麻烦了,更好的方法是用命令行,我找到了这个:0xc010d/mobileprovision-read 使用方法: 在Termina

iOS命令行获取工程内所有的国际化资源并且整合

之前做国际化的时间,每次写一个国际化的字符串都要到.strings里面对应添加资源,这样很容易遗忘掉一些国际换的资源. 而且写起来很麻烦.一个个写太费劲了. 今天心血来潮,正好考虑到国际化的事情,就上网查看了下资料,果然有利器啊. 之前写国际化的方法 太low了. 下面一一道来. test: 1,首先,新建一个工程,随意写. 2,随意写几个国际化的资源,类似: NSLocalizedString(@"你好1", nil); 3,关键来了 打开控制台 进入你当前工程的目录下 进入目录后,

nodejs 命令行获取入参

安装:npm install yargs --save-dev Example index.js const argv = yargs.alias('n', 'name').alias('p', 'port').alias('t', 'tpl').argv console.log(argv.name); 命令行输入: $ node index.js -n fuck -p 8080 -tpl components 原文地址:https://www.cnblogs.com/CyLee/p/84344

go语言 从命令行获取参数解析

go语言内置的flag包实现了命令行参数的解析,flag包使得开发命令行工具更为简单. os.Args 如果你只是简单的想要获取命令行参数,可以像下面的示例代码一样使用os.Args来获取命令行参数 package main import "fmt" import "os" func main(){ if len(os.Args)>0 { for index , value := range os.Args { fmt.Println(index, value

命令行获取苹果电脑的主要硬件配置

最近公司想收一台苹果电脑,Boss在他的朋友圈里喊话,收到了一大堆响应.我发觉很多买了Mac的人用不来,电脑都闲置着,甚至连自己的电脑的基本配置都不知道,所以这个事情也挺困难…… 如何最方便地获取到电脑的主要硬件配置?要是在Windows下,肯定是到控制面板的设备管理器里去看看,或者安装一个第三方的软件,国产的貌似就很多,但我现在对装软件很有抵触,原因是这些免费的国产软件里捆绑了太多太多让你一不小心就装上的流氓软件,设备管理器嘛,信息很多,可有用的却不多,像显卡这种东西又不太看得出来型号,点来点

win下命令行获取管理员权限

在win下运行npm install安装依赖出现错误: Error: EBUSY, resource busy or locked 搜索错误信息后发现是由于没有管理员权限,在bash中输入以下命令后运行,再次安装依赖就完全成功啦.   runas /noprofile /user:Administrator cmd 最后吐槽一下,windows就是破事多...