搜索
查看: 15305|回复: 31

详细讲解如何抽取fasta文件里面指定序列名的序列

[复制链接]

634

主题

1182

帖子

4030

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4030
发表于 2016-9-8 11:50:00 | 显示全部楼层 |阅读模式
偶尔写一个小白文吧!
要从fasta序列里面抽取指定序列,那么必须要首先了解什么是fasta格式
很简单,比如有下面这个fasta文件
>ID_1
abcd
abcd
abcd
>ID_2
efgh
efgh
efgh
>ID_3
ijkl
ijkl
ijkl
~~~~~~~~~~
如果我想提取序列名为 ID_2的序列,那么我肯定需要写代码。用什么编程语言无所谓啦,我这里就用perl给大家讲解一下。我们先来回顾一下生信菜鸟团群里朋友的问题:



perl问题:大家好,我有个fasta文件,我想从里边把chrM匹配并把匹配的行打印出来 :试写脚本如下                                                                             
[Perl] 纯文本查看 复制代码
#!usr/bin/perl -w
open B,"fasta" or die "$!";
while (my $line = <B>) {
        if($line=~/chrM/i)
                {print $line}
} 

没结果输出是何问题呢????







他的代码错误在打印出来的就是chrM这个序列名而已,而不是序列。要做到提取序列名+序列,首先需要用编程的思想来解析一下整个任务!
首先,要提取文件里面的序列,必须要用代码读取文件。
然后,我们要根据某个特定的序列名来提取,所以我们碰到了>开头的序列名,就需要判断一下,是不是我们想要的。
最后,我们对每一行都需要判断,是否需要打印出来,所以肯定有一个控制变量来决定是否要打印!
我最讨厌的就是写正规代码和注释,但是为了讲解清楚这个最基础代码,我破例一次,就在上面问问题的同学代码基础上面修改一下:
[Perl] 纯文本查看 复制代码
#!usr/bin/perl -w
open B,"fasta" or die "$!"; ##  要保证你的fasta文件的文件名就是fasta
my $sign=0; ## 这个是控制是否打印的变量
while (my $line = <B>) {  ## 一行行的读取
         chomp $line;  ## 去掉末尾的换行符
        if($line =~ /^>/){ # 判断是否是fasta序列名
                if ($line eq '>ID_2' ){  ## 如果是你需要的序列名
                        $sign=1;  ## 就把控制变量设为 1 , 意思是可以输出啦
                }else{
                        $sign=0;  ## 如果不是你需要的序列名,就不输出
                }
        }
        print "$line\n" if $sign==1; ## 根据控制变量的情况来决定是否打印输出
}


如果是我自己写代码,我一般是这样写
[Perl] 纯文本查看 复制代码
$sign=0;
while(<>){
        chomp;
        if(/^>/){
                $_ eq '>ID_2' ? ($sign=1):($sign=0);
        }
        print "$_\n" if $sign==1;
}
其实,这么简单的需求,我一般是不写代码的, 就用perl的单行命令即可!都能很轻松的得到我们想要的序列:
>ID_2
efgh
efgh
efgh














上一篇:生信刷题打卡贴
下一篇:annovar对人类基因组variants注释
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复

使用道具 举报

634

主题

1182

帖子

4030

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4030
 楼主| 发表于 2016-10-31 00:46:19 | 显示全部楼层
这其实是一个很经典的问题的衍生!
就是如何根据一系列ID文件从序列文件里面提取序列!
比如我有ID文件如下:
head longest.transcript.ID
TR3|c0_g1_i1
TR6|c0_g1_i1
TR6|c0_g1_i2
TR6|c0_g1_i3
TR6|c0_g1_i4
TR6|c0_g1_i4
TR6|c0_g1_i5
TR19|c0_g1_i2
TR19|c0_g1_i3
TR19|c0_g1_i3
~~~~~~~~~~~~~~~~~~~~此处省略一万字
然后我有序列文件Trinity.fasta,需要提取最长转录本的序列,所以就需要写程序。
perl choose.pl longest.transcript.ID Trinity.fasta >longest.transcript.fasta
程序如下:
[Perl] 纯文本查看 复制代码
open FH, $ARGV[0];
while(<FH>){
	chomp;
	$_=">".$_ if $_ !~ /^>/;
	$h{$_}=1;
}
close FH;
open FH, $ARGV[1];
$sign=0;
while(<FH>){
        chomp;
        if(/^>/){
			@F=split;
            exists $h{$F[0]} ? ($sign=1):($sign=0);
        }
        print "$_\n" if $sign==1;
}
close FH;


亲测,可以用,请不要再质疑了!
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复 支持 2 反对 0

使用道具 举报

4

主题

24

帖子

184

积分

注册会员

Rank: 2

积分
184
发表于 2017-2-22 00:33:30 | 显示全部楼层
其实最好的方案是用samtools。
首先是对fasta文件建立索引: samtools faidx test.fa.  
产生的索引文件里有每个sequence的名字和长度以及硬盘保存地址的offset。

然后取序列: samtools faidx test.fa seqname 就可以了。
回复 支持 1 反对 0

使用道具 举报

8

主题

52

帖子

324

积分

版主

Rank: 7Rank: 7Rank: 7

积分
324
发表于 2016-9-8 12:32:12 | 显示全部楼层
[Shell] 纯文本查看 复制代码
perl -e 'open ia,"hg19.fasta";$con = 0;while(<ia>){chomp;if(/^>/){($_ eq ">chr1")?($con = 1)$con = 0);}print "$_\n" if ($con eq 1);}'

嘎嘎,我以前是使用temp来做的,就是名字的继承,那样可以做做其他的·~
我的微博:dulunar
回复 支持 1 反对 0

使用道具 举报

0

主题

3

帖子

47

积分

新手上路

Rank: 1

积分
47
发表于 2016-9-8 16:04:24 | 显示全部楼层
楼主这个,还是只是打印了ID2这个名字吧。
回复 支持 1 反对 0

使用道具 举报

634

主题

1182

帖子

4030

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4030
 楼主| 发表于 2016-9-8 17:27:34 | 显示全部楼层
Huffy 发表于 2016-9-8 16:04
楼主这个,还是只是打印了ID2这个名字吧。

你测试了?
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复 支持 反对

使用道具 举报

0

主题

3

帖子

47

积分

新手上路

Rank: 1

积分
47
发表于 2016-9-9 15:39:13 | 显示全部楼层

不好意思,弄错了,是可以打印的。这个方法比我之前做的要简单些,我还设置了>作为while的分行符。赞!
回复 支持 反对

使用道具 举报

0

主题

8

帖子

140

积分

注册会员

Rank: 2

积分
140
发表于 2016-9-12 20:15:48 | 显示全部楼层
我用的是python,也写过类似的一个,原理是一样的。从命令行里读取要提取的序列,如Chr1,然后就会输出到屏幕,linux下可以重定向到文件里。

[Python] 纯文本查看 复制代码
import sys
name = sys.argv[2]
flag = 0
with open(sys.argv[1]) as f:
    for line in f:
        line = line.strip()
        if line[0] == '>':
            if line.split()[0][1:] == name:
                flag = 1
                print(line)
            else :
                flag = 0
        else:
            if flag == 0:
                pass
            else:
                print(line)


回复 支持 反对

使用道具 举报

0

主题

8

帖子

140

积分

注册会员

Rank: 2

积分
140
发表于 2016-9-12 20:20:03 | 显示全部楼层
但是这种思路会把文件全部遍历一遍,最好是能输出完后就停止遍历,效率可以高一点,因为python里面一行行读着好像还挺慢的,但是具体做法一时没想出来,也懒得去搞了,快也快不了一分钟。
回复 支持 反对

使用道具 举报

1

主题

15

帖子

126

积分

注册会员

Rank: 2

积分
126
发表于 2016-9-12 21:46:49 | 显示全部楼层
本帖最后由 starhqy2 于 2016-9-12 21:48 编辑

用awk 更简单点, 把xxx替换成名字就可以了
[Shell] 纯文本查看 复制代码
awk -v RS='>' 'match($0, "xxx"){print ">"$0}' file.fa




回复 支持 反对

使用道具 举报

634

主题

1182

帖子

4030

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4030
 楼主| 发表于 2016-10-31 00:44:35 | 显示全部楼层
starhqy2 发表于 2016-9-12 21:46
用awk 更简单点, 把xxx替换成名字就可以了
[mw_shl_code=shell,true]awk -v RS='>' 'match($0, "xxx"){pri ...

你回答的是对的,awk也是比较好用的命令。但是需求没有那么简单,这其实是一个很经典的问题的衍生!
就是如何根据一系列ID文件从序列文件里面提取序列,请看下面的更贴
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树    

GMT+8, 2019-2-21 22:52 , Processed in 0.082427 second(s), 33 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.