zoukankan      html  css  js  c++  java
  • 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;
  • 相关阅读:
    TCP与UDP的差别以及TCP三次握手、四次挥手
    MAC帧格式、IPV4数据报格式、TCP报文格式、UDP数据报格式
    维特比算法(Viterbi)-实例讲解(暴力破解+代码实现)
    对ajax的理解
    get与post两种方式的优缺点
    什么是Ajax和JSON,它们的优缺点
    浅谈一下如何避免用户多次点击造成的多次请求
    ajax是什么?
    同步和异步的区别?
    如何解决跨域问题
  • 原文地址:https://www.cnblogs.com/lyon2014/p/4022489.html
Copyright © 2011-2022 走看看