搜索
查看: 3707|回复: 2

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

[复制链接]

13

主题

66

帖子

276

积分

中级会员

Rank: 3Rank: 3

积分
276
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。
    代码如下(我这里省略了对引物序列的反向互补序列的生成):
  [mw_shl_code=perl,true]
  1 #!/home/huizhen/bin/perl/bin/perl
  2
  3 use strict;
  4 use warnings;
  5 use Getopt:ong;
  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(<R>){
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 readread\n";
87 print "The 5 end trimmedplus\n";
88 print "The 3 end trimmedreverse\n";
[/mw_shl_code]




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

使用道具 举报

13

主题

66

帖子

276

积分

中级会员

Rank: 3Rank: 3

积分
276
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, 2020-4-2 11:47 , Processed in 0.023432 second(s), 25 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.