搜索
查看: 768|回复: 7

从fasta序列中随机提取部分序列

[复制链接]

20

主题

84

帖子

635

积分

高级会员

Rank: 4

积分
635
发表于 2017-1-9 00:04:56 | 显示全部楼层 |阅读模式
本帖最后由 anlan 于 2017-1-10 00:36 编辑

我有一个大约有300万条序列的fasta文件,我想从中提取10%的序列,我的思路如下:
1. 用希哈将序列ID与序列对应起来;
2. 将序列ID写入一个数组中;
3. 自定义一个函数生成不重复的随机数;
4. 利用随机数来对应数组的下标,提取出10%的序列ID,然后用哈希生成对应序列。

我对上述步骤3用perl写了个代码,总的序列条数为2779988,我想提取的条数为277998,即从0-2779988中生成277998个随机数,如下:
[Perl] 纯文本查看 复制代码
#!/usr/bin/perl -w
use strict;

my $percent_num = 277998;
my $all_num = 2779988;

my @id_list = &random_num($all_num,$percent_num);
open my $fh_out, ">count_id.txt" or die "Cannot write count_id.txt";
map{
        print $fh_out "$_\n";                
}@id_list;

sub random_num {
        my ($k,$count) = @_;
        my @tmp;
        my %repeat;
        while(1){                
                my $num = int(rand($k));
                my $i = 0;
                foreach my $info (@tmp){
                        if ($info == $num){
                                $i = 1;
                                last;
                        }
                }                
                if ($i == 0){
                        push @tmp, $num;
                }
                my $j = $#tmp + 1;
                if ($j == $count){
                        return(@tmp);
                        last;
                }
        }
}

这是我最开始写的代码,总运行时间为45分钟,后来发现在自定义函数中,去除重复的随机数那步用foreach循环太消耗时间了。。

接下来是我改进的代码,用希哈来去除重复的随机数,如下:
[Perl] 纯文本查看 复制代码
#!/usr/bin/perl -w
use strict;

my $percent_num = 277998;
my $all_num = 2779988;

my @id_list = &random_num($all_num,$percent_num);
open my $fh_out, ">count_id.txt" or die "Cannot write count_id.txt";
map{
        print $fh_out "$_\n";                
}@id_list;


sub random_num {
        my ($k,$count) = @_;
        my @array;
        my %repeat;
        while(1){
                my $num;
                my $tmp = 0;
                do{
                        $num = int(rand($k));
                        $tmp = ++$repeat{$num};
                }while($tmp > 1);
                @array = keys %repeat;
                my $j = $#array + 1;
                if ($j == $count){
                        return(@array);
                        last;
                }
        }
}

总运行时间跟上面一个方法一样长。。说明检查是否重复那步消耗太多时间了,因此还是不行。

另外个方法,完全跟前面的思路不一样:
先初始化一个1-n的数组,然后每次循环取一个随机位置,将这个位置的元素和第一个位置的元素交换!

[Perl] 纯文本查看 复制代码
sub random_num {
	my ($k,$count) = @_;
	my $i;
	my @array = 1..$k;
	for ($i=0;$i<$count;$i++){
		my $index = int(rand($k-1));
		my $a = $array[$i];
		my $b = $array[$index];
		$array[$index] = $a;
		$array[$i] = $b;
	}
	my $j = 0;
	my @tmp;
	while($j<$count){
		push @tmp, $array[$j];
	 	$j++;
	}
	return @tmp;
}


这个方式能生成n个不重复随机数,消除了循环嵌套,10S内就能完成。。我的思路是先生成1到n的n个数字的数组,然后从中随机挑选一个替换数组中的第一个元素,然后再从中随机挑选一个替换数组中第二个元素,依次运行n次,然后取前m个作为我所需要的m个随机数。。。不知道对不对。。




上一篇:GO,Kegg,GSEA富集分析神包clusterprofiler
下一篇:Python (excel )一个问题
回复

使用道具 举报

11

主题

44

帖子

213

积分

中级会员

Rank: 3Rank: 3

积分
213
发表于 2017-1-9 14:58:43 | 显示全部楼层
首先指出你代码里的错误
[Perl] 纯文本查看 复制代码
do{
            $num = int(rand($k));
            $tmp = $repeat{$num}++;
        }while($tmp > 1);


这里应该是
$tmp = ++$repeat{$num};


否则第二次判断tmp还是1
这个需要去查自增计算相关知识

然后是你代码里的缺陷
你用hash去记是能接受的,但既然用了hash为什么还要用array来再存一次~直接keys就好了
然后你这里影响效率的地方是当 rand($k) 重复多次出现时你的计算量会增加非常多

最后说一个高效的方法
惯例,使用伪代码
m = 100
n = 10

for i in 0..m
    list = i
for i in 0..n
    index = random(m-1)
    swap(list[index], list)
return list[1..n]

回复 支持 反对

使用道具 举报

20

主题

84

帖子

635

积分

高级会员

Rank: 4

积分
635
 楼主| 发表于 2017-1-10 00:31:34 | 显示全部楼层
dongye 发表于 2017-1-9 14:58
首先指出你代码里的错误[mw_shl_code=perl,true]do{
            $num = int(rand($k));
            $tmp  ...

我确实写错了,已修改。
但没看懂python的swap函数,网上也没查到,是自定义函数吗?
我按照自己的理解也写了个perl的
回复 支持 反对

使用道具 举报

0

主题

12

帖子

101

积分

注册会员

Rank: 2

积分
101
发表于 2017-1-10 14:43:25 | 显示全部楼层
一秒钟以内完成随机数生成,利用random模块
[Python] 纯文本查看 复制代码
import random
numList=list(range(3000000))
randNum=random.sample(numList,300000)
回复 支持 反对

使用道具 举报

0

主题

12

帖子

101

积分

注册会员

Rank: 2

积分
101
发表于 2017-1-10 14:43:57 | 显示全部楼层

................

.....................
回复

使用道具 举报

20

主题

84

帖子

635

积分

高级会员

Rank: 4

积分
635
 楼主| 发表于 2017-1-11 12:17:37 | 显示全部楼层
sxuan 发表于 2017-1-10 14:43
一秒钟以内完成随机数生成,利用random模块
[mw_shl_code=python,true]import random
numList=list(range ...

是不重复的随机数吗
回复 支持 反对

使用道具 举报

0

主题

12

帖子

101

积分

注册会员

Rank: 2

积分
101
发表于 2017-1-11 20:12:50 | 显示全部楼层
anlan 发表于 2017-1-11 12:17
是不重复的随机数吗

恩  你试试看不就知道了   
回复 支持 反对

使用道具 举报

2

主题

3

帖子

26

积分

新手上路

Rank: 1

积分
26
发表于 2017-3-2 13:04:03 | 显示全部楼层
大家去读读李恒的 samtools 源码。他不是这么做的。我这里就直接用了 samtools 源码 https://github.com/huboqiang/NOM ... nzymeSites/main.cpp
回复 支持 反对

使用道具 举报

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

本版积分规则

QQ|关于我们|手机版|小黑屋|生信技能树    

GMT+8, 2017-8-24 05:13 , Processed in 0.031734 second(s), 29 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.