注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

wangyufeng的博客

祝愿BB 健康开心快乐每一天

 
 
 

日志

 
 

如何向NCBI提交测序数据(转)  

2013-08-01 19:40:07|  分类: 生物信息分析 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |

一、submit reads

1、prepare  .fq data

准备.fa格式的数据,这个其实平常用的就是,不是的话估计还是比较麻烦的。其格式如下:

@FC61K87AAXX:6:1:1035:17969#0/2
NNNNNNACAGCACAGNNNNNNNGNNNNNNAACAGGTGNNNNNAA
+
BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB
@FC61K87AAXX:6:1:1035:2669#0/2
NNNNNNCAAATTCCANNNNNNNGNNNNNNNAACGTGANNNNNGN
+
BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB
2、NCBI网页注册和准备

进入主页(http://www.ncbi.nlm.nih.gov/guide/)

1) Submissions:

2) Sequence Read Archive (SRA)

3) Submit to SRA

4) NCBI Primary Data Archive Submitters

5)sign in

Username:bgi

Password:*******

在这儿最好不记住密码,有时候合作方会先注册,你再提交后有时会出不来。

6)create new submission

7)后面的就要自己创建任务。

注意:在申请的过程可能还需要申请Project ID。这首先要确认自己的项目是否已经注册过了,一个项目一个ID,申请多了的话不能自己注销,和NCBI的人员交流需要很长时间,而且只能用邮件。本身提交数据就是一个拉锯战。

3、上传数据。

到此为止,准备就做完了,之后就要上传数据了。自己上传不了,需要用ftp上传,在华大就是向3811发任务。

这儿有具体的发任务格式,按要求进行。

猪蛔虫实例:

http://trace.ncbi.nlm.nih.gov/Traces/sra_sub/sub.cgi?subid=69656&from=list&action=show:submission

 

 


二、submit assembly to NCBI

1、prepare data

首先要具有fasta格式数据(NO .gz),这是处理的基础,具体格式如下:

  1. >Scaffold633  
  2.   
  3. TCATTTCTCCACTCTCGATGAACAAATCTGGAGGGATTTTTTTTCATTCC  
  4.   
  5. ACTCAATAGGTTGTCTATAAAGGTGTGATTCGTGGAACTTCTTCACACAG  
  6.   
  7. CAGCTAGTCTATATAATACAGAAGATCG  
  8.   
  9. >Scaffold553  
  10.   
  11. AAAAAATTTTTTTTTTAAACTATCATCTCATGGATCAGCAGCAATTCTGA  
  12.   
  13. GTGTAACGTCTTCATTAAATGCGTATATAAATTTGCATAAAGATATGCGA  
  14.   
  15. CCAATATTGAGCCTGGAAATATATGCGCAGAGTGCAAAATTGTGTTTTTT  
  16.   
  17. GATCGGTTAATTAAAGG  
  18.   
  19. >Scaffold641  
  20.   
  21. GTTTCCCAGTAGGTCTCTCCCGCTACGGCGTCCGCACGAACGCGATCTGC  
  22.   
  23. CCTCGTGCCCGCACCGCCATGACGGCAGAAGCCTTCGGCGAGAACAACAC  
  24.   
  25. CGGCGTCGTCGGCCTCGATCCGCTTGCACCCGAGCGCGTCGCGACCCTGG  
  26.   
  27. TCAGCTACCTCGCATCCCCCGATTCCGACGAGATCAACGGACAGGTCTTC  
  28.   
  29. GTCGTCTACGGCAAGATGGTGGCGTTGATGGAAGCACCCAAGGTCGAGAA  
  30.   
  31. CCGTTTCGACGCAGCCGGATCCGCGTTCACCGTCGAAGAACTCGGTGGCC  
  32.   
  33. AGCTCTCGTCTTACTTCTCCGGCCGTGGGCCGTACGAGACCTACTGGGAA  
  34.   
  35. AC  

2、处理数据

分为几步:

(1)生成.greater, short.list和ZERO_BASE_COUNT文件

perl  ../ scaf_filter_2k.pl  Ascaris_suum.scaf.fa

scaf_filter_2k.pl代码

  1. #!/usr/bin/perl  
  2.   
  3. use strict;  
  4.   
  5. use warnings;  
  6.   
  7.    
  8.   
  9. my $file=shift;  
  10.   
  11. #my $cutoff=shift;  
  12.   
  13. my $outfile="short.list";  
  14.   
  15. my $outfile2="$file.greater";  
  16.   
  17. my $outfile3="ZERO_BASE_COUNT";  
  18.   
  19.    
  20.   
  21. open IN,"< $file" or die $!;  
  22.   
  23. open OUT,"> $outfile" or die $!;  
  24.   
  25. open OUT1,"> $outfile2" or die $!;  
  26.   
  27. open OUT2, "> $outfile3" or die $!;  
  28.   
  29. $/='>';<IN>;$/="\n";  
  30.   
  31. while(<IN>){  
  32.   
  33.         chomp;  
  34.   
  35.         my $id=$1 if(/^(\S+)/);  
  36.   
  37.         $/='>';  
  38.   
  39.         my $seq=<IN>;  
  40.   
  41.         chomp($seq);  
  42.   
  43.         $/="\n";  
  44.   
  45.         $seq=~s/\s//g;  
  46.   
  47.         my $len=length($seq);  
  48.   
  49.         if ($len < 200){  
  50.   
  51.                 print  OUT "$id\t$len\n";  
  52.   
  53.                 next;  
  54.   
  55.         }  
  56.   
  57.         else{  
  58.   
  59.                 my $a=$seq=~tr/aA/aA/;  
  60.   
  61.                 my $t=$seq=~tr/tT/tT/;  
  62.   
  63.                 my $c=$seq=~tr/cC/cC/;  
  64.   
  65.                 my $g=$seq=~tr/gG/gG/;  
  66.   
  67.                 if($a==0 || $c==0 || $g==0 || $t==0){  
  68.   
  69.                         print OUT2 "$id\t$len\n";  
  70.   
  71.                         next;  
  72.   
  73.                 }  
  74.   
  75.                 print OUT1 ">$id\n$seq\n";  
  76.   
  77.         }  
  78.   
  79. }  
  80.   
  81. close IN;  
  82.   
  83. close OUT;  
  84.   
  85. close OUT1;  
  86.   
  87. close OUT2;  

(2)生成.Nchange文件。

perl ../ info_N_plus.pl  Ascaris_suum.scaf.fa.greater  > Ascaris_suum.scaf.fa.greater.Nchange

info_N_plus.pl代码:

  1. #!/usr/bin/perl -w  
  2.   
  3. use strict;  
  4.   
  5. #use Getopt::Long;  
  6.   
  7. sub usage{  
  8.   
  9.   print STDERR <<USAGE;  
  10.   
  11.       ############################################  
  12.   
  13.             Version 1.1 by Wing-L 2011.07.15  
  14.   
  15.    
  16.   
  17.       usage: $0 <sequence.fa> <len> >STDOUT  
  18.   
  19.    
  20.   
  21.       ############################################  
  22.   
  23. USAGE  
  24.   
  25. exit;  
  26.   
  27. }  
  28.   
  29. &usage if(@ARGV <1);  
  30.   
  31. my ($fa,$len)=@ARGV;  
  32.   
  33. $len||=10;  
  34.   
  35. open IN,"<$fa" or die("$!\n");  
  36.   
  37. $/='>';<IN>;$/="\n";  
  38.   
  39. while(my $line=<IN>){  
  40.   
  41.     my @block;  
  42.   
  43.   $line=~/^\S+/;  
  44.   
  45.   my $tag={1};  
  46.   
  47.   $/='>';  
  48.   
  49.   my $seq=<IN>;  
  50.   
  51.   chomp $seq;  
  52.   
  53.   $/="\n";  
  54.   
  55.   $seq=~s/\s//g;  
  56.   
  57.   my $chr_length=length $seq;  
  58.   
  59.   while($seq=~/[^N]N{1,9}[^N]/g){  
  60.   
  61.     substr ($seq,$-[0]+1,$len)=~s/\S/N/g;  
  62.   
  63.   }  
  64.   
  65.   while($seq=~/N([^N]{1,49})N/g){  
  66.   
  67.       my $tlen=length $1;  
  68.   
  69.       substr ($seq,$-[0]+1,$tlen)=~s/\S/N/g;  
  70.   
  71.   }  
  72.   
  73.   if($seq=~/^N?[^N]{0,49}N+/){  
  74.   
  75.       print STDERR "$tag\t1\t{1}[0]\n";  
  76.   
  77.     substr($seq,0,{1}[0]-$-[0])='';  
  78.   
  79.   }  
  80.   
  81.   if($seq=~/N+[^N]{0,49}N{0,}$/){  
  82.   
  83.     print STDERR "$tag\t$-[0]\t$chr_length\n";  
  84.   
  85.     substr($seq,$-[0],$chr_length-$-[0])='';  
  86.   
  87.   }  
  88.   
  89.   print ">$tag\n$seq\n";  
  90.   
  91. }  
  92.   
  93. close IN;  
  94.   
  95. #open IN,"<" or die("$!\n");  
  96.   
  97. #while(my $line=<IN>){}  
  98.   
  99. #foreach my $e (){}  
  100.   
  101. #(split /\s+/,$line)[0]  
  102.   
  103. #open OUT,">" or die("$!\n");  
  104.   
  105. ###############  sub  ###############  

(3)生成分割文件

perl  ../get_scaftig.pl  Ascaris_suum.scaf.fa.greater.Nchange  > Ascaris_suum.scaf.fa.greater.Nchange.split

get_scaftig.pl代码:

  1. #!/usr/bin/perl -w  
  2.   
  3. #  
  4.   
  5. #Author: Ruan Jue <ruanjue@genomics.org.cn>  
  6.   
  7. #  
  8.   
  9. use warnings;  
  10.   
  11. use strict;  
  12.   
  13.    
  14.   
  15. my $min_length = 0;  
  16.   
  17. my $name = '';  
  18.   
  19. my $seq = '';  
  20.   
  21.    
  22.   
  23. while(<>){  
  24.   
  25.    if(/^>(\S+)/){  
  26.   
  27.       &print_scafftig($name, $seq) if($seq);  
  28.   
  29.       $name = $1;  
  30.   
  31.       $seq  = '';  
  32.   
  33.    } else {  
  34.   
  35.       chomp;  
  36.   
  37.       $seq .<pre class="html" name="code">{1}</pre><br>  
  38. ; }}&print_scafftig($name, $seq) if($seq); 1; sub print_scafftig { my $name = shift; my $seq = shift; my $temp = $seq; my $id = 1; my $flag = 0; my $pos = 1; while($seq=~/([ATGCatgc]+)/g){ my $s = $1; if($flag==1){ if($temp=~/([ATGCatgc]+[Nn]+)/g){ my $g =  
  39.  $1; $pos+=length($g); } } else{$flag=1;} next if(length($s) < $min_length); my $length=length($s); my $end=$pos+$length-1;# print ">$name\_$id start=$pos length=".length($s)."\n"; print ">$name\_$id\t$pos\t$end\t$length\n"; while($s=~/(.{1,60})/g){ print "$1\n";  
  40.  } $id++;}}  
  41. <pre></pre>  
  42. <p align="left">(4)生成.agp文件</p>  
  43. <p align="left">perl ../AGP.pl Ascaris_suum.scaf.fa.greater.Nchange.split > Ascaris_suum.scaf.fa.agp</p>  
  44. <p align="left">AGP.pl 代码:</p>  
  45. <pre class="html" name="code">#!/usr/bin/perl  
  46.   
  47. use strict;  
  48.   
  49.    
  50.   
  51. my $fasta=shift;  
  52.   
  53.    
  54.   
  55. my %Scaf;  
  56.   
  57. open IN,$fasta or die "$!";  
  58.   
  59. while(<IN>){  
  60.   
  61.         if (/^>(\S+)(_\d+)\s+(\d+)\s+(\d+)\s+(\d+)/){  
  62.   
  63.                 my $scaf=$1;  
  64.   
  65.                 my $contig=$1.$2;  
  66.   
  67.                 my $start=$3;  
  68.   
  69.                 my $end=$4;  
  70.   
  71.                 my $len=$5;  
  72.   
  73.                 push @{$Scaf{$scaf}},[$contig,$start,$end,$len];  
  74.   
  75.         }  
  76.   
  77. }  
  78.   
  79. close IN;  
  80.   
  81.    
  82.   
  83. foreach my $id ( sort keys %Scaf ){  
  84.   
  85.         @{$Scaf{$id}} = sort {$a->[1]<=>$b->[1]}  @{$Scaf{$id}};  
  86.   
  87. }  
  88.   
  89.    
  90.   
  91. foreach my $id ( sort { $Scaf{$b}[-1][2] <=> $Scaf{$a}[-1][2] } keys %Scaf ){  
  92.   
  93.         my $p=$Scaf{$id};  
  94.   
  95.         my $idx=1;  
  96.   
  97.         for(my $i=0;$i<@$p;$i++){  
  98.   
  99.                 if ($i>0){  
  100.   
  101.                         print join("\t",$id,$p->[$i-1][2]+1,$p->[$i][1]-1,$idx,'N',($p->[$i][1]-1)-($p->[$i-1][2]+1)+1,'fragment','yes',"\t")."\n";  
  102.   
  103.                         $idx++;  
  104.   
  105.                 }  
  106.   
  107.                 print join("\t",$id,$p->[$i][1],$p->[$i][2],$idx,'W',$p->[$i][0],1,$p->[$i][3],'+')."\n";  
  108.   
  109.                 $idx++;  
  110.   
  111.         }  
  112.   
  113. }  
  114.   
  115. </pre>  
  116. <p align="left">至此为止,第一步准备数据的过程已经完成。下一步就是向3811提交任务。</p>  
  117. <p align="left">这个过程需要的东西:</p>  
  118. <p align="left">l  Ascaris_suum.scaf.fa.greater.Nchange.gz</p>  
  119. <p align="left">l  Ascaris_suum.agp.gz</p>  
  120. <p align="left">和NCBI沟通的文件,主要是确定的Project ID的准确性。</p>  
  121. <p align="left">未完待续.........</p>  
  122. <p> </p>  
  123. <pre></pre>  
以上内容转自:http://blog.csdn.net/bangemantou/article/details/7294911
  评论这张
 
阅读(2282)| 评论(0)
推荐 转载

历史上的今天

在LOFTER的更多文章

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2017