获取序列信息的数据库有很多,首先介绍从NCBI获取数据:
想从NCBI上获取数据,学会使用Entrez是必不可少的!什么是Entrez?
Entrez 是美国国家生物技术信息中心所提供的在线资源检索器。该资源将GenBank序列与其原始文献出处链接在一起。 Entrez 是由NCBI主持的一个数据库检索系统。
Entrez中核酸数据库为:GenBank,EMBL,DDBJ
蛋白质数据库为:Swiss-Prot,PIR,PFR,PDB
PubMed
基因组和染色体图谱资料
NCBI上是如何整合数据的呢?
LOCUS AF086833 18959 bp cRNA linear VRL 13-FEB-2012
DEFINITION Ebola virus - Mayinga, Zaire, 1976, complete genome.
ACCESSION AF086833
VERSION AF086833.2 GI:10141003
LOCUS:这一行里包括基因座的名字,核酸序列长度,分子的类别,拓扑类型,原核生物的基因拓扑类型都是线性的,最后是更新日期。
DEFINITION:对序列的简短定义。
ACCESSION:就是在搜索条中输入的那个数据库编号,也叫做检索号,每条记录的检索号在数据库中是唯一且不变的。即使数据提交者改变了数据内容,Accession 也不会变。
Version:版本号和 Locus,Accession 长得差不多。版本号的格式是“检索号点上一个数字”。版本号于 1999 年 2 月由三大数据库采纳使用。主要用于识别数据库中一条单一的特定核苷酸序列。在数据库中,如果某条序列发生了改变,即使是单碱基的改变,它的版本号都将增加,而它的 Accession 也就是检索号保持不变。比如,版本号由 U12345.1 变为 U12345.2,而检索号依然是 U12345。版本号后面还有个 GI 号。
GI号:与前面的版本号系统是平行运行的。当一条序列改变后,它将被赋予一个新的 GI 号,同时它的版本号将增加。
Locus 是一个同学的真实姓名,而 Accession 是这个同学的学号。同一个人在不同的学校里会有不同的学号,而名字只有一个。基因也是一样,同一个基因在不同的数据库中会有不同的检索号,而基因的名字只有一个。
如何安装和使用Entrez? 参见:Entrez安装和使用说明
使用示例:
#查看其序列信息(前10行)
efetch -db=nuccore -format=gb -id=AF086833 | head
#导入Genbank format
efetch -db=nuccore -format=gb -id=AF086833 > AF086833.gb
#导入Fasta format
efetch -db=nuccore -format=fasta -id=AF086833 > AF086833.fa
#查看序列的前三个碱基(区域可选)
efetch -db=nuccore -format=gb -id=AF086833 -seq_start=1 -seq_stop=3
#查看序列的反向互补序列
efetch -db=nuccore -format=fasta -id=AF086833 -seq_start=1 -seq_stop=5 -strand=1
efetch -db=nuccore -format=fasta -id=AF086833 -seq_start=1 -seq_stop=5 -strand=2
如何用Entrez搜索?
#在核酸与蛋白质数据库中搜索‘PRJNA257197’
esearch -db nucleotide -query PRJNA257197
esearch -db protein -query PRJNA257197
#结果如下:
#将搜索到的核酸与蛋白质序列全部抓取
esearch -db nucleotide -query PRJNA257197 | efetch -format=fasta >genomes.fa
esearch -db protein -query PRJNA257197 | efetch -format=fasta > proteins.fa
如何从Short Read Archive(SRA)获取数据?
利用SRA Toolkit,使用与安装见以下链接:
#将sra数据转换成fastq格式
fastq-dump
#下载sam格式数据
sam-dump
#产生一个SRR1553607.fastq
fastq-dump SRR1553607
fastq -dump --split-files SRR1553607 # --split-files:将每一个read读到单独的文件中。文件将接收到读取相应的后缀数。针对pairend 测序结果,将原文件转换成两个fastq文件,里面的reads是成对的,文件名其中一个是*_1.fastq,一个是*_2.fastq 得到:
SRR1553607_1.fastq
SRR1553607_2.fastq
fastq-dump --split-files -X 10000 SRR1553607
#提取需要的子集(前10000个读长)
#简要查看数据信息
sra-stat --xml --quick SRR1553607
#查看一个数据的读长统计
sra-stat --xml --statistics SRR4237168
如何批量获取SRA数据?
以 bioproject number PRJNA257197为例:
#首先搜索数据库并下载runinfo并导入到一个csv文本中:
esearch -db sra -query PRJNA257197 | efetch -format runinfo > info.csv
#查看第一列信息:
cat info.csv|cut -d ',' -f 1|head
#将第一列SRR开头的id写入一个文本:
cat info.csv | cut -f 1 -d ',' | grep 'SRR' | head > ids.txt
#获取ids.txt中id的前10000个读长
cat ids.txt | xargs -n 1 fastq-dump -X 10000 --split-files $1
#还可以通过筛选日期来得到想要的序列
cat info.csv | grep '2014-08-19' | cut -f 1 -d ',' | grep SRR | head -10 > ids.txt
下面说说如何查看测序数据信息
#根据发布日期获取测序数据信息:
esearch -db sra -query '"2015/05"[Publication Date]'|efetch -format runinfo>2015-05.cs
#收集SRR开头为编号的数据:
cat sra/*.csv | grep 'SRR' > allruns.csv
下面利用csvcut查看csv文件:
#统计一共测序了多少个碱基:
cat allruns.csv | csvcut -c 5 | datamash sum
1 4.966817753558e+15 #如果报错可删除对应的行
#查看测序物种并统计频次:
cat allruns.csv|csvcut -c 29|sort|uniq -c|sort -nr|head -25
#对使用的测序平台进行统计:
cat allruns.csv|csvcut -c 19|sort|uniq -c|sort -rn|head -25
#查看测序仪型号:
cat allruns.csv|csvcut -c 20|sort|uniq -c|sort -rn|head -25
#查看哪个数据容量最大:
cat allruns.csv|csvcut -c 1,8|datamash -f -t ',' mean 2 sum 2 max 2