搜索
查看: 3343|回复: 3

把含有简并碱基的引物序列还原成多条序列-高级难度

[复制链接]

634

主题

1182

帖子

4030

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4030
发表于 2017-1-15 21:14:40 | 显示全部楼层 |阅读模式
新手切勿入坑!
引物是 ATGCVHGTAACTAGCT,其中VH是简并碱基,需要还原,如果就这么一条,很容易还原:

如果是批量的,就很头疼了,我写了两次才写对,你们可以试一下,用perl/python/R均可!

R:agY:CTM:ACK:GTS:gcW:ATH:atcB:gtcV:gacD:GATN:ATgc把这个文本拷贝到txt文件里面保存好,或者直接写入hash里面也行!
我处理过程,代码十分之复杂:
[Perl] 纯文本查看 复制代码
while(<DATA>){
chomp;
@F=split/:/;
$hash{$F[0]}=uc $F[1];
} ##这里记录简并碱基的对应关系
## %hash stored the tables;
sub primer2multiple{
$primer=$_[0];
$prod=1;
$primer_len=length $primer ;
foreach $i (0..$primer_len-1){
$char=substr($primer,$i,1);
#$prod*=length $hash{$char} if ($char !~/[ATCG]/) ;
if ($char !~/[ATCG]/) {
push @pos_list,$i;
push @char_list,$hash{$char};
##首先找出所有的不是ATCG的碱基位置以及它对应的碱基
## record all of the positions which are not ATCG;
}
}
@out_list=($primer);
##循环处理每个不是ATCG的碱基位置,让它们根据对应关系扩展
foreach my $i (0..scalar(@pos_list)-1){
@out_list=&new_out_list(\@out_list,$pos_list[$i],$char_list[$i]);
} ##&new_out_list 这个函数非常重要,会把数组不停的扩展,最终达到应该有的个数!
print join"\n",@out_list;
print "\n";
}
sub new_out_list{
my @array = @{$_[0]};
my $pos = $_[1];
my $char = $_[2];
my @new_array=();
foreach my $i (@array){
foreach my $j (0..length($char)-1){
substr($i,$pos,1,substr($char,$j,1));
push @new_array,$i;
}
}
return(@new_array);
}
primer2multiple('ATGCVCGCDCTNCCTGAB');
__DATA__
R:ag
Y:CT
M:AC
K:GT
S:gc
W:AT
H:atc
B:gtc
V:gac
D:GAT
N:ATgc



http://www.bio-info-trainee.com/926.html
http://www.bio-info-trainee.com/1528.html





上一篇:对snp进行注释并格式化输出-高级难度题目
下一篇:实现模糊匹配-高级难度
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复

使用道具 举报

11

主题

52

帖子

279

积分

中级会员

Rank: 3Rank: 3

积分
279
发表于 2017-1-15 21:29:48 | 显示全部楼层
占楼
这里会有一个优雅的解法
回复 支持 反对

使用道具 举报

1

主题

16

帖子

120

积分

注册会员

Rank: 2

积分
120
发表于 2017-9-13 17:36:40 | 显示全部楼层
[Perl] 纯文本查看 复制代码
### 将DATA 转为HASH ###
while(<DATA>){
    chomp;
    @a=split ":";
    $base{$a[0]}=uc($a[1])
}
### 输入序列生成最终结果 ###
sub seq2result{
    $raw=shift;
    $count=$raw=~tr/RYMKSWHBVDN/RYMKSWHBVDN/;
    for $new(1..$count){
       if ($new==1){
        push @{"tmp$new"},&one2more($raw);
        next;
       }
       $old=$new-1;
       for (@{"tmp$old"}){
        push @{"tmp$new"},&one2more($_);
       }

    }
    return @{"tmp$count"};

}
### 输入序列,将序列中的第一个简并碱基转化返回2条或多条序列 ###
sub one2more{
    my @result;
    $seq=shift;
    if($seq=~/(R|Y|M|K|S|W|H|B|V|D|N)/i){
        my $single_base=$1;
        for (split //,$base{$single_base}){
            my $tem=$seq;
            $tem=~s/$single_base/$_/i;
            push @result,$tem;
        }
    }
    return @result;

}
print "$_\n" for &seq2result(shift);

__DATA__
R:ag
Y:CT
M:AC
K:GT
S:gc
W:AT
H:atc
B:gtc
V:gac
D:GAT
N:ATgc



做批量的话,如果是输入文件,读取文件,加个循环即可,若是命令行,也是对@_进行循环即可。

perl reduction.pl ACRMKGC

ACAAGGC
ACAATGC
ACACGGC
ACACTGC
ACGAGGC
ACGATGC
ACGCGGC
ACGCTGC

有什么问题还请多多指教!
回复 支持 反对

使用道具 举报

0

主题

1

帖子

41

积分

新手上路

Rank: 1

积分
41
发表于 2018-12-12 09:48:43 | 显示全部楼层
本帖最后由 bagadesu 于 2018-12-12 10:03 编辑

[Python] 纯文本查看 复制代码
import re
ambi_base_map = {'R':['A', 'G'],
                 'Y':['C', 'T'],
                 'M':['A', 'C'],
                 'K':['G', 'T'],
                 'S':['C', 'G'],
                 'W':['A', 'T'],
                 'H':['A', 'C', 'T'],
                 'B':['C', 'G', 'T'],
                 'V':['A', 'C', 'G'],
                 'D':['A', 'G', 'T'],
                 'N':['A', 'C', 'G', 'T']}


#基本思路:
def ambiguous_to_unambiguous(ambiguous_seq): #假设输入的序列为'ATMCYN'
    ambiguous_base = re.findall(r'[RYMKSWHBVDN]', ambiguous_seq) #查找出所有简并碱基,存放于一个列表中,这里是['M', 'Y', 'N']。
    out = [ambiguous_seq] #把原始序列存放在一个列表当中。
    for ambi_base in ambiguous_base: #遍历搜寻到的所有简并碱基,每次把原始序列最左边的一个简并碱基转化成非简并碱基,
                                     #把得到的序列存放在变量out当中。如第一次循环后,out应该变成:['ATACYN', 'ATCCYN']。
        result = []
        for i in ambi_base_map[ambi_base]:
            for j in out:
                result.append(j.replace(ambi_base, i, 1))
        out = result
    return out

#测试
ambiguous_to_unambiguous('ATMCYN')

#输出结果:
#['ATACCA', 'ATCCCA', 'ATACTA', 'ATCCTA', 'ATACCC', 'ATCCCC', 'ATACTC', 'ATCCTC', 'ATACCG', 'ATCCCG', 'ATACTG', 'ATCCTG', 'ATACCT', 'ATCCCT', 'ATACTT', 'ATCCTT']
回复 支持 反对

使用道具 举报

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

本版积分规则

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

GMT+8, 2019-9-17 07:15 , Processed in 0.029667 second(s), 25 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.