根据序列ID提取fasta序列

根据序列ID快速抽提fasta序列

在进行完序列比对以及，在微生物多样性分析中常常要根据物种信息抽提特定ID的fasta格式的序列文件。传统方法你们确定是一个个ID去查找，复制，粘贴。这种方法虽然很精确可是效率实在是过低了。若是能用几行命令解决不只简单高效，并且还能一劳永逸，妈妈不再用担忧我花几天时间在查找序列上了。linux

下面就正式教你们，怎么把两个文件名字不同的文件中特定想要的序列文件一一匹配并提取出来。小程序

# 用grep命令匹配全部有>号的行，即全部序列名字的行，而后另存为文件 ITS_all.name

grep ">" ITS_all.fas > ITS_all.name

## 用sed查找>，并把全部>号删除
sed 's/>//g' ITS_all.name > ITS_all.txt

# 在windows下的格式在Linux环境下可能不识别须要转换格式

dos2unix ITS_all.name
dos2unix ITS.name

# 经过grep -f能够在前者的文件中查找并匹配后者文件中的全部内容，而后经过sort找出unique的序列ID，另存为ITS.out文件
grep -f ITS.name ITS_all.name | sort -u > ITS.out
## 用perl小程序进行提取
perl ./extract.pl ITS.out ITS_all.fas > ITS_extract.fa

那么这个extract.pl

这个小程序的具体代码是什么呢？windows

你们用的时候把它包装成extract.pl的脚本形式就行其内部代码以下：

#!/opt/perl/bin/perl

use warnings;
use strict;
use Bio::SeqIO;
use File::Basename;
use Getopt::Long;
use Pod::Usage;
use Bio::Index::Fasta;
use File::Find;
use Data::Dumper;

my $seqID=shift or die "sequence ID?";
my $fastafile=shift or die "fasta file?";

my $s_wrote1;
my $s_read1;
my $print_it1;
$s_read1 = $s_wrote1 = $print_it1 = 0; 

open(ID,"<$seqID"); 

#open(FASTA,">$seqID.fasta");
my %ids1=();
while (<ID>) { 
  s/\r?\n//; /^>?(\S+)/; $ids1{$1}++; 
} 
my $num_ids1 = keys %ids1; 

open(F, "<$fastafile");
while (<F>) { 
  if (/^>(\S+)/) { 
    $s_read1++; 
    if ($ids1{$1}) { 
      $s_wrote1++; 
      $print_it1 = 1; 
      $_ = join("\n", $_);
      #$_=join("\n", $new, $_);
      #$_=join($new, $_);
      delete $ids1{$1} 
    } else { 
      $print_it1 = 0  
    }
  }
  if ($print_it1) { 
     print $_;
#    print FASTA $_; 
  }
} 
warn "Searched $s_read1 FASTA records.\nFound $s_wrote1 IDs out of $num_ids1 in the ID list.\n";

grep命令经常使用参数介绍

linux的使用过程当中会常常使用到匹配两个文件的操做markdown

grep文件匹配时用到的参数

grep进行文件内容匹配工做是用到的参数主要有两个，分别是 **1. 取出两个文件中的相同部份内容“-wf”参数. app

1. 取出两个文件中的不一样部份内容“-wvf”参数 **

取出两个文件中相同的行

grep -wf AAA.txt BBB.txt
111
222
ddd
fff
asda
1
2
3
4
5
6
7

取出两个文件中不一样的行

grep -wvf aaa.txt bbb.txt
233
ccc
1
2
3
4

此时取出的是bbb.txt文件中存在而aaa.txt文件中没有的内容spa

grep -wvf bbb.txt aaa.txt
333
1
2
3

此时取出的是aaa.txt文件中存在而bbb.txt文件中没有的内容 unix

grep -wvf firstfile secondfile规则就是取出secondfile中存在可是在firstfile中不存在的内容并输出到屏幕上code