primer漏配问题解决

在对之前的ITS数据(454数据)做split时,发现有一些reads没有被匹配上,但是barcode能够完全匹配,虽然之后的primer在中间漏了一个碱基,导致后面的碱基全部误匹配,从而导致这条reads没有被匹配上的问题。

终于解决Qiime的问题后,使用 split_libraries.py 做切分,发现同样有这样的问题,Qiime并没有解决漏匹配的问题。

考虑如果用正常方法去做的话,对较小的异常数据需要花费N倍于正常数据的计算资源(包括硬件资源和运行时间),对于这个问题来说是非常不明智的。

由于primer长度在20左右,我的解决办法是,取末端6个碱基,在reads的这6个碱基相应位置左移1-2位,作为漏匹配1-2位的替代处理,这样既解决了漏匹配的问题,而且还能够使原正常的数据匹配速度进一步加快。

 1 #!/usr/bin/perl
 2 use strict;
 3 
 4 my $usage = "usage:
split.pl	mapfile	.fa_file	outprefix
";
 5 die $usage unless @ARGV==3;
 6 
 7 my $mapfile = shift @ARGV;
 8 my $fafile = shift @ARGV;
 9 my $outprefix = shift @ARGV;
10 
11 my %barcode;my $barcode_length = 0;
12 my %primer;my $primer_length = 0;
13 open MAP,$mapfile or die $!;
14 while(<MAP>){
15     chomp;
16     next if /^#/;
17     my @a = split /s+/;
18     $barcode{$a[0]} = $a[1];
19     $barcode_length = length($a[1]) unless $barcode_length;
20     die "barcode length do not match!" unless ($barcode_length == length($a[1]));
21     $primer{$a[0]} = $a[2];
22     $primer_length = length($a[2]) unless $primer_length;
23     die "primer length do not match!" unless ($primer_length == length($a[2]));
24 
25     print "$barcode_length	$primer_length
";
26 }
27 close MAP;
28 
29 my %fa;
30 open FA,$fafile or die $!;
31 $/ = ">";
32 <FA>;
33 while(<FA>){
34     chomp;
35     my @a = split /
/;
36     my $id = shift @a;
37     my $seq = join ("",@a);
38     @a = split (/s+/,$id);
39     $id = shift @a;
40     $fa{$id} = $seq;
41 }
42 $/ = "
";
43 close FA;
44 
45 open OUT,">$outprefix.fna" or die $!;
46 foreach my $id (sort keys %fa){
47     foreach my $sample (sort keys %barcode){
48         my $seq = substr($fa{$id},0,$barcode_length);# print "$seq
";
49         if ($barcode{$sample} eq $seq){
50     #        print "barcode matched
";
51             my $pri0 = substr($fa{$id},$barcode_length+$primer_length-6,6);
52             my $pri1 = substr($fa{$id},$barcode_length+$primer_length-7,6);
53             my $pri2 = substr($fa{$id},$barcode_length+$primer_length-8,6);
54             my $pri = substr($primer{$sample},$primer_length-6,6);
55             if ($pri0 eq $pri){
56                 my $s = substr($fa{$id},$barcode_length+$primer_length,length($fa{$id})-$barcode_length-$primer_length);
57                 print OUT ">$sample	$id
$s
";
58                 last;
59             }
60             if($pri1 eq $pri){
61                 
62                 my $s = substr($fa{$id},$barcode_length+$primer_length-1,length($fa{$id})-$barcode_length-$primer_length+1);
63                 print OUT ">$sample	$id
$s
";
64                 last;
65             }
66             if($pri2 eq $pri){
67                 my $s = substr($fa{$id},$barcode_length+$primer_length-2,length($fa{$id})-$barcode_length-$primer_length+2);
68                 print OUT ">$sample	$id
$s
";
69                 last; 
70             }
71         } 
72     }
73 }
74 close OUT;
原文地址:https://www.cnblogs.com/lyon2014/p/4022489.html