搜索
查看: 2992|回复: 2

去除靶向测序fastq文件中的扩增引物序列

[复制链接]

11

主题

54

帖子

242

积分

中级会员

Rank: 3Rank: 3

积分
242
QQ
发表于 2016-12-25 21:22:51 | 显示全部楼层 |阅读模式
本帖最后由 惠真-市二医院 于 2016-12-25 21:25 编辑

    针对靶向扩增测序数据,这类数据往往是几个基因的基因组外显子区域扩增序列!实验上面就是在建库前做多重PCR扩增,然后对扩增产物构建二代测序文库,所以文库中就会代入扩增引物序列!在下机数据中,在我们比对数据到refseq( 目的基因的扩增序列)时,如果不提前去除扩增引物序列,会导致该引物序列区域的假阴性结果,因为引物序列肯定是没有突变的,但是该区域的样本基因组序列上就有可能发生突变,所以为了排除引物序列的干扰,我们要在比对前去除引物序列!
    引物序列可以在fastq文件时去除,或者在bam文件时去除;该引物序列的存在是有助于序列map到refseq的,所以在bam文件时去除引物序列是理想情况,QIAGEN的网页分析就是在bam文件时,对引物序列做soft clipped处理!只是我还没有找到方法或软件实现这一soft clipped过程,所以我这里是在fastq格式时做的引物去除!
   思路就是:
    1.引物序列只可能出现在read序列的两端(当然扩增靶向测序都会采用PE测序,且测通);
    2.fastq序列都是5端到3端的顺序,厂家给的引物序列也是5端到3端的顺序;
    3.一个read的两端理论上都是引物序列,所以要针对引物序列的正向,反向以及正向互补和反向互补来做trim。
    代码如下(我这里省略了对引物序列的反向互补序列的生成):
  
[Perl] 纯文本查看 复制代码
  1 #!/home/huizhen/bin/perl/bin/perl
  2
  3 use strict;
  4 use warnings;
  5 use Getopt::Long;
  6 use List::Util qw(max min sum);
  7
  8 my ($fq,$primer,$out,$help);
  9 GetOptions(
 10     "fastq=s" => \$fq,
 11     "primer=s" => \$primer,
 12     "out=s" => \$out,
 13     "help|?" => \$help,
 14     );
 15
 16 my $usage=<<END;
 17 Usages:
 18     perl $0 [Options]
 19 Options:
 20     -fastq <fastq.file>
 21     -primer <primer.file>
 22     -out <out.fastq.file>
 23 END
 24 die $usage if (! $fq||! $primer||! $out|| $help);
 25
 26 open PR,"<",$primer or die "cant open:$!";
 27 open FQ,"<",$fq or die "cant open:$!";
 28 open OUT,">",$out or die "cant open:$!";
 29
 30 my (@start,@end);
 31 while(<PR>){
 32     chomp;
 33     my ($plus_5,$plus_3) = (split /\t/,$_)[3,4];
 34     push @start,$plus_5;
 35     push @start,$plus_3;
 36     $plus_5 =~ tr/ATCG/TAGC/;
 37     $plus_3 =~ tr/ATCG/TAGC/;
 38     $plus_5 = reverse $plus_5;
 39     $plus_3 = reverse $plus_3;
 40     push @end,$plus_5;
 41     push @end,$plus_3;
 42 }
 43 close PR;
 44
 45 my %hash;
 46 while(<FQ>){
 47     chomp;
 48     $hash{$.} = $_;
 49 }
 50 close FQ;
 51
 52 my ($plus,$reverse) = (0,0);
 53 foreach my $num(sort {$a <=> $b} keys %hash){
 54     next unless $num%4 == 0;
 55     my $seq = $hash{$num-2};
 56     my $qual = $hash{$num};
 57     foreach my $p(@start){
 58         chomp $p;
 59         if($seq =~ /^$p/){
 60             my $len = length($p);
 61             $seq = substr($seq,$len);
 62             $qual = substr($qual,$len);
 63             $plus ++;
 64             last;
 65         }
 66     }
 67     foreach my $r(@end){
 68         chomp $r;
 69         if($seq =~ /$r$/){
 70             my $len_r = length($r);
 71             my $len_seq = length($seq);
 72             $seq = substr($seq,0,$len_seq - $len_r);
 73             $qual = substr($qual,0,$len_seq - $len_r);
 74             $reverse ++;
 75             last;
 76         }
 77     }
 78     print OUT "$hash{$num-3}\n";
 79     print OUT "$seq\n";
 80     print OUT "+\n";
 81     print OUT "$qual\n";
 82 }
 83 close OUT;
 84
 85 my $read = (max keys %hash)/4;
 86 print "The number of read:$read\n";
 87 print "The 5 end trimmed:$plus\n";
 88 print "The 3 end trimmed:$reverse\n";





上一篇:根据PubMed ID批量提取文献摘要
下一篇:基因融合结果展示
苛求远离完美
回复

使用道具 举报

11

主题

54

帖子

242

积分

中级会员

Rank: 3Rank: 3

积分
242
QQ
 楼主| 发表于 2016-12-25 21:27:59 | 显示全部楼层
    谁要是有 对bam文件trim引物序列的思路或者软件,请告诉我,感激啊!针对bam的trim就是觉得对CIGAR做了soft clipped标记,感觉这个很复杂!
苛求远离完美
回复 支持 反对

使用道具 举报

0

主题

2

帖子

55

积分

注册会员

Rank: 2

积分
55
QQ
发表于 2017-8-29 15:25:37 | 显示全部楼层
python中的pysam可以处理bam文件
回复 支持 反对

使用道具 举报

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

本版积分规则

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

GMT+8, 2019-9-21 09:33 , Processed in 0.032336 second(s), 28 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.