搜索
楼主: Jimmy

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

[复制链接]

634

主题

1182

帖子

4030

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4030
 楼主| 发表于 2016-10-31 22:11:10 | 显示全部楼层
这其实是一个很经典的问题的衍生!
就是如何根据一系列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 进行打赏,谢谢
回复 支持 反对

使用道具 举报

1

主题

17

帖子

268

积分

中级会员

Rank: 3Rank: 3

积分
268
发表于 2016-12-29 22:26:45 | 显示全部楼层
最近在联系时发现一个问题。如果我用perl单行处理,同样是需要提取多条序列。我想使用shell的while循环,该怎么进行套用呢?
[Perl] 纯文本查看 复制代码
le longest.transcript.ID|while read aa ;do perl -e 'open dd,"Trinity.fasta";$sign = 0;while(<dd>){chomp;if(/^>/){($_=~/$aa/)?($sign = 1):($sign = 0);}print "$_\n" if ($sign eq 1);}'

也就是该如何去调用读入的没一个参数$aa? 望解答!

回复 支持 反对

使用道具 举报

634

主题

1182

帖子

4030

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4030
 楼主| 发表于 2016-12-30 08:53:39 | 显示全部楼层
btrainee 发表于 2016-12-29 22:26
最近在联系时发现一个问题。如果我用perl单行处理,同样是需要提取多条序列。我想使用shell的while循环,该 ...

好蛋疼,我都懒得看!
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复 支持 反对

使用道具 举报

1

主题

17

帖子

268

积分

中级会员

Rank: 3Rank: 3

积分
268
发表于 2016-12-30 23:11:49 | 显示全部楼层
Jimmy 发表于 2016-12-30 08:53
好蛋疼,我都懒得看!

回复 支持 反对

使用道具 举报

1

主题

41

帖子

285

积分

中级会员

Rank: 3Rank: 3

积分
285
发表于 2017-1-17 15:05:38 | 显示全部楼层
本帖最后由 learnyoung 于 2017-1-17 15:26 编辑

对于fasta和gennebank文件的解析可以用biopython这个包,我选取ls_orchid 前十个gene的id存入gene_id 文件中做测试
[Python] 纯文本查看 复制代码
#coding=utf-8
from Bio import SeqIO
#需要提取的序列的基因ID在gene_id文件中,到ls_orchid 中提取对应的序列
with open('gene_id.txt')as f:
    for gen_id in f:
        gen_id=gen_id.rstrip()
        seq = SeqIO.parse('ls_orchid.fasta.txt', 'fasta')
        for line in seq:
            if line.id==gen_id:
                print gen_id+':'+'\n'+line.seq

输出结果为:

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
回复 支持 反对

使用道具 举报

29

主题

131

帖子

1208

积分

金牌会员

Rank: 6Rank: 6

积分
1208
发表于 2017-1-19 10:19:12 | 显示全部楼层
前段时间刚好学习做无参转录组,对于从Trinity.fasta结果中提取最长转录本也写了个代码,思路可能不是最优,代码如下:
[Perl] 纯文本查看 复制代码
#!/usr/bin/perl -w
use strict;

my $in_file;
my $out_file = "unigene.fasta";

$in_file = "trinity_out_dir.Trinity.fasta";

open my $fh_in, $in_file or die "Cannot open $in_file";
my $header;
my $seq;
my %hash;
my $hash_rev;

while(<$fh_in>){
	chomp;
	if ($_ =~ /^>(.*)/){
		$header = $1;
		$seq = "";
	}else{
		$seq .= $_;
		$hash{$header} = $seq;
	}
}
close $fh_in;

my %hash_rev = reverse %hash;
my %comp_longest_trans;
my %longest_seq;

open my $fh_out, ">$out_file";
foreach my $comp (keys %hash){
	if ($comp =~ /^(TRINITY.*?g\d+).*?len=(\d+)/){
		my $trans = $1;
		my $len = $2;
		if ((! exists $comp_longest_trans{$trans}) 
			 || $comp_longest_trans{$trans} < $len){
			$comp_longest_trans{$trans} = $len;
			$longest_seq{$trans} = $hash{$comp};
		}
	}
}

foreach (keys %longest_seq){
	my $seq_trans = $longest_seq{$_};
	print $fh_out ">$hash_rev{$seq_trans}\n$seq_trans\n";
}
回复 支持 反对

使用道具 举报

11

主题

54

帖子

242

积分

中级会员

Rank: 3Rank: 3

积分
242
QQ
发表于 2017-2-7 13:21:18 | 显示全部楼层
好贴!!!!!
苛求远离完美
回复

使用道具 举报

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 就可以了。
回复 支持 2 反对 0

使用道具 举报

634

主题

1182

帖子

4030

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4030
 楼主| 发表于 2017-2-22 08:27:54 | 显示全部楼层
tsznxx 发表于 2017-2-22 00:33
其实最好的方案是用samtools。
首先是对fasta文件建立索引: samtools faidx test.fa.  
产生的索引文件里有 ...

最好的方案?
很抱歉,我不这么认为,一万个生信工程师有一万个需求,难道我们都要把samtools软件的细枝末节全部背诵下来才能做数据分析吗?
你这个解决方案的确很好,略掉了自己编程的过程,但初学者最好不要这样,误入歧途。
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复 支持 反对

使用道具 举报

5

主题

37

帖子

485

积分

中级会员

Rank: 3Rank: 3

积分
485
发表于 2017-2-22 09:23:01 | 显示全部楼层
好贴,学习了两个新方法,biopython和samtools
回复 支持 反对

使用道具 举报

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

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-7-21 04:57 , Processed in 0.045011 second(s), 24 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.